## Abstract

Precise and credible runoff forecasting is extraordinarily vital for various activities of water resources deployment and implementation. The neoteric contribution of the current article is to develop a hybrid model (ANFIS-GPR) based on adaptive neuro-fuzzy inference system (ANFIS) and Gaussian process regression (GPR) for monthly runoff forecasting in the Beiru river of China, and the optimal input schemes of the models are discussed in detail. Firstly, variables related to runoff are selected from the precipitation, soil moisture content, and evaporation as the first set of input schemes according to correlation analysis (CA). Secondly, principal component analysis (PCA) is used to eliminate the redundant information between the original input variables for forming the second set of input schemes. Finally, the runoff is predicted based on different input schemes and different models, and the prediction performance is compared comprehensively. The results show that the input schemes jointly established by CA and PCA (CA-PCA) can greatly improve the prediction accuracy. ANFIS-GPR displays the best forecasting performance among all the peer models. In the single models, the performance of GPR is better than that of ANFIS.

## HIGHLIGHTS

Principal component analysis is used to simplify the prediction factors and improve the prediction accuracy.

The novel contribution of the article is to develop a new hybrid model, i.e., ANFIS-GPR.

The new model (ANFIS-GPR) is more accurate and reliable in predicting the peak discharge.

The research results are more reliable because there is no large water conservancy project in the target basin.

## INTRODUCTION

Improving the accuracy and reliability of runoff prediction plays a key role in various activities of water resources deployment and implementation, such as allocation of water, drought control, flood mitigation, dam planning, and designing for navigation (Allawi & El-Shafie 2016; Yaseen *et al.* 2016; Yuan *et al.* 2018; Liu *et al.* 2021). In particular, runoff prediction at monthly time scales is very significant for reservoir operation, water planning, and irrigation management (Chen *et al.* 2014; Yuan *et al.* 2015; Buizer *et al.* 2016; Liang *et al.* 2017; Hadi & Tombul 2018). Moreover, runoff prediction is an important non-engineering measure to realize the efficient utilization of water resources, so it has always been a hot and core issue in hydrology.

*et al.*2007). Generalizing the basic conditions of the basin and collating the original data are required to apply these models, which greatly increases the workload of runoff prediction and simulation (Yuan

*et al.*2014). What's more, the spatiotemporal heterogeneity of different watersheds greatly influences the applicability of these models, which is not conducive to their popularization. To overcome the shortcomings of the above process-driven models, data-driven models have been widely used in runoff prediction against the background of the rapid development of computer technology (Yuan

*et al.*2018; Achieng & Zhu 2019). The advantage of these data-driven models is that the quantitative relationships between inputs and outputs can be determined by the analysis of historical data without elaborating on the complex physical mechanism. In summary, these models can be roughly divided into two categories: time series (TS) models and machine learning (ML) techniques (He

*et al.*2020). A stable system, which assumes no significant changes in the laws of history and the future, is a prerequisite for applying TS models in prediction (Box

*et al.*2015). Therefore, TS models, e.g., autoregressive (AR), moving average (MA), AR moving average (ARMA), AR integrated moving average (ARIMA), vector autoregression (VAR), and structural VAR (SVAR), are difficult to capture the nonlinear characteristics of runoff TS. Nevertheless, ML models, e.g., artificial neural networks (ANNs) (Abbot & Marohasy 2012; Gholami

*et al.*2015), support vector machine (SVM) (Yu

*et al.*2006; Huang

*et al.*2014), deep belief networks (DBNs) (Xie

*et al.*2019), etc., possess the strong ability of nonlinear function and have been successfully applied to runoff prediction (Maity

*et al.*2010; He

*et al.*2019). Due to the excellent prediction performance of ML models, modified ML models, such as genetic algorithm-tuned backpropagation (GA-BP) (Wang

*et al.*2016), GA-SVM (Su

*et al.*2014), and adaptive neuro-fuzzy inference system (ANFIS) (Jang 1993; Kim & Kasabov 1999), have also been put forward in turn. It is proved that these modified models have further improved the accuracy and stability of runoff prediction (Hui & Ding 2010; Lin & Xiao 2012; Liang

*et al.*2018; Nguyen

*et al.*2019; Feng

*et al.*2020; Jing

*et al.*2020; Yue

*et al.*2020; Ershadi

*et al.*2021). It can be seen from the above-modified models that there are mainly two revised strategies based on ML models: intelligent algorithm tuned neural network (IATNN) (e.g., GA-BP, GA-SVM, and PSO-ELM) and multi-model coupling scheme (MMCS) (e.g., ANFIS). IATNN aims to enhance the models’ learning ability by improving the performance of the intelligent algorithm. However, these models are prone to overfitting in the training period. However, the models built according to MMCS have the advantage of different models being more stable and reliable, except for the poor processing ability of data noise. For example, ANFIS has both the self-learning ability of neural networks and the recognition ability of fuzzy inference systems but no capacity to handle noise. Recently, Gaussian process regression (GPR) has attracted more attention because it has several advantages, such as flexibility, adaptability, and strong generality (Hultquist

*et al.*2014; Sun

*et al.*2014). GPR has a strong noise processing ability, which makes it better compensate for the shortcomings of the ANFIS model in this respect (Garnelo

*et al.*2018; Inyurt

*et al.*2020). Thus, a hybrid model (ANFIS-GPR) is proposed by combining ANFIS and GPR to further improve the prediction accuracy. The ANFIS-GPR is compared with ANFIS and GPR to check its feasibility and effectiveness. In the past, the research on model improvement focused on the adjustment of model structure, and there is little research on the coupling construction of a hybrid model from the perspective of multi-data fusion. Therefore, the current study proposes a conceptual framework for building the ANFIS-GPR model. This model's output is obtained from optimizing the ANFIS and GPR model's output via the Lagrange multiplier method. Three hydro-meteorological variables (i.e., precipitation, soil moisture content (SMC), evaporation) are selected for the candidate inputs based on the physical mechanism of runoff generation. These candidate variables are screened by correlation analysis (CA) and CA–principal component analysis (PCA) to form a variety of input schemes, and the best input schemes are determined after comparative analysis. Finally, a system integrating variable screening and runoff prediction proposed in this paper is verified and discussed based on the Beiru River of Henan province in China (Figure 1). Accordingly, the objectives of this study are as follows:

A novel hybrid ANFIS-GPR model is developed by coupling ANFIS and GPR models based on the Lagrange multiplier method.

The influence of different inputs on the prediction accuracy of the model is evaluated.

Performance of ANFIS-GPR in forecasting runoff is assessed using RMSE, NSC, and CORR.

## METHODOLOGY

### Selection method of predictors

#### Correlation analysis

*et al.*2018). Suppose that two data series and . And the formula for calculating the correlation coefficient (CC) in CA is as follows (Lee & Nicewander 1988; Stigler 1989):where is the average of

*x*, is the average of

_{i}*y*

_{in}, and

*n*is the length of the data series. Values of CC are between [−1,1]. The closer the value of is to 1, the closer the relationship between the two variables is. Otherwise, the contrary is true.

#### Principal component analysis

Principal component analysis (PCA) is a multivariate statistical analysis tool that selects a small number of important variables by a linear transformation of multiple variables. This method can produce several independent principal components (PCs), which can approximately replace the original variables to express the main information (Wold *et al.* 1987).

*P*into

*P*

_{A}and

*P*

_{p}_{-A},

*X*can be separated into two components as follows:where and stand for the principal subspace composed of the first

*A*PCs and the residual subspace composed of the last

*p*−

*A*PCs.

*X*and

_{A}*X*

_{p}_{-A}are the modeled variations and non-modeled variations of

*X*by projection on the

*P*and

_{A}*P*

_{p}_{-A}, respectively.

### Construction of coupled prediction model

#### Adaptive neural-fuzzy inference system

*et al.*1997). Figure 2 shows the structure of ANFIS. Square and circle node properties represent different modules, and the parameters are all located in the square nodes. Each node comes with a node function that changes from node to node. The connection line between two nodes only represents the direction of data transmission.

*x*

_{1},

*x*

_{2},

*x*

_{3}) with the addition of five layers, and there are two fuzzy if-then rules. Normally, the fuzzy rules can be expressed as follows:where

*f*

_{1}and

*f*

_{2}denote the outputs of fuzzy rules, and

*p*

_{1},

*q*

_{1},

*r*

_{1},

*s*

_{1},

*p*

_{2},

*q*

_{2},

*r*

_{2}, and

*s*

_{2}represent the result parameters in the training period, and

*A*

_{1},

*B*

_{1},

*C*

_{1},

*A*

_{2},

*B*

_{2}, and

*C*

_{2}express the membership functions (MFs) for inputs (

*x*

_{1},

*x*

_{2},

*x*

_{3}). A Gaussian MF is a curve that evaluates the strength of a point belonging to the set by distributing a membership degree between 0 and 1, and the calculation formula is shown in the following equation.where

*x*denotes the input value to note

*i*, and

*c*and

*σ*express the center and the width of the Gaussian curve, respectively.

*A*

_{i},

*B*

_{i−2}, and

*C*

_{i−4}are linguistic variables, which are denoted as ‘low’, ‘medium’, and ‘high’ associated with that node.

*O*

_{1,i}is the membership grade of a fuzzy set Z (including

*A*

_{i},

*B*

_{i}, and

*C*

_{i}).

#### Gaussian process regression model

Gaussian process regression (GPR) is an ML tool based on Bayesian theory (Quiñonero-Candela & Rasmussen 2005). The main principle of GPR can be summarized as follows. Assume that the joint distribution of finite random variables obeys Gaussian distribution in the sample space, and the random process satisfies the prior probability of the Gaussian process. Therefore, GPR can follow Bayesian theory to estimate the posterior probability of the random process, and the maximum likelihood method is used to estimate the optimal hyperparameters in this paper.

*x*is the input vector of the predictors, and

_{i}*y*is the output scalar of runoff, a value with Gaussian white noise. The calculation relationship between

_{i}*x*and

_{i}*y*is as follows.with , and is the variance of the distribution.

_{i}*x*is input to GPR, the corresponding output

_{*}*y*can be calculated by Equation (15).

_{*}*f*(

*x*) is the conditional distribution about

_{*}*x*, and

_{*}*y*is the expectation of the conditional distribution, i.e.,where

_{*}*K*(

*x*,

_{i}*x*) is the symmetric positive definite covariance matrix of order

_{i}*n*×

*n*,

*I*is the

_{n}*n*-dimensional identity matrix. Therefore,

*y*and

_{i}*y*must satisfy the joint prior distribution, i.e.,with

_{*}*K*(

*x*,

_{*}*x*) is the symmetric covariance matrix of order

_{i}*n*× 1, and

*k*(

*x*,

_{*}*x*) is the covariance of

_{*}*x*. According to Bayes’ theorem, the posterior distribution of

_{*}*y*can be calculated as:

_{*}*K*. In other words, the form of the kernel function directly determines the prediction performance of GPR. After repeated trials, the mean square exponential expressed in Equation (21) is used as the kernel function in this paper.where

*σ*and

*l*are the hyperparameters that can be determined by the maximum likelihood estimation (MLE) (Perrin

*et al.*2017). In concrete terms, the logarithmic likelihood function of the conditional probability of training samples is established to find the partial derivatives of the hyperparameters, and then the conjugate gradient optimization method is used to get the optimal solution of the hyperparameters. The logarithmic natural function is expressed as

#### ANFIS-GPR model

The construction idea of the ANFIS-GPR model is to minimize the prediction error using optimization. The specific approach can be briefly described as follows. Firstly, ANFIS and GPR are used for runoff prediction to obtain the predicted values, respectively. Secondly, the weight parameters (*w*_{1} and *w*_{2}) of the predicted values of the two models are calculated. Finally, the final predicted values are combined and calculated.

*x*represents the observed value at time

_{t}*t*,

*x*represents the predicted value of the

_{m,t}*m*model at time

*t*, and

*Err*

_{m}_{, t}represents the prediction error of the

*m*model at time

*t*. When

*m*is replaced by

*A*,

*G*, and

*A*-

*G*, it means ANFIS, GPR, and ANFIS-GPR models, respectively. Therefore, there are the following series of expressions.

The key to improving the prediction accuracy of ANFIS-GPR lies in finding appropriate weights *w*_{1} and *w*_{2} to minimize the error from Equation (24), and this is a typical optimization issue. The Lagrange multiplier method is used to solve it (Bertsekas & Dimitri 1982), and the steps are as follows.

*w*

_{1}and

*w*

_{2}are obtained based on Equation (26), respectively, i.e.,

Obviously, the method of constructing the ANFIS-GPR model is different from IATNN and MMCS, but the model coupling problem is transformed into an optimization problem, and the classical mathematical method – Lagrange multiplier method is introduced to solve it, realizing the organic combination of mathematical theory and neural network model.

## STUDY AREA AND DATA

### Case study site

^{2}(Figure 3). There are 16 rainfall stations and 1 hydrographic station named Ruzhou hydrographic station, and it has a typical continental monsoon region climate. The annual average temperature and evaporation in the period 1985–2016 are about 14 °C and 1,032 mm, respectively. The annual average precipitation in the period 1985–2016 is 744.7 mm, and rainfall is concentrated mainly in May, June, July, August, and September, comprising 62% of the average annual rainfall. The total length of the Beiru river is about 250 km, and Ruzhou hydrographic station is located at its outlet. There is no large water conservancy project upstream, and runoff is significantly affected by natural environmental factors, which provides a good physical environment foundation for this study.

### Data collection and preparation

Hydro-meteorological variables from 1/1985 to 12/2016, including the monthly precipitation (P), SMC, evaporation (E), and runoff (R), are gathered from Ruzhou hydrological station. According to the production and confluence theory, P, SMC, and E are the key factors affecting R. The potential runoff-driving factors can be established with 12 months as the maximum lead time, including 36 (12 × 3) variables. The details can be described as follows: The runoff-driving factors of P, SMC, and E are expressed as P(*t* − *i*), SMC(*t* − *i*), and E(*t* − *i*), respectively. Among them, the index signs *i* = 1, 2, ···, 12 indicate lead time. For example, P(*t* − *i*) indicates P 1 month ahead of R, and P(*t* − 12) indicates P 12 months ahead of *R*, so 12 alternative predictors can be formed based on P. If E and SMC are all counted together, there are 36 predictors.

## DEVELOPMENT OF THE MODELS FOR RUNOFF PREDICTION

### Selection of forecast factors

*t*− 1), SMC(

*t*− 1), P(

*t*− 2), P(

*t*− 12), E(

*t*− 7), E(

*t*− 8), E(

*t*− 9), E(

*t*− 6), E(

*t*− 2), E(

*t*− 3), and E(

*t*− 12) are selected in sequence. Depending on the theory of physical mechanism in hydrology, we know that the greater the P or SMC is, the bigger the R is, and the greater the E is, the smaller the R is. Therefore, P and SMC should positively correlate with R, while E and R have a negative relationship. That is, CC between P (or SMC) and R is positive, and CC between E and R is negative. Based on the above analysis, P(

*t*− 1), SMC(

*t*− 1), P(

*t*− 2), P (

*t*− 12), E(

*t*− 7), E(

*t*− 8), E(

*t*− 9), and E(

*t*− 6) are selected eventually as the suitable predictors in this article. Besides, there is a reasonable characteristic that P and SMC in the preceding month have a more significant influence on R than those in other months, which shows that CA is effective in selecting runoff prediction variables.

### Construction of the input schemes

The structure of the prediction models is determined by the input schemes because these variables can affect the debugging of model parameters. Therefore, the best input schemes are confirmed by trial and error. The specific approach is divided into the following three steps. Firstly, as shown in Table 1, seven input schemes are constructed based on the conclusions in section 4.1. ANFIS or GPR is used to evaluate the prediction performance of different input schemes to select better schemes. Secondly, PCA is used to conduct linear dimensional-reduction recombination (the cumulative contribution rate is set as 0.9) based on the better schemes to form the new input schemes. Finally, all the input schemes are compared to determine the final model input schemes.

Input scheme . | Input structure . | Number of variables . |
---|---|---|

M1 | P(t − 1), SMC(t − 1) | 2 |

M2 | P(t − 1), SMC(t − 1), P(t − 2) | 3 |

M3 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12) | 4 |

M4 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7) | 5 |

M5 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7), E(t − 8) | 6 |

M6 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7), E(t − 8), E(t − 9) | 7 |

M7 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7), E(t − 8), E(t − 9), E(t − 6) | 8 |

Input scheme . | Input structure . | Number of variables . |
---|---|---|

M1 | P(t − 1), SMC(t − 1) | 2 |

M2 | P(t − 1), SMC(t − 1), P(t − 2) | 3 |

M3 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12) | 4 |

M4 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7) | 5 |

M5 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7), E(t − 8) | 6 |

M6 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7), E(t − 8), E(t − 9) | 7 |

M7 | P(t − 1), SMC(t − 1), P(t − 2), P (t − 12), E(t − 7), E(t − 8), E(t − 9), E(t − 6) | 8 |

### Data used and normalization

As mentioned above in 3.2, the data series from 1985 to 2016 is divided into two parts, i.e., the training datasets (1985–2008) and the test datasets (2009–2016). The training datasets are applied to train the structure of the models, and the testing datasets serve to verify the performance of the models.

*x*,

*x*

_{max},

*x*

_{min}, and

*x*

_{nor}denote the raw data, the maximum value of data, the minimum value of data, and the normalized data. Finally, the model outputs need to be reversed to the raw value before the model performance evaluation.

### Performance indices

*et al.*2017):where

*x*,

_{i}*y*, , , and

_{i}*n*are the predicted value, measured value, the average value of

*x*, the average value of

_{i}*y*, and the length of the data.

_{i}### Parameter settings

The ANFIS, GPR, and ANFIS-GPR models are developed in a Matlab environment. The computations related to all models, including ANFIS, GPR, and ANFIS-GPR, are implemented in a computer with Inter core i5, 2.3-GHz CPU, and 16GB of RAM. The parameters of ANFIS are optimized by the backpropagation algorithm in the first and fourth layers. The Gaussian function is used as MFs. Fuzzy C-means clustering algorithm (FCM) is used to generate the primordial FIS, which can overcome the dimension disaster caused by excessive fuzzy rules in the ANFIS. The training and testing procedures of the ANFIS model are carried out from scratch for the mentioned datasets (Fattahi 2016). The maximum likelihood is calculated by directly maximizing the posterior distribution of the hyperparameters (*σ* and *l*) of GPR. The remaining key parameter settings are shown in Tables 2–4.

Parameter . | Description/value . |
---|---|

Number of clusters | 10 |

Partition matrix exponent | 2 |

Maximum number of iterations | 100 |

Minimum improvement | 0.00001 |

Parameter . | Description/value . |
---|---|

Number of clusters | 10 |

Partition matrix exponent | 2 |

Maximum number of iterations | 100 |

Minimum improvement | 0.00001 |

Parameter . | Description/value . |
---|---|

Maximum number of epochs | 300 |

Error goal | 0 |

Initial step size | 1.1 |

Step size decrease rate | 0.9 |

Step size increase rate | 1.1 |

Parameter . | Description/value . |
---|---|

Maximum number of epochs | 300 |

Error goal | 0 |

Initial step size | 1.1 |

Step size decrease rate | 0.9 |

Step size increase rate | 1.1 |

Parameter . | Description/value . |
---|---|

Kernel function | Mean excess function |

Covariance function | Rational quadratic covariance function |

Likelihood function | Gaussian likelihood function |

Parameter . | Description/value . |
---|---|

Kernel function | Mean excess function |

Covariance function | Rational quadratic covariance function |

Likelihood function | Gaussian likelihood function |

## RESULTS AND DISCUSSION

### Determination of the best input schemes

ANFIS and GPR are applied to evaluate the availability of seven different input schemes, i.e., M1, M2, M3, M4, M5, M6, and M7 are used as the inputs of ANFIS and GPR to predict R in turn, to obtain the performance indices (RMSE, NSC, and CORR) for comparative analysis. Due to the instability of ANFIS prediction, the average prediction statistics of ANFIS are obtained after running 10 times, while the stability of GPR prediction is good, so the prediction statistics of GPR are collected after running once. The RMSE, NSC, and CORR of ANFIS and GPR using M1, M2, M3, M4, M5, M6, and M7 in the testing period are shown in Table 5, and the best results are significantly signed. There is a significant trend between the performance of models and the input schemes. That is, the prediction accuracy of the models continues to increase as the number of variables in input schemes increases. In other words, models with more input variables usually show higher performance, such as M5, M6, and M7. Moreover, the results also reflect that the most satisfactory models’ prediction accuracy cannot be obtained by simply increasing the number of predictors. We can find that the prediction performance of models with M6 is better than that of models with M7, indicating that too many predictors cannot improve the prediction accuracy, but reduces the prediction performance of the models due to a large amount of redundant information.

Input scheme . | ANFIS . | GPR . | ||||
---|---|---|---|---|---|---|

RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . | |

M1 | 8.4766 | 0.4605 | 0.4578 | 8.1814 | 0.5954 | 0.5329 |

M2 | 8.4247 | 0.3217 | 0.3959 | 7.6445 | 0.5651 | 0.4744 |

M3 | 5.7224 | 0.5972 | 0.4665 | 7.1064 | 0.6006 | 0.4994 |

M4 | 6.8289 | 0.6474 | 0.5946 | 6.7069 | 0.7082 | 0.5468 |

M5 | 7.1524 | 0.5266 | 0.5529 | 6.6052 | 0.7284 | 0.5858 |

M6 | 7.4154 | 0.6784 | 0.5979 | 6.5725 | 0.7113 | 0.5918 |

M7 | 7.8122 | 0.4177 | 0.5285 | 6.4937 | 0.6705 | 0.5886 |

Input scheme . | ANFIS . | GPR . | ||||
---|---|---|---|---|---|---|

RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . | |

M1 | 8.4766 | 0.4605 | 0.4578 | 8.1814 | 0.5954 | 0.5329 |

M2 | 8.4247 | 0.3217 | 0.3959 | 7.6445 | 0.5651 | 0.4744 |

M3 | 5.7224 | 0.5972 | 0.4665 | 7.1064 | 0.6006 | 0.4994 |

M4 | 6.8289 | 0.6474 | 0.5946 | 6.7069 | 0.7082 | 0.5468 |

M5 | 7.1524 | 0.5266 | 0.5529 | 6.6052 | 0.7284 | 0.5858 |

M6 | 7.4154 | 0.6784 | 0.5979 | 6.5725 | 0.7113 | 0.5918 |

M7 | 7.8122 | 0.4177 | 0.5285 | 6.4937 | 0.6705 | 0.5886 |

The bold values represent where the optimal performance indicator values occur.

*t*− −1) and SMC(

*t*− 1) (between P(

*t*− 2) and E(

*t*− 9)) is the strongest, and the CC is 0.6611(−0.4956). Thus, it is worth using PCA for performing dimensionality reduction to construct the new input schemes (M5-PCA, M6-PCA, and M7-PCA), and the cumulative contribution rate is set as 0.9. Due to space reasons, the specific predictors of M5-PCA, M6-PCA, and M7-PCA are not listed here, and the main parameters of these schemes are shown in Table 6.

Input scheme . | Types of variables . | Number of variables . |
---|---|---|

M5-PCA | Recombinational data | 5 |

M6-PCA | Recombinational data | 6 |

M7-PCA | Recombinational data | 7 |

Input scheme . | Types of variables . | Number of variables . |
---|---|---|

M5-PCA | Recombinational data | 5 |

M6-PCA | Recombinational data | 6 |

M7-PCA | Recombinational data | 7 |

### Runoff prediction based on ANFIS, GPR, and ANFIS-GPR

Since the input schemes (i.e., M5-PCA, M6-PCA, and M7-PCA) have been proven to be effective in improving prediction accuracy, they served as the preferred input schemes for the models in this part. Because of the instability of ANFIS, the ANFIS-GPR model is also run 10 times in the same way. The average performances of ANFIS, GPR, and ANFIS-GPR for runoff prediction in the testing phase are shown in Table 7, respectively. It is seen that the prediction performances of ANFIS-GPR are better than ANFIS and GPR in terms of RMSE and NSC, and only CORR shows that ANFIS-GPR has a slightly lower performance than GPR. The results show that the proposed prediction technique is superior to other methods. Besides, we can find that the prediction performance of GPR is better than that of ANFIS, indicating that the probabilistic prediction model is more suitable than the coupled model for runoff prediction. Meanwhile, Table 8 shows the statistical performance of ANFIS and GPR models during the training period. It can be seen that GPR outperforms ANFIS in terms of RMSE, NSC, and CORR, demonstrating that GPR's robustness is better than that of ANFIS.

Input scheme . | M5-PCA . | M6-PCA . | M7-PCA . | ||||||
---|---|---|---|---|---|---|---|---|---|

Model . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . |

ANFIS | 1.0336 | 0.9966 | 0.9985 | 1.5985 | 0.9907 | 0.9964 | 1.2116 | 0.9952 | 0.9979 |

GPR | 0.8286 | 0.9978 | 0.9990 | 0.8843 | 0.9976 | 0.9989 | 0.8130 | 0.9979 | 0.9991 |

ANFIS-GPR | 0.7637 | 0.9991 | 0.9982 | 0.8453 | 0.9989 | 0.9978 | 0.7776 | 0.9991 | 0.9981 |

Input scheme . | M5-PCA . | M6-PCA . | M7-PCA . | ||||||
---|---|---|---|---|---|---|---|---|---|

Model . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . |

ANFIS | 1.0336 | 0.9966 | 0.9985 | 1.5985 | 0.9907 | 0.9964 | 1.2116 | 0.9952 | 0.9979 |

GPR | 0.8286 | 0.9978 | 0.9990 | 0.8843 | 0.9976 | 0.9989 | 0.8130 | 0.9979 | 0.9991 |

ANFIS-GPR | 0.7637 | 0.9991 | 0.9982 | 0.8453 | 0.9989 | 0.9978 | 0.7776 | 0.9991 | 0.9981 |

The bold values represent where the optimal performance indicator values occur.

Input scheme . | M5-PCA . | M6-PCA . | M7-PCA . | ||||||
---|---|---|---|---|---|---|---|---|---|

Model . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . |

ANFIS | 0.3593 | 0.9994 | 0.9997 | 0.4423 | 0.9991 | 0.9995 | 0.3829 | 0.9993 | 0.9997 |

GPR | 0.2861 | 0.9996 | 0.9997 | 0.3087 | 0.9993 | 0.9997 | 0.2788 | 0.9996 | 0.9998 |

Input scheme . | M5-PCA . | M6-PCA . | M7-PCA . | ||||||
---|---|---|---|---|---|---|---|---|---|

Model . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . | RMSE . | NSC . | CORR . |

ANFIS | 0.3593 | 0.9994 | 0.9997 | 0.4423 | 0.9991 | 0.9995 | 0.3829 | 0.9993 | 0.9997 |

GPR | 0.2861 | 0.9996 | 0.9997 | 0.3087 | 0.9993 | 0.9997 | 0.2788 | 0.9996 | 0.9998 |

The bold values represent where the optimal performance indicator values occur.

*R*

^{2}performance improvement of 0.073 and 0.03% present over ANFIS and GPR in the test phase, respectively, which indicates that the ANFIS-GPR is better at capturing the long-term characteristics of runoff series. Figure 8 shows the degree of deviation between the observed and the best-predicted values of different models, and it is difficult to compare the runoff prediction accuracy of each part intuitively. Therefore, we can make a detailed analysis of the prediction of the maximum runoff, because the flood peak discharge is very important for the guidance of water resources management and water conservancy project construction. The predicted values of ANFIS-GPR with M5-PCA and M6-PCA are closer to the observed, indicating that ANFIS-GPR can give a satisfactory result on the prediction accuracy of the key runoff.

## CONCLUSIONS

Improving runoff prediction accuracy is critical to water resources management and water conservancy project planning and design. A neoteric combinatorial model ANFIS-GPR for runoff prediction is proposed by merging ANFIS with GPR in this article. Firstly, the runoff forecasting outputs by ANFIS and GPR models are collected separately. Secondly, *w*_{1} and *w*_{2} are given to the different outputs from ANFIS and GPR. Finally, the issue of solving the weights is transformed into an optimization problem. The Lagrange multiplier method is used to calculate the weights, and the predicted values of ANFIS-GPR are obtained.

Meanwhile, the influence of different inputs on prediction accuracy is also studied. The results demonstrated that the relationship between prediction accuracy and model inputs is sensitive. It is not that the more input variables of the models are, the more prediction accuracy will be. That is, the prediction accuracy is closely related to the effective information provided by the input variables. The input schemes M5-PCA, M6-PCA, and M7-PCA based on the principal component extraction strategy are proven to be able to improve the prediction accuracy effectively. Thus, ANFIS, GPR, and ANFIS-GPR are used for runoff prediction with them (i.e., M5-PCA, M6-PCA, and M7-PCA). Results show that the ANFIS-GPR outperformed the ANFIS and GPR for runoff prediction in terms of RMSE, NSC, and *R*^{2}, and GPR outperformed the ANFIS-GPR in terms of CORR. In the single models, the performance of GPR is better than that of ANFIS.

Comparing the measures indicates that the coupling strategy of the ANFIS and GPR combination has more advantages in improving the accuracy of runoff prediction. Significantly, the ANFIS-GPR is more accurate than others in predicting the runoff peak values. Therefore, the runoff prediction project combining CA-PCA and ANFIS-GPR is the first choice in the study area.

## AVAILABILITY OF DATA AND MATERIALS

The data that support the funding of this study are available from the corresponding author upon reasonable request.

## FUNDING

This work is supported by the National Natural Science Foundation of China (No. 52069005 and 62163008) and Guizhou Provincial Science and Technology Projects (Guizhou Science Foundation-ZK[2021] General 295 and [2020]1Y248),and Science and Technology Special Funds of Guizhou Water Resources Department (KT202232), and Startup Project for High-level Talents of Guizhou Institute of Technology (XJGC20210425) and special thanks are given to the anonymous reviewers and editors for their constructive comments.

## AUTHOR CONTRIBUTIONS

Zhennan Liu conceptualized the whole article and developed the methodology. Jingnan Zhou wrote the original draft. Xianzhong Zeng conducted data curation. Xiaoyu Wang wrote the program. Weiguo Jiao wrote the review. Min Xu edited the article. Anjie Wu visualized the data.

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

**50**(