# Evaluation of Ecological Water Consumption in Yanhe River Basin Based on Big Data

Authors: Ting Guo; Huiying Yu
Journal: Computational Intelligence and Neuroscience
 (2021)

---

## Abstract

Starting from the main eco-environmental problems faced by water environment, taking Yanhe River Basin as an example, this paper discusses the theoretical connotation and evaluation calculation method of eco-environmental water consumption. In order to study the eco-environmental water consumption of Yanhe River Basin, a runoff driving factor mining method based on big data analysis is established in this paper. Aiming at the problem that the statistical law and genetic law of runoff change frequently in changing environment, the mining technology method of runoff key driving factors is proposed by combining traditional methods with big data analysis. The characteristic factors that have no significant impact on runoff change are removed, the implicit characteristic factors affecting runoff change are extracted, the driving relationship of hydrological, meteorological, and vegetation characteristic factors on ecological water consumption change is identified, and the key driving factors of ecological water consumption change are extracted, which lays a data foundation for ecological water consumption prediction based on machine learning. The economic water consumption based on eco-environmental water consumption in Yanhe River Basin in the future is predicted (including water demand in three aspects of industry, agriculture, and life); that is, the prediction is to meet the economic water demand on the basis of ensuring that the water consumption of ecological environment will not be occupied, which can effectively ensure the improvement of ecological environment function in Yanhe River Basin and is conducive to the sustainable utilization of water resources in Yanhe River Basin. The research is only based on a small watershed such as Yanhe River Basin, and the purpose of the research is to provide a reference for ecological environment protection and sustainable utilization of water resources in the Loess Plateau, even in the arid, semiarid, and semihumid areas of North China.


---

## Body

## 1. Introduction

Facing the resources and environmental problems in China in the 21st century, especially the water problem, is the most serious. Less water and waste of water coexist. Water abundance coexists with ecological imbalance. Water dirty and water management coexist. In fact, water resources, water disasters, and water environment are interrelated and mutually transformed. As the core element of ecology and environment, the utilization of water is a positive benefit of water resources [1, 2]: the change of its surplus and shortage causes negative effects of floods and droughts; the pollution caused by its quality evolution is a negative impact on the environment. Ecological water consumption is the most active and sensitive central factor in Yanhe River Basin, and around this center, there are many water consumption contradictions [3]. At present, the contradiction between industrial and agricultural production water and ecological environment water is the most urgent and has the greatest impact on the sustainable development of local economy and society [4]. Because of the lag effect of eco-environmental problems, this problem has not attracted people’s attention before, especially the lag of theoretical research on eco-environmental water consumption [5, 6], which directly affects the solution of this problem. This paper takes Yanhe River Basin as the research object, studies the ecological environmental water consumption of Yanhe River Basin, and evaluates the ecological water consumption of Yanhe River Basin based on big data analysis [7, 8].


## 2. Research on Runoff Driving Factor Mining Based on Big Data Analysis

Ecological water consumption prediction based on machine learning needs data-driven, which is a typical supervised learning problem [9]. It needs to rely on a large number of data samples to mine and establish the implicit relationship between driving factors and runoff. Therefore, screening runoff driving factors and constructing a runoff prediction sample set are the basis of constructing a runoff prediction model in response to changes based on machine learning technology. Based on big data analysis, this paper studies the driving factor mining for runoff prediction in Yanhe River Basin. Big data analysis is aimed at a large number of structured, semistructured, or unstructured data sets from different sources [7]. The advanced analysis technology used can mine the laws within or between data.


### 2.1. Pretreatment of Characteristic Factors of Runoff Prediction

When extracting characteristic data from natural process data (including climate change process and underlying surface evolution process) and human activity process data, there are several problems in runoff prediction modeling [10]:(1)
Some data are missing.(2)
Feature sequences are characterized by complex nonlinearity and highly irregular and multiscale variation and contain a large amount of hidden information, which is difficult to be identified.(3)
There is a linear correlation between features. Therefore, it is necessary to preprocess the characteristic data in order to extract the key driving factors of runoff.On the basis of feature extraction [11], feature data preprocessing first deals with missing values to meet the modeling requirements. Then, based on complex nonlinear feature decomposition and image data dimension reduction, lumped transformation is carried out to obtain lumped feature factor information of the whole basin. Finally, the characteristic factor set is cleaned to exclude the characteristic factors that have no significant influence on runoff change. By decomposing complex nonlinear features, the hidden features in the original time series can be obtained. Lumping transformation is applied to the same feature of different longitude and latitude in the catchment area controlled by the target section, which can evaluate the overall change of the feature in the catchment area. Feature cleaning eliminates the driving factors that have no significant influence on runoff and can reduce the risk of overfitting of the model [12].The methods of lumping transformation include the arithmetic average method [13], Thiessen polygon method, and isoline method [14, 15], whose process and principle are relatively simple, so we will not repeat them here. The following studies are carried out on feature extraction and the preprocessing of missing values, complex nonlinear feature factors, and linear correlation feature factors in the obtained feature sequence.


#### 2.1.1. Feature Factor Extraction

According to the theory of runoff generation and confluence [16], the precipitation, evaporation, and environmental characteristics of the underlying surface in the upstream catchment area controlled by a certain section of a river have a direct impact on the runoff formation process of the section. In addition, direct human activities, such as water intake, water use, water consumption, drainage, and water transfer, can change the temporal and spatial distribution of water resources, thus affecting the runoff formation of target sections. Therefore, according to the information of catchment area, distribution of hydrometeorological stations, and distribution of water intakes and outlets, based on the data sets related to climate change, underlying surface environmental evolution and human activities, the hydrometeorological time series, and underlying surface environmental evolution time series and direct human activity time series driving runoff formation are extracted to extract runoff driving factors and construct runoff prediction sample set.


#### 2.1.2. Processing of Missing Value of Characteristic Factor

Interpolation schematic diagram of piecewise cubic spline method for missing daily flow is shown in Figure1.Figure 1
Interpolation schematic diagram of piecewise cubic spline method for missing daily flow.There are generally two ways to deal with missing values [17]:(1)
When the continuous deletion period of the feature sequence is too long or the number of missing values reaches more than 50% of the total period length, the feature sequence is directly removed.(2)
When there are scattered missing values, an interpolation method is used to make up. Commonly used time series interpolation methods include mean interpolation, correlation analysis (unary linear regression, multiple linear regression, and iterative regression), and spline interpolation. Under the disturbance of changing environment, the statistical characteristics of hydrometeorological time series are constantly changing, so the mean interpolation method is no longer applicable. When the correlation between interpolation sequence and other reference sequences is small, the effect of correlation analysis is poor. The runoff data (see Table1) and meteorological data (see Table 2) used in this study are partially missing or missing, so it is necessary to interpolate the missing values.Table 1
Basic information of hydrological stations in the study basin.
WatershedWatershed areaAverage annual runoffAverage sediment amount for many yearsYanhe River Basin7725 square kilometers2.93 billion cubic meters244∼311 kg/m3Table 2
Basic information of meteorological observation variables.
Variable nameVariable abbreviationUnitAverage station pressureAVG_PShPaDaily maximum station pressureMAX_PShPaDaily minimum station pressureMIN_PShPaAverage temperatureAVG_T°CDaily maximum temperatureMAX_T°CDaily minimum temperatureMIN_T°CRelative humidityRHU%Minimum relative humidityMIN_RHU%Precipitation at 20-8 o’clockP208mmPrecipitation at 8-20 o’clockP820mmPrecipitation at 20-20 o’clockP2020mmSmall evaporation capacityS_EVPmmLarge evaporation capacityL_EVPmmAverage wind speedAVG_Wm/sMaximum wind speedMAX_Wm/sAverage surface temperatureAVG_ST°CMaximum daily surface temperatureMAX_ST°CDaily minimum surface temperatureMIN_ST°C


#### 2.1.3. Complex Nonlinear Eigenfactor Decomposition

(1) Commonly Used Signal Decomposition Algorithms. Decomposition and noise reduction of complex nonlinear, highly nonstationary and multiscale variation feature series, extraction of hidden features in the original series, and construction of runoff prediction model can significantly improve the accuracy of runoff prediction. Commonly used signal decomposition algorithms include Empirical Mode Decomposition (EMD), Wavelet Transform (WT), and Variational Mode Decomposition (VMD) [18]. The following is a brief description of how to use various decomposition algorithms to decompose features and extract hidden features.(2) Empirical Mode Decomposition. EMD is an adaptive time-frequency signal processing algorithm proposed by Huang E in 1998 [19]. A series of applications have proved that EMD is especially suitable for nonlinear and nonstationary signal processing with noise. EMD decomposes the original signal x (t) into several Intrinsic Mode Functions (IMF) and trend terms, (1)xt=∑j=1ncj+rn,where cj is IMF and rn is trend term. IMF is a random oscillation function with different amplitudes and frequencies and meets the following two characteristics: (1) the number of extreme points must be equal to the number of zero crossings, or the maximum difference is 1; (2) the average of the upper envelope (defined by the maximum) and the lower envelope (defined by the minimum) is 0. The EMD calculation process is as follows:(1)
Identify the local maxima and minima in the original sequencex(t), connect the local maxima points by cubic spline function, and obtain the upper and lower envelopes xmax(t) and xmin(t).(2)
Calculate the differencec1 (t) between the original sequence x(t) and the mean sequence m(t) of the upper and lower envelopes(2)mt=xmaxt+xmint2,c1t=xt−mt.(3)
Replace the original sequence withc1(t), and repeat steps 1-2 until the envelope is symmetric, the mean value is zero (i.e., the above two characteristics of IMF are satisfied), the residual residuals are monotonous, and IMF cannot be separated from the residual sequence.The disadvantage of the EMD algorithm is that modal aliasing often occurs, which leads to incomplete IMF separation; that is, multiple IMF components contain repeated information and have a high correlation among components. In addition, because EMD is adaptive, the decomposition level is determined by the algorithm according to the data characteristics. For time series decomposition, the decomposition level is also variable, so it cannot be used for hydrometeorological time series prediction.(3) Wavelet Transform. WT is a mathematical operation that can perform convolution operation on time series or signals in time domain and frequency domain at the same time. Discrete wavelet transform (DWT) is commonly used in hydrology [20]. Given the original signal x(t), DWT can be defined as(3)DWTj,k=2−j/2∫−∞+∞xtψ∗2−jt−kdt,where k is the position index, j is the decomposition level, and ψ∗ is the wavelet function. DWT operates a high-pass filter and a low-pass filter to decompose the original time series into detail subseries representing high-frequency components and approximation subseries representing trends (low-frequency). The specific calculation process is as follows:(1)
The original time seriesx(t) and the decomposition layer j in the wavelet function ψ∗ are given.(2)
DWT decomposesx(t) into a high-frequency detail component (D1(t)) and a low-frequency approximation component (A1(t)) at the first decomposition level, and(4)xt=D1t+A1t.(3)
At the second decomposition level, DWT continues to decomposeA1(t) into a high-frequency detail component (D2(t)) and a low-frequency approximation component (A2(t)), and(5)At=D2t+A2t.(4)
Similar to Step 3, at each subsequent decomposition level, the low-frequency approximate component obtained at the previous decomposition level is decomposed into the high-frequency detail component and the low-frequency approximate component until a preset decomposition level is reached.The original time series is expressed as the sum of approximate components and all detail components; that is,(6)xt=D1t+⋯+Djt+Ajt.The disadvantage of DWT is that there is no uniform standard and theory to determine the wavelet function and decomposition level, and a large number of experiments are needed.(4) Variational Mode Decomposition. VMD is a signal processing algorithm proposed in 2014, which can decompose the original signal into K IMF at one time [21]. The basic idea of VMD is to construct a variational problem and solve the variational problem to obtain decomposed subsignals. Given the input signal f(t), the variational problem can be defined as (7)minukwk∑k∂tδt+jπt∗ukte−jωkt22s.t.∑kukt=ft,where {uk} and {wk} represent the set of modes and their corresponding center frequencies. T represents time, j2 = −1 represents the square of imaginary units, ∗ represents the convolution operator, and σ represents the Dirac trigonometric function.For VMD, the alternating direction multiplier method is used to solve the formula, and the updated formulas of modaluk(w), center frequency wk, and Lagrangian multiplier λ in frequency domain are obtained as follows:(8)u^kn+1ω=f^ω−∑i<ku^kn+1ω+λ^nω/21+2aω−ωk2,ω^kn+1=∫0∞ωu^kn+1ω2dω∫0∞u^kn+1ω2dω,λ^n+1ω=λ^nω+τf^ω−∑ku^kn+1ω,where n is the iterative counter, τ is the noise tolerance, and u^kn+1ω, f^ω, and λ^nω are the Fourier transforms of ukn+1ω, uω, and unω, respectively. Variables u^kn+1ω, f^ω, and λ^nω are continuously updated until they are less than the convergence error ε:(9)∑kμkn+1−μkn22μkn22<ε.The decomposition effect of VMD is affected byK. A small K value cannot extract IMF effectively from the original signal, and a large K value may lead to redundancy of IMF information. A smaller alpha value may lead to larger bandwidth, redundant information expression, and increased additional noise. On the contrary, a larger alpha value may lead to smaller bandwidth and loss of effective information. However, after a large number of experiments, it is determined that when VMD is used to decompose runoff time series, when α = 2000, τ = 0, and ε = 1 × 10−9, it can ensure a good denoising effect and effective IMF separation. In this study, K is optimized by observing whether there is center frequency aliasing, that is, starting from K = 2 (K value increases by 1 every time), and stopping optimization until the last VMD component shows center frequency aliasing.


### 2.2. Characteristic Factor Cleaning

The purpose of feature factor cleaning is (1) to exclude the feature factors that contribute 0 to the prediction target and (2) to reconstruct the relevant feature factors. Feature cleaning can effectively remove the insignificant features, reduce the burden for the subsequent quantitative calculation of driving factors, and effectively avoid the risk of overfitting the machine learning model of runoff prediction. There are two ways to realize feature cleaning [22]: (1) directly remove the features with a variance of 0, analyze the linear correlation or multicollinearity between the remaining features, and reconstruct the related features; (2) principal component analysis or cluster analysis is directly used to reduce the dimension of features and realize feature cleaning. The variance is 0, which means that the value of the feature does not change with time, so it has no influence on the prediction target. The smaller the characteristic variance, the smaller the influence on the prediction target and the smaller the contribution. Figures 2 and 3, respectively, give the variance information of meteorological characteristics and ERA5L characteristics in Yanhe River Basin. It can be seen that although there is no characteristic with a variance of 0, there are a large number of small variance characteristics, which need further cleaning.Figure 2
Variance of meteorological characteristics from the source of the Yellow River Basin to Longyangxia.Figure 3
ERA5L characteristic variance from the source of the Yellow River Basin to Longyangxia.VIF first establishes a regression equation between one feature and other features and then calculates based on the following formula:(10)VIF=11−R2.The original feature set is reorganized into a main component set. Principal components are defined as(11)Zi=ai1∗X1+ai2∗X2+⋯+aiP∗Xp,where Zi represents principal components, ai1, ai2, …, aiP represent feature vectors, X1, X2, …, Xp represent original features, and p represents the number of features obtained from Zi.


### 2.3. Extraction of Driving Factors Based on Genetic Contribution Analysis

Under the changing environment, the factors driving runoff change are constantly changing, so it is necessary to analyze the contribution of driving factors to runoff and determine the dominant factors. In addition, containing a large number of driving factors in the sample set will increase the risk of model overfitting, reduce the generalization prediction ability of the model, and also increase the consumption of modeling resources and time. Therefore, it is necessary to screen key driving factors and exclude driving factors that have less impact on runoff. In this paper, the key driving factors are screened based on the contribution of driving factors to the prediction target. The following is a study on how to quantitatively calculate the contribution of driving factors.Mutual information is a measure of the interdependence between two random variables, which is closely related to the origin of random variables. For two joint discrete random variablesX and Y, their mutual information is defined as(12)IX;Y=∑y∈Y∑x∈XpX,Yx,ylogpX,Yx,ypXxPYy,where P(X,Y)(x, y) is the joint distribution function of X and Y and PX(x) and PY(y) are the edge distributions of X and Y, respectively.


## 2.1. Pretreatment of Characteristic Factors of Runoff Prediction

When extracting characteristic data from natural process data (including climate change process and underlying surface evolution process) and human activity process data, there are several problems in runoff prediction modeling [10]:(1)
Some data are missing.(2)
Feature sequences are characterized by complex nonlinearity and highly irregular and multiscale variation and contain a large amount of hidden information, which is difficult to be identified.(3)
There is a linear correlation between features. Therefore, it is necessary to preprocess the characteristic data in order to extract the key driving factors of runoff.On the basis of feature extraction [11], feature data preprocessing first deals with missing values to meet the modeling requirements. Then, based on complex nonlinear feature decomposition and image data dimension reduction, lumped transformation is carried out to obtain lumped feature factor information of the whole basin. Finally, the characteristic factor set is cleaned to exclude the characteristic factors that have no significant influence on runoff change. By decomposing complex nonlinear features, the hidden features in the original time series can be obtained. Lumping transformation is applied to the same feature of different longitude and latitude in the catchment area controlled by the target section, which can evaluate the overall change of the feature in the catchment area. Feature cleaning eliminates the driving factors that have no significant influence on runoff and can reduce the risk of overfitting of the model [12].The methods of lumping transformation include the arithmetic average method [13], Thiessen polygon method, and isoline method [14, 15], whose process and principle are relatively simple, so we will not repeat them here. The following studies are carried out on feature extraction and the preprocessing of missing values, complex nonlinear feature factors, and linear correlation feature factors in the obtained feature sequence.


### 2.1.1. Feature Factor Extraction

According to the theory of runoff generation and confluence [16], the precipitation, evaporation, and environmental characteristics of the underlying surface in the upstream catchment area controlled by a certain section of a river have a direct impact on the runoff formation process of the section. In addition, direct human activities, such as water intake, water use, water consumption, drainage, and water transfer, can change the temporal and spatial distribution of water resources, thus affecting the runoff formation of target sections. Therefore, according to the information of catchment area, distribution of hydrometeorological stations, and distribution of water intakes and outlets, based on the data sets related to climate change, underlying surface environmental evolution and human activities, the hydrometeorological time series, and underlying surface environmental evolution time series and direct human activity time series driving runoff formation are extracted to extract runoff driving factors and construct runoff prediction sample set.


### 2.1.2. Processing of Missing Value of Characteristic Factor

Interpolation schematic diagram of piecewise cubic spline method for missing daily flow is shown in Figure1.Figure 1
Interpolation schematic diagram of piecewise cubic spline method for missing daily flow.There are generally two ways to deal with missing values [17]:(1)
When the continuous deletion period of the feature sequence is too long or the number of missing values reaches more than 50% of the total period length, the feature sequence is directly removed.(2)
When there are scattered missing values, an interpolation method is used to make up. Commonly used time series interpolation methods include mean interpolation, correlation analysis (unary linear regression, multiple linear regression, and iterative regression), and spline interpolation. Under the disturbance of changing environment, the statistical characteristics of hydrometeorological time series are constantly changing, so the mean interpolation method is no longer applicable. When the correlation between interpolation sequence and other reference sequences is small, the effect of correlation analysis is poor. The runoff data (see Table1) and meteorological data (see Table 2) used in this study are partially missing or missing, so it is necessary to interpolate the missing values.Table 1
Basic information of hydrological stations in the study basin.
WatershedWatershed areaAverage annual runoffAverage sediment amount for many yearsYanhe River Basin7725 square kilometers2.93 billion cubic meters244∼311 kg/m3Table 2
Basic information of meteorological observation variables.
Variable nameVariable abbreviationUnitAverage station pressureAVG_PShPaDaily maximum station pressureMAX_PShPaDaily minimum station pressureMIN_PShPaAverage temperatureAVG_T°CDaily maximum temperatureMAX_T°CDaily minimum temperatureMIN_T°CRelative humidityRHU%Minimum relative humidityMIN_RHU%Precipitation at 20-8 o’clockP208mmPrecipitation at 8-20 o’clockP820mmPrecipitation at 20-20 o’clockP2020mmSmall evaporation capacityS_EVPmmLarge evaporation capacityL_EVPmmAverage wind speedAVG_Wm/sMaximum wind speedMAX_Wm/sAverage surface temperatureAVG_ST°CMaximum daily surface temperatureMAX_ST°CDaily minimum surface temperatureMIN_ST°C


### 2.1.3. Complex Nonlinear Eigenfactor Decomposition

(1) Commonly Used Signal Decomposition Algorithms. Decomposition and noise reduction of complex nonlinear, highly nonstationary and multiscale variation feature series, extraction of hidden features in the original series, and construction of runoff prediction model can significantly improve the accuracy of runoff prediction. Commonly used signal decomposition algorithms include Empirical Mode Decomposition (EMD), Wavelet Transform (WT), and Variational Mode Decomposition (VMD) [18]. The following is a brief description of how to use various decomposition algorithms to decompose features and extract hidden features.(2) Empirical Mode Decomposition. EMD is an adaptive time-frequency signal processing algorithm proposed by Huang E in 1998 [19]. A series of applications have proved that EMD is especially suitable for nonlinear and nonstationary signal processing with noise. EMD decomposes the original signal x (t) into several Intrinsic Mode Functions (IMF) and trend terms, (1)xt=∑j=1ncj+rn,where cj is IMF and rn is trend term. IMF is a random oscillation function with different amplitudes and frequencies and meets the following two characteristics: (1) the number of extreme points must be equal to the number of zero crossings, or the maximum difference is 1; (2) the average of the upper envelope (defined by the maximum) and the lower envelope (defined by the minimum) is 0. The EMD calculation process is as follows:(1)
Identify the local maxima and minima in the original sequencex(t), connect the local maxima points by cubic spline function, and obtain the upper and lower envelopes xmax(t) and xmin(t).(2)
Calculate the differencec1 (t) between the original sequence x(t) and the mean sequence m(t) of the upper and lower envelopes(2)mt=xmaxt+xmint2,c1t=xt−mt.(3)
Replace the original sequence withc1(t), and repeat steps 1-2 until the envelope is symmetric, the mean value is zero (i.e., the above two characteristics of IMF are satisfied), the residual residuals are monotonous, and IMF cannot be separated from the residual sequence.The disadvantage of the EMD algorithm is that modal aliasing often occurs, which leads to incomplete IMF separation; that is, multiple IMF components contain repeated information and have a high correlation among components. In addition, because EMD is adaptive, the decomposition level is determined by the algorithm according to the data characteristics. For time series decomposition, the decomposition level is also variable, so it cannot be used for hydrometeorological time series prediction.(3) Wavelet Transform. WT is a mathematical operation that can perform convolution operation on time series or signals in time domain and frequency domain at the same time. Discrete wavelet transform (DWT) is commonly used in hydrology [20]. Given the original signal x(t), DWT can be defined as(3)DWTj,k=2−j/2∫−∞+∞xtψ∗2−jt−kdt,where k is the position index, j is the decomposition level, and ψ∗ is the wavelet function. DWT operates a high-pass filter and a low-pass filter to decompose the original time series into detail subseries representing high-frequency components and approximation subseries representing trends (low-frequency). The specific calculation process is as follows:(1)
The original time seriesx(t) and the decomposition layer j in the wavelet function ψ∗ are given.(2)
DWT decomposesx(t) into a high-frequency detail component (D1(t)) and a low-frequency approximation component (A1(t)) at the first decomposition level, and(4)xt=D1t+A1t.(3)
At the second decomposition level, DWT continues to decomposeA1(t) into a high-frequency detail component (D2(t)) and a low-frequency approximation component (A2(t)), and(5)At=D2t+A2t.(4)
Similar to Step 3, at each subsequent decomposition level, the low-frequency approximate component obtained at the previous decomposition level is decomposed into the high-frequency detail component and the low-frequency approximate component until a preset decomposition level is reached.The original time series is expressed as the sum of approximate components and all detail components; that is,(6)xt=D1t+⋯+Djt+Ajt.The disadvantage of DWT is that there is no uniform standard and theory to determine the wavelet function and decomposition level, and a large number of experiments are needed.(4) Variational Mode Decomposition. VMD is a signal processing algorithm proposed in 2014, which can decompose the original signal into K IMF at one time [21]. The basic idea of VMD is to construct a variational problem and solve the variational problem to obtain decomposed subsignals. Given the input signal f(t), the variational problem can be defined as (7)minukwk∑k∂tδt+jπt∗ukte−jωkt22s.t.∑kukt=ft,where {uk} and {wk} represent the set of modes and their corresponding center frequencies. T represents time, j2 = −1 represents the square of imaginary units, ∗ represents the convolution operator, and σ represents the Dirac trigonometric function.For VMD, the alternating direction multiplier method is used to solve the formula, and the updated formulas of modaluk(w), center frequency wk, and Lagrangian multiplier λ in frequency domain are obtained as follows:(8)u^kn+1ω=f^ω−∑i<ku^kn+1ω+λ^nω/21+2aω−ωk2,ω^kn+1=∫0∞ωu^kn+1ω2dω∫0∞u^kn+1ω2dω,λ^n+1ω=λ^nω+τf^ω−∑ku^kn+1ω,where n is the iterative counter, τ is the noise tolerance, and u^kn+1ω, f^ω, and λ^nω are the Fourier transforms of ukn+1ω, uω, and unω, respectively. Variables u^kn+1ω, f^ω, and λ^nω are continuously updated until they are less than the convergence error ε:(9)∑kμkn+1−μkn22μkn22<ε.The decomposition effect of VMD is affected byK. A small K value cannot extract IMF effectively from the original signal, and a large K value may lead to redundancy of IMF information. A smaller alpha value may lead to larger bandwidth, redundant information expression, and increased additional noise. On the contrary, a larger alpha value may lead to smaller bandwidth and loss of effective information. However, after a large number of experiments, it is determined that when VMD is used to decompose runoff time series, when α = 2000, τ = 0, and ε = 1 × 10−9, it can ensure a good denoising effect and effective IMF separation. In this study, K is optimized by observing whether there is center frequency aliasing, that is, starting from K = 2 (K value increases by 1 every time), and stopping optimization until the last VMD component shows center frequency aliasing.


## 2.1.1. Feature Factor Extraction

According to the theory of runoff generation and confluence [16], the precipitation, evaporation, and environmental characteristics of the underlying surface in the upstream catchment area controlled by a certain section of a river have a direct impact on the runoff formation process of the section. In addition, direct human activities, such as water intake, water use, water consumption, drainage, and water transfer, can change the temporal and spatial distribution of water resources, thus affecting the runoff formation of target sections. Therefore, according to the information of catchment area, distribution of hydrometeorological stations, and distribution of water intakes and outlets, based on the data sets related to climate change, underlying surface environmental evolution and human activities, the hydrometeorological time series, and underlying surface environmental evolution time series and direct human activity time series driving runoff formation are extracted to extract runoff driving factors and construct runoff prediction sample set.


## 2.1.2. Processing of Missing Value of Characteristic Factor

Interpolation schematic diagram of piecewise cubic spline method for missing daily flow is shown in Figure1.Figure 1
Interpolation schematic diagram of piecewise cubic spline method for missing daily flow.There are generally two ways to deal with missing values [17]:(1)
When the continuous deletion period of the feature sequence is too long or the number of missing values reaches more than 50% of the total period length, the feature sequence is directly removed.(2)
When there are scattered missing values, an interpolation method is used to make up. Commonly used time series interpolation methods include mean interpolation, correlation analysis (unary linear regression, multiple linear regression, and iterative regression), and spline interpolation. Under the disturbance of changing environment, the statistical characteristics of hydrometeorological time series are constantly changing, so the mean interpolation method is no longer applicable. When the correlation between interpolation sequence and other reference sequences is small, the effect of correlation analysis is poor. The runoff data (see Table1) and meteorological data (see Table 2) used in this study are partially missing or missing, so it is necessary to interpolate the missing values.Table 1
Basic information of hydrological stations in the study basin.
WatershedWatershed areaAverage annual runoffAverage sediment amount for many yearsYanhe River Basin7725 square kilometers2.93 billion cubic meters244∼311 kg/m3Table 2
Basic information of meteorological observation variables.
Variable nameVariable abbreviationUnitAverage station pressureAVG_PShPaDaily maximum station pressureMAX_PShPaDaily minimum station pressureMIN_PShPaAverage temperatureAVG_T°CDaily maximum temperatureMAX_T°CDaily minimum temperatureMIN_T°CRelative humidityRHU%Minimum relative humidityMIN_RHU%Precipitation at 20-8 o’clockP208mmPrecipitation at 8-20 o’clockP820mmPrecipitation at 20-20 o’clockP2020mmSmall evaporation capacityS_EVPmmLarge evaporation capacityL_EVPmmAverage wind speedAVG_Wm/sMaximum wind speedMAX_Wm/sAverage surface temperatureAVG_ST°CMaximum daily surface temperatureMAX_ST°CDaily minimum surface temperatureMIN_ST°C


## 2.1.3. Complex Nonlinear Eigenfactor Decomposition

(1) Commonly Used Signal Decomposition Algorithms. Decomposition and noise reduction of complex nonlinear, highly nonstationary and multiscale variation feature series, extraction of hidden features in the original series, and construction of runoff prediction model can significantly improve the accuracy of runoff prediction. Commonly used signal decomposition algorithms include Empirical Mode Decomposition (EMD), Wavelet Transform (WT), and Variational Mode Decomposition (VMD) [18]. The following is a brief description of how to use various decomposition algorithms to decompose features and extract hidden features.(2) Empirical Mode Decomposition. EMD is an adaptive time-frequency signal processing algorithm proposed by Huang E in 1998 [19]. A series of applications have proved that EMD is especially suitable for nonlinear and nonstationary signal processing with noise. EMD decomposes the original signal x (t) into several Intrinsic Mode Functions (IMF) and trend terms, (1)xt=∑j=1ncj+rn,where cj is IMF and rn is trend term. IMF is a random oscillation function with different amplitudes and frequencies and meets the following two characteristics: (1) the number of extreme points must be equal to the number of zero crossings, or the maximum difference is 1; (2) the average of the upper envelope (defined by the maximum) and the lower envelope (defined by the minimum) is 0. The EMD calculation process is as follows:(1)
Identify the local maxima and minima in the original sequencex(t), connect the local maxima points by cubic spline function, and obtain the upper and lower envelopes xmax(t) and xmin(t).(2)
Calculate the differencec1 (t) between the original sequence x(t) and the mean sequence m(t) of the upper and lower envelopes(2)mt=xmaxt+xmint2,c1t=xt−mt.(3)
Replace the original sequence withc1(t), and repeat steps 1-2 until the envelope is symmetric, the mean value is zero (i.e., the above two characteristics of IMF are satisfied), the residual residuals are monotonous, and IMF cannot be separated from the residual sequence.The disadvantage of the EMD algorithm is that modal aliasing often occurs, which leads to incomplete IMF separation; that is, multiple IMF components contain repeated information and have a high correlation among components. In addition, because EMD is adaptive, the decomposition level is determined by the algorithm according to the data characteristics. For time series decomposition, the decomposition level is also variable, so it cannot be used for hydrometeorological time series prediction.(3) Wavelet Transform. WT is a mathematical operation that can perform convolution operation on time series or signals in time domain and frequency domain at the same time. Discrete wavelet transform (DWT) is commonly used in hydrology [20]. Given the original signal x(t), DWT can be defined as(3)DWTj,k=2−j/2∫−∞+∞xtψ∗2−jt−kdt,where k is the position index, j is the decomposition level, and ψ∗ is the wavelet function. DWT operates a high-pass filter and a low-pass filter to decompose the original time series into detail subseries representing high-frequency components and approximation subseries representing trends (low-frequency). The specific calculation process is as follows:(1)
The original time seriesx(t) and the decomposition layer j in the wavelet function ψ∗ are given.(2)
DWT decomposesx(t) into a high-frequency detail component (D1(t)) and a low-frequency approximation component (A1(t)) at the first decomposition level, and(4)xt=D1t+A1t.(3)
At the second decomposition level, DWT continues to decomposeA1(t) into a high-frequency detail component (D2(t)) and a low-frequency approximation component (A2(t)), and(5)At=D2t+A2t.(4)
Similar to Step 3, at each subsequent decomposition level, the low-frequency approximate component obtained at the previous decomposition level is decomposed into the high-frequency detail component and the low-frequency approximate component until a preset decomposition level is reached.The original time series is expressed as the sum of approximate components and all detail components; that is,(6)xt=D1t+⋯+Djt+Ajt.The disadvantage of DWT is that there is no uniform standard and theory to determine the wavelet function and decomposition level, and a large number of experiments are needed.(4) Variational Mode Decomposition. VMD is a signal processing algorithm proposed in 2014, which can decompose the original signal into K IMF at one time [21]. The basic idea of VMD is to construct a variational problem and solve the variational problem to obtain decomposed subsignals. Given the input signal f(t), the variational problem can be defined as (7)minukwk∑k∂tδt+jπt∗ukte−jωkt22s.t.∑kukt=ft,where {uk} and {wk} represent the set of modes and their corresponding center frequencies. T represents time, j2 = −1 represents the square of imaginary units, ∗ represents the convolution operator, and σ represents the Dirac trigonometric function.For VMD, the alternating direction multiplier method is used to solve the formula, and the updated formulas of modaluk(w), center frequency wk, and Lagrangian multiplier λ in frequency domain are obtained as follows:(8)u^kn+1ω=f^ω−∑i<ku^kn+1ω+λ^nω/21+2aω−ωk2,ω^kn+1=∫0∞ωu^kn+1ω2dω∫0∞u^kn+1ω2dω,λ^n+1ω=λ^nω+τf^ω−∑ku^kn+1ω,where n is the iterative counter, τ is the noise tolerance, and u^kn+1ω, f^ω, and λ^nω are the Fourier transforms of ukn+1ω, uω, and unω, respectively. Variables u^kn+1ω, f^ω, and λ^nω are continuously updated until they are less than the convergence error ε:(9)∑kμkn+1−μkn22μkn22<ε.The decomposition effect of VMD is affected byK. A small K value cannot extract IMF effectively from the original signal, and a large K value may lead to redundancy of IMF information. A smaller alpha value may lead to larger bandwidth, redundant information expression, and increased additional noise. On the contrary, a larger alpha value may lead to smaller bandwidth and loss of effective information. However, after a large number of experiments, it is determined that when VMD is used to decompose runoff time series, when α = 2000, τ = 0, and ε = 1 × 10−9, it can ensure a good denoising effect and effective IMF separation. In this study, K is optimized by observing whether there is center frequency aliasing, that is, starting from K = 2 (K value increases by 1 every time), and stopping optimization until the last VMD component shows center frequency aliasing.


## 2.2. Characteristic Factor Cleaning

The purpose of feature factor cleaning is (1) to exclude the feature factors that contribute 0 to the prediction target and (2) to reconstruct the relevant feature factors. Feature cleaning can effectively remove the insignificant features, reduce the burden for the subsequent quantitative calculation of driving factors, and effectively avoid the risk of overfitting the machine learning model of runoff prediction. There are two ways to realize feature cleaning [22]: (1) directly remove the features with a variance of 0, analyze the linear correlation or multicollinearity between the remaining features, and reconstruct the related features; (2) principal component analysis or cluster analysis is directly used to reduce the dimension of features and realize feature cleaning. The variance is 0, which means that the value of the feature does not change with time, so it has no influence on the prediction target. The smaller the characteristic variance, the smaller the influence on the prediction target and the smaller the contribution. Figures 2 and 3, respectively, give the variance information of meteorological characteristics and ERA5L characteristics in Yanhe River Basin. It can be seen that although there is no characteristic with a variance of 0, there are a large number of small variance characteristics, which need further cleaning.Figure 2
Variance of meteorological characteristics from the source of the Yellow River Basin to Longyangxia.Figure 3
ERA5L characteristic variance from the source of the Yellow River Basin to Longyangxia.VIF first establishes a regression equation between one feature and other features and then calculates based on the following formula:(10)VIF=11−R2.The original feature set is reorganized into a main component set. Principal components are defined as(11)Zi=ai1∗X1+ai2∗X2+⋯+aiP∗Xp,where Zi represents principal components, ai1, ai2, …, aiP represent feature vectors, X1, X2, …, Xp represent original features, and p represents the number of features obtained from Zi.


## 2.3. Extraction of Driving Factors Based on Genetic Contribution Analysis

Under the changing environment, the factors driving runoff change are constantly changing, so it is necessary to analyze the contribution of driving factors to runoff and determine the dominant factors. In addition, containing a large number of driving factors in the sample set will increase the risk of model overfitting, reduce the generalization prediction ability of the model, and also increase the consumption of modeling resources and time. Therefore, it is necessary to screen key driving factors and exclude driving factors that have less impact on runoff. In this paper, the key driving factors are screened based on the contribution of driving factors to the prediction target. The following is a study on how to quantitatively calculate the contribution of driving factors.Mutual information is a measure of the interdependence between two random variables, which is closely related to the origin of random variables. For two joint discrete random variablesX and Y, their mutual information is defined as(12)IX;Y=∑y∈Y∑x∈XpX,Yx,ylogpX,Yx,ypXxPYy,where P(X,Y)(x, y) is the joint distribution function of X and Y and PX(x) and PY(y) are the edge distributions of X and Y, respectively.


## 3. Study on Ecological Water Consumption Evaluation Model Based on Machine Learning

Based on the adaptive prediction of ecological water consumption [3], this paper studies the construction method and model interpretation of ecological water consumption prediction model based on machine learning on the basis of runoff statistics and causes based on big data and runoff adaptive prediction model [23].


### 3.1. Sample Preprocessing

#### 3.1.1. Sample Set Partition

Ideally, the sample set should be divided into three subsets: training set, development set, and test set. The training set is used to correct model parameters, and the development set is used to optimize model superparameters, screen features, make decisions, and so on. The test set tests the final optimized model and provides confidence level estimation for the use of the model.In addition, the partition ratio of the training set, development set, and test set also affects the prediction effect of the model. At present, there is no uniform standard for the division ratio of the training set, development set, and test set. When the sample space is larger, the sample is more representative of reality, and only a small number of development and test samples are needed at this time. When the sample space is smaller, the actual uncertainty of the sample description is greater. At this time, it is necessary to increase the proportion of development set and test set to improve the persuasiveness of the model. For runoff prediction, this paper gives an approximate empirical range: when the sample size is less than or equal to 500, the ratio of the training set, development set, and test set can be set to 50%, 25%, and 25%; when the sample size is greater than 500 and less than or equal to 1000, the ratio of the training set, development set, and test set can be set to 70%, 15%, and 15%; when the sample size is greater than 1000 and less than or equal to 5000, the ratio of the training set, development set, and test set can be set to 80%, 10%, and 10%; when the sample size is greater than 5000 and less than 10000, the ratio of the training set, development set, and test set can be set to 90%, 5%, and 5%; when the data sample size continues to increase, the proportion of development set and test set can continue to decrease.


#### 3.1.2. Normalization of Sample Set

The magnitude or value range of predictor sequence is far from each other, which leads to the failure of the objective function optimization algorithm in the machine learning model. Sample normalization makes the prediction factor sequence have the same influence on the objective function, which can accelerate the optimal convergence. At present, the commonly used normalization methods include linear normalization and mean normalization. Linear normalization can scale the input feature and output target to [−1, 1] or [0, 1], and the calculation formulas are as follows:(13)x′=2∗x−xminxmax−xmin−1,(14)x′=x−xminxmax−xmin,where x and x′ are the original value and normalized value, respectively, and xmax and xmin are the maximum value and minimum value of the original sequence, respectively. Generally, when the original sequence contains negative values, formula (13) is used; otherwise, formula (14) is used. There is no definite range of mean value normalization, and its calculation formula is(15)x′=x−xmeanxsd,where xmean and xsd are the mean and standard deviation of the original sequence, respectively.


### 3.2. Construction of Machine Learning Model for Ecological Water Consumption Assessment

Gradient enhancement (GB) is a powerful machine learning strategy, which can efficiently build robust and competitive models for classification or regression problems. The basic idea of Boosting is to combine many weak learner outputs into a powerful integration model. GB builds the model through a forward phase-by-phase approach. The integration modelFm(x) of GB stage m (1 ≤ m ≤ M, M is the number of stages and the number of weak learners) is defined as(16)Fmx=Fm−1x+hmx,where hm(x) represents a weak learner, and in the Gradient Boosting Regression Tree (GBRT), it is a single decision regression tree. At the same time, the prediction target y of the training sample is evaluated based on hm(x),(17)Fmx=Fm−1x+hmx=y.It can be further converted to(18)hmx=y−Fm−1x.Therefore,hm(x) is used to fit the residual γm = y −Fm−1(x) of the current stage model, and the residual is also the square error loss function of the model:(19)−a1/2y−Fm−1x2∂Fm−1x=y−Fm−1x.


#### 3.2.1. Model Hyperparameter Optimization

There are many methods to optimize the superparameters of machine learning models, such as trial-and-error method, grid search, Bayesian optimization, genetic algorithm, and gravity search algorithm. In this paper, the trial-and-error method and Bayesian optimization are used as two superparametric optimization methods, and their implementation ideas are introduced, respectively, in the following.Bayesian optimization (BO) is a serialized model optimization method, which is usually applied to the optimization of the black-box objective function with unknown real distribution or is very difficult to solve. When optimizing the superparameters of the machine learning model, BO firstly sets a priori belief for the loss function and continuously optimizes the model by continuously calculating the loss function value and updating the Bayesian posterior estimation in Figure4.Figure 4
Flowchart of BO algorithm.


## 3.1. Sample Preprocessing

### 3.1.1. Sample Set Partition

Ideally, the sample set should be divided into three subsets: training set, development set, and test set. The training set is used to correct model parameters, and the development set is used to optimize model superparameters, screen features, make decisions, and so on. The test set tests the final optimized model and provides confidence level estimation for the use of the model.In addition, the partition ratio of the training set, development set, and test set also affects the prediction effect of the model. At present, there is no uniform standard for the division ratio of the training set, development set, and test set. When the sample space is larger, the sample is more representative of reality, and only a small number of development and test samples are needed at this time. When the sample space is smaller, the actual uncertainty of the sample description is greater. At this time, it is necessary to increase the proportion of development set and test set to improve the persuasiveness of the model. For runoff prediction, this paper gives an approximate empirical range: when the sample size is less than or equal to 500, the ratio of the training set, development set, and test set can be set to 50%, 25%, and 25%; when the sample size is greater than 500 and less than or equal to 1000, the ratio of the training set, development set, and test set can be set to 70%, 15%, and 15%; when the sample size is greater than 1000 and less than or equal to 5000, the ratio of the training set, development set, and test set can be set to 80%, 10%, and 10%; when the sample size is greater than 5000 and less than 10000, the ratio of the training set, development set, and test set can be set to 90%, 5%, and 5%; when the data sample size continues to increase, the proportion of development set and test set can continue to decrease.


### 3.1.2. Normalization of Sample Set

The magnitude or value range of predictor sequence is far from each other, which leads to the failure of the objective function optimization algorithm in the machine learning model. Sample normalization makes the prediction factor sequence have the same influence on the objective function, which can accelerate the optimal convergence. At present, the commonly used normalization methods include linear normalization and mean normalization. Linear normalization can scale the input feature and output target to [−1, 1] or [0, 1], and the calculation formulas are as follows:(13)x′=2∗x−xminxmax−xmin−1,(14)x′=x−xminxmax−xmin,where x and x′ are the original value and normalized value, respectively, and xmax and xmin are the maximum value and minimum value of the original sequence, respectively. Generally, when the original sequence contains negative values, formula (13) is used; otherwise, formula (14) is used. There is no definite range of mean value normalization, and its calculation formula is(15)x′=x−xmeanxsd,where xmean and xsd are the mean and standard deviation of the original sequence, respectively.


## 3.1.1. Sample Set Partition

Ideally, the sample set should be divided into three subsets: training set, development set, and test set. The training set is used to correct model parameters, and the development set is used to optimize model superparameters, screen features, make decisions, and so on. The test set tests the final optimized model and provides confidence level estimation for the use of the model.In addition, the partition ratio of the training set, development set, and test set also affects the prediction effect of the model. At present, there is no uniform standard for the division ratio of the training set, development set, and test set. When the sample space is larger, the sample is more representative of reality, and only a small number of development and test samples are needed at this time. When the sample space is smaller, the actual uncertainty of the sample description is greater. At this time, it is necessary to increase the proportion of development set and test set to improve the persuasiveness of the model. For runoff prediction, this paper gives an approximate empirical range: when the sample size is less than or equal to 500, the ratio of the training set, development set, and test set can be set to 50%, 25%, and 25%; when the sample size is greater than 500 and less than or equal to 1000, the ratio of the training set, development set, and test set can be set to 70%, 15%, and 15%; when the sample size is greater than 1000 and less than or equal to 5000, the ratio of the training set, development set, and test set can be set to 80%, 10%, and 10%; when the sample size is greater than 5000 and less than 10000, the ratio of the training set, development set, and test set can be set to 90%, 5%, and 5%; when the data sample size continues to increase, the proportion of development set and test set can continue to decrease.


## 3.1.2. Normalization of Sample Set

The magnitude or value range of predictor sequence is far from each other, which leads to the failure of the objective function optimization algorithm in the machine learning model. Sample normalization makes the prediction factor sequence have the same influence on the objective function, which can accelerate the optimal convergence. At present, the commonly used normalization methods include linear normalization and mean normalization. Linear normalization can scale the input feature and output target to [−1, 1] or [0, 1], and the calculation formulas are as follows:(13)x′=2∗x−xminxmax−xmin−1,(14)x′=x−xminxmax−xmin,where x and x′ are the original value and normalized value, respectively, and xmax and xmin are the maximum value and minimum value of the original sequence, respectively. Generally, when the original sequence contains negative values, formula (13) is used; otherwise, formula (14) is used. There is no definite range of mean value normalization, and its calculation formula is(15)x′=x−xmeanxsd,where xmean and xsd are the mean and standard deviation of the original sequence, respectively.


## 3.2. Construction of Machine Learning Model for Ecological Water Consumption Assessment

Gradient enhancement (GB) is a powerful machine learning strategy, which can efficiently build robust and competitive models for classification or regression problems. The basic idea of Boosting is to combine many weak learner outputs into a powerful integration model. GB builds the model through a forward phase-by-phase approach. The integration modelFm(x) of GB stage m (1 ≤ m ≤ M, M is the number of stages and the number of weak learners) is defined as(16)Fmx=Fm−1x+hmx,where hm(x) represents a weak learner, and in the Gradient Boosting Regression Tree (GBRT), it is a single decision regression tree. At the same time, the prediction target y of the training sample is evaluated based on hm(x),(17)Fmx=Fm−1x+hmx=y.It can be further converted to(18)hmx=y−Fm−1x.Therefore,hm(x) is used to fit the residual γm = y −Fm−1(x) of the current stage model, and the residual is also the square error loss function of the model:(19)−a1/2y−Fm−1x2∂Fm−1x=y−Fm−1x.


### 3.2.1. Model Hyperparameter Optimization

There are many methods to optimize the superparameters of machine learning models, such as trial-and-error method, grid search, Bayesian optimization, genetic algorithm, and gravity search algorithm. In this paper, the trial-and-error method and Bayesian optimization are used as two superparametric optimization methods, and their implementation ideas are introduced, respectively, in the following.Bayesian optimization (BO) is a serialized model optimization method, which is usually applied to the optimization of the black-box objective function with unknown real distribution or is very difficult to solve. When optimizing the superparameters of the machine learning model, BO firstly sets a priori belief for the loss function and continuously optimizes the model by continuously calculating the loss function value and updating the Bayesian posterior estimation in Figure4.Figure 4
Flowchart of BO algorithm.


## 3.2.1. Model Hyperparameter Optimization

There are many methods to optimize the superparameters of machine learning models, such as trial-and-error method, grid search, Bayesian optimization, genetic algorithm, and gravity search algorithm. In this paper, the trial-and-error method and Bayesian optimization are used as two superparametric optimization methods, and their implementation ideas are introduced, respectively, in the following.Bayesian optimization (BO) is a serialized model optimization method, which is usually applied to the optimization of the black-box objective function with unknown real distribution or is very difficult to solve. When optimizing the superparameters of the machine learning model, BO firstly sets a priori belief for the loss function and continuously optimizes the model by continuously calculating the loss function value and updating the Bayesian posterior estimation in Figure4.Figure 4
Flowchart of BO algorithm.


## 4. Experiment

### 4.1. Example Verification and Comparative Evaluation of Ecological Water Consumption Prediction in Yanhe River Basin

#### 4.1.1. Comparison of Screening Thresholds for Different Driving Factors

In order to explain the effect of screening thresholds of different key driving factors on runoff prediction, based on the CD1 sample set, this paper optimizes the LSTM model by BO and evaluates the water consumption of related cities in Yanhe River Basin, with a forecast period of one month. In the CD1 sample set, according to the normalized mutual information between predictors and runoff, 10 screening thresholds were set; that is, the normalized mutual information was 0.0∼0.9, and the interval was 0.1.The results in Figure5 show that, before screening key driving factors, it is necessary to test the screening threshold to avoid losing runoff driving information. At the same time, human activities have a great influence on the formation and evolution of ecological water consumption. This is the key problem to accurately assess ecological water consumption at present: whether it is to restore the measured ecological water consumption, exclude the influence of human activities, or incorporate the influence of human activities into the construction of ecological water consumption prediction sample set, it is limited by objective factors such as limited data, inaccurate observation, and mismatch between time and space scales.Figure 5
LSTM model NSE, NRMSE, and PPTS equation (5) values based on CD1 sample set.


#### 4.1.2. Comparison of Forecast Effects of Meteorological and ERA5L Elements

In order to illustrate the correction effect of ecological water consumption data on the model in the same period of history and compare the runoff prediction effect of meteorological data, ERA5L data, and mixed data of meteorological data and ERA5L, this paper optimizes the LSTM model by BO based on sample set and forecasts it with a forecast period of one month. Part of the sample set uses the historical ecological water consumption participation model to correct the prediction results, while the other part does not use the historical ecological water consumption correction model. Figure6 gives the NSE, NRMSE, and PPTS equation (5) box diagrams of the LSTM model under Yan’an and Dingbian one-month forecast period, and the box diagrams are drawn based on the corresponding evaluation indexes of the training set, development set, and test set.Figure 6
LSTM model NSE, NRMSE, and PPTS equation (5) box diagram of Yan’an and fixed edge one-month forecast period.The results in Figure6 show that the prediction effect of ecological water consumption by linear correlation feature reconstruction is better than that by feature dimension reduction as a whole. The reasons for this result may be as follows: (1) Feature dimension reduction excludes some principal components with small variance, resulting in information loss to a certain extent. (2) Feature dimension reduction selects more meteorological data principal components and fewer ERA5L data principal components, while linear correlation feature reconstruction results in screening fewer meteorological features and more ERA5L features. The prediction effect of ERA5L is slightly better than that of meteorological data, so a higher ERA5L information ratio may lead to relatively better prediction accuracy.


### 4.2. Assessment of Ecological Environment Water Consumption in Yanhe River Basin

According to the comparison in Table3, this subsection is also based on the ecological water consumption of the basin obtained by big data analysis. In order to evaluate the ecological water consumption of Yanhe River Basin, the data needed are not only deep groundwater recharge, long-distance water transfer (Yellow River Diversion), domestic sewage and industrial sewage discharge, and industrial and domestic water intake but also the total amount of local water resources.Table 3
Urban ecological water consumption in Yanhe River Basin (ten thousand cubic meters).
YearZhi DanYan’anYan ChangAn SaiDing BianHeng ShanSui DeLuo Chuan20001189.21250.3763.2644.71307.41274.41249.81231.520011171.61270.5784.76771287.21306.81246.11248.520021141.51348.3792.2652.81278.11304.41252.61242.920031276.21236.2719.5645.91397.81381.61342.81322.5200412851162.9753.3683.91402.61406.913701389.120051285.51274.1731.8701.41404.31376.61326.6137520061277.81233.8700.9701.21377.41337.81309.41356.320071166.21360.7701.5661.91280.31307.11266.91222.420081138.91271694.2689.51288.31290.71256.71208.420091221.81217.6688.4622.91313.71340.71293.61205.8201012691230730.8620.41414.713701324.81289.820111231.31232.1659.7639.31365.41408.413831289.120121305.71270.5666.1658.91423.91422.113851290201312391232.6700.3670.813611341.61314.71320.320141273.61158.4670688.71400.21405.71338.81294.320151278.91228.3698.6702.91389.11417.81375.31393.220161275.31217.9701.97011392.41394.31317.71211.420171220.41227.2716.8696.31354.11344.51298.91226.120181270.21201.1694.7682.31354.61394.81334.81357.920191163.31198.3679.4649.31270.21270.21251.11147.5202012641161.6810.2686.91373.51377.51325.11260.5Total25944.425983.415058.21407828436.228473.927563.726882.5The results in Figure7 show that the ecological water consumption of the basin defined based on big data analysis in this paper belongs to the broad ecological water consumption, which depicts the total amount of water consumed from rainfall in order to realize the ecological functions of the basin and represents the natural attributes of the basin. Due to the large buffer capacity of nature, the change of the total amount is insensitive to the change of environment.Figure 7
Local water resources calculated after standardization.
(a)(b)(c)(d)


## 4.1. Example Verification and Comparative Evaluation of Ecological Water Consumption Prediction in Yanhe River Basin

### 4.1.1. Comparison of Screening Thresholds for Different Driving Factors

In order to explain the effect of screening thresholds of different key driving factors on runoff prediction, based on the CD1 sample set, this paper optimizes the LSTM model by BO and evaluates the water consumption of related cities in Yanhe River Basin, with a forecast period of one month. In the CD1 sample set, according to the normalized mutual information between predictors and runoff, 10 screening thresholds were set; that is, the normalized mutual information was 0.0∼0.9, and the interval was 0.1.The results in Figure5 show that, before screening key driving factors, it is necessary to test the screening threshold to avoid losing runoff driving information. At the same time, human activities have a great influence on the formation and evolution of ecological water consumption. This is the key problem to accurately assess ecological water consumption at present: whether it is to restore the measured ecological water consumption, exclude the influence of human activities, or incorporate the influence of human activities into the construction of ecological water consumption prediction sample set, it is limited by objective factors such as limited data, inaccurate observation, and mismatch between time and space scales.Figure 5
LSTM model NSE, NRMSE, and PPTS equation (5) values based on CD1 sample set.


### 4.1.2. Comparison of Forecast Effects of Meteorological and ERA5L Elements

In order to illustrate the correction effect of ecological water consumption data on the model in the same period of history and compare the runoff prediction effect of meteorological data, ERA5L data, and mixed data of meteorological data and ERA5L, this paper optimizes the LSTM model by BO based on sample set and forecasts it with a forecast period of one month. Part of the sample set uses the historical ecological water consumption participation model to correct the prediction results, while the other part does not use the historical ecological water consumption correction model. Figure6 gives the NSE, NRMSE, and PPTS equation (5) box diagrams of the LSTM model under Yan’an and Dingbian one-month forecast period, and the box diagrams are drawn based on the corresponding evaluation indexes of the training set, development set, and test set.Figure 6
LSTM model NSE, NRMSE, and PPTS equation (5) box diagram of Yan’an and fixed edge one-month forecast period.The results in Figure6 show that the prediction effect of ecological water consumption by linear correlation feature reconstruction is better than that by feature dimension reduction as a whole. The reasons for this result may be as follows: (1) Feature dimension reduction excludes some principal components with small variance, resulting in information loss to a certain extent. (2) Feature dimension reduction selects more meteorological data principal components and fewer ERA5L data principal components, while linear correlation feature reconstruction results in screening fewer meteorological features and more ERA5L features. The prediction effect of ERA5L is slightly better than that of meteorological data, so a higher ERA5L information ratio may lead to relatively better prediction accuracy.


## 4.1.1. Comparison of Screening Thresholds for Different Driving Factors

In order to explain the effect of screening thresholds of different key driving factors on runoff prediction, based on the CD1 sample set, this paper optimizes the LSTM model by BO and evaluates the water consumption of related cities in Yanhe River Basin, with a forecast period of one month. In the CD1 sample set, according to the normalized mutual information between predictors and runoff, 10 screening thresholds were set; that is, the normalized mutual information was 0.0∼0.9, and the interval was 0.1.The results in Figure5 show that, before screening key driving factors, it is necessary to test the screening threshold to avoid losing runoff driving information. At the same time, human activities have a great influence on the formation and evolution of ecological water consumption. This is the key problem to accurately assess ecological water consumption at present: whether it is to restore the measured ecological water consumption, exclude the influence of human activities, or incorporate the influence of human activities into the construction of ecological water consumption prediction sample set, it is limited by objective factors such as limited data, inaccurate observation, and mismatch between time and space scales.Figure 5
LSTM model NSE, NRMSE, and PPTS equation (5) values based on CD1 sample set.


## 4.1.2. Comparison of Forecast Effects of Meteorological and ERA5L Elements

In order to illustrate the correction effect of ecological water consumption data on the model in the same period of history and compare the runoff prediction effect of meteorological data, ERA5L data, and mixed data of meteorological data and ERA5L, this paper optimizes the LSTM model by BO based on sample set and forecasts it with a forecast period of one month. Part of the sample set uses the historical ecological water consumption participation model to correct the prediction results, while the other part does not use the historical ecological water consumption correction model. Figure6 gives the NSE, NRMSE, and PPTS equation (5) box diagrams of the LSTM model under Yan’an and Dingbian one-month forecast period, and the box diagrams are drawn based on the corresponding evaluation indexes of the training set, development set, and test set.Figure 6
LSTM model NSE, NRMSE, and PPTS equation (5) box diagram of Yan’an and fixed edge one-month forecast period.The results in Figure6 show that the prediction effect of ecological water consumption by linear correlation feature reconstruction is better than that by feature dimension reduction as a whole. The reasons for this result may be as follows: (1) Feature dimension reduction excludes some principal components with small variance, resulting in information loss to a certain extent. (2) Feature dimension reduction selects more meteorological data principal components and fewer ERA5L data principal components, while linear correlation feature reconstruction results in screening fewer meteorological features and more ERA5L features. The prediction effect of ERA5L is slightly better than that of meteorological data, so a higher ERA5L information ratio may lead to relatively better prediction accuracy.


## 4.2. Assessment of Ecological Environment Water Consumption in Yanhe River Basin

According to the comparison in Table3, this subsection is also based on the ecological water consumption of the basin obtained by big data analysis. In order to evaluate the ecological water consumption of Yanhe River Basin, the data needed are not only deep groundwater recharge, long-distance water transfer (Yellow River Diversion), domestic sewage and industrial sewage discharge, and industrial and domestic water intake but also the total amount of local water resources.Table 3
Urban ecological water consumption in Yanhe River Basin (ten thousand cubic meters).
YearZhi DanYan’anYan ChangAn SaiDing BianHeng ShanSui DeLuo Chuan20001189.21250.3763.2644.71307.41274.41249.81231.520011171.61270.5784.76771287.21306.81246.11248.520021141.51348.3792.2652.81278.11304.41252.61242.920031276.21236.2719.5645.91397.81381.61342.81322.5200412851162.9753.3683.91402.61406.913701389.120051285.51274.1731.8701.41404.31376.61326.6137520061277.81233.8700.9701.21377.41337.81309.41356.320071166.21360.7701.5661.91280.31307.11266.91222.420081138.91271694.2689.51288.31290.71256.71208.420091221.81217.6688.4622.91313.71340.71293.61205.8201012691230730.8620.41414.713701324.81289.820111231.31232.1659.7639.31365.41408.413831289.120121305.71270.5666.1658.91423.91422.113851290201312391232.6700.3670.813611341.61314.71320.320141273.61158.4670688.71400.21405.71338.81294.320151278.91228.3698.6702.91389.11417.81375.31393.220161275.31217.9701.97011392.41394.31317.71211.420171220.41227.2716.8696.31354.11344.51298.91226.120181270.21201.1694.7682.31354.61394.81334.81357.920191163.31198.3679.4649.31270.21270.21251.11147.5202012641161.6810.2686.91373.51377.51325.11260.5Total25944.425983.415058.21407828436.228473.927563.726882.5The results in Figure7 show that the ecological water consumption of the basin defined based on big data analysis in this paper belongs to the broad ecological water consumption, which depicts the total amount of water consumed from rainfall in order to realize the ecological functions of the basin and represents the natural attributes of the basin. Due to the large buffer capacity of nature, the change of the total amount is insensitive to the change of environment.Figure 7
Local water resources calculated after standardization.
(a)(b)(c)(d)


## 5. Conclusion

Based on big data analysis, this paper proposes a method to evaluate ecological water consumption in river basins by machine learning. Because this method starts from precipitation, which is the original water source of the basin, it avoids the calculation of water resources in the basin which is easy to cause confusion. It is suitable for the calculation of ecological water consumption in various river basins. Taking Yanhe River Basin as an example, the data of annual rainfall, water diversion from Yellow River, groundwater overexploitation, industrial water consumption, agricultural water consumption, domestic water consumption, industrial wastewater discharge, and domestic sewage discharge in Yanhe River Basin are used for calculation. The results show that the restoration of the ecological environment in Yanhe River Basin should not only depend on interbasin water transfer but also strengthen ecological environment protection and rational water allocation. In the future, those factors in the ecological environment that need to be solved have a direct impact on river water consumption, and the internal causes existing in a large number of data can be found. It is necessary to collect more data from different weather and geography in order to analyze more direct reasons affecting water consumption.



---

*Source: 2201964-2021-11-30.xml*