# Coupling a Bat Algorithm with XGBoost to Estimate Reference Evapotranspiration in the Arid and Semiarid Regions of China

Authors: Yixiu Han; Jianping Wu; Bingnian Zhai; Yanxin Pan; Guomin Huang; Lifeng Wu; Wenzhi Zeng
Journal: Advances in Meteorology
 (2019)

---

## Abstract

Accurate estimation of reference evapotranspiration (ETo) is key to agricultural irrigation scheduling and water resources management in arid and semiarid areas. This study evaluates the capability of coupling a Bat algorithm with the XGBoost method (i.e., the BAXGB model) for estimating monthly ETo in the arid and semiarid regions of China. Meteorological data from three stations (Datong, Yinchuan, and Taiyuan) during 1991–2015 were used to build the BAXGB model, the multivariate adaptive regression splines (MARS), and the gaussian process regression (GPR) model. Six input combinations with different sets of meteorological parameters were applied for model training and testing, which included mean air temperature (Tmean), maximum air temperature (Tmax), minimum air temperature (Tmin), wind speed (U), relative humidity (RH), and solar radiation (Rs) or extraterrestrial radiation (Ra, MJ m−2·d−1). The results indicated that BAXGB models (RMSE = 0.114–0.412 mm·d−1, MAE = 0.087–0.302 mm·d−1, and R2 = 0.937–0.996) were more accurate than either MARS (RMSE = 0.146–0.512 mm·d−1, MAE = 0.112–0.37 mm·d−1, and R2 = 0.935–0.994) or GPR (RMSE = 0.289–0.714 mm·d−1, MAE = 0.197–0.564 mm·d−1, and R2 = 0.817–0.980) model for estimating ETo. Findings of this study would be helpful for agricultural irrigation scheduling in the arid and semiarid regions and may be used as reference in other regions where accurate models for improving local water management are needed.


---

## Body

## 1. Introduction

There are two ways for water entering the atmosphere from the land surface. One is the direct evaporation of water from soils and plant surfaces, while the other one is transpiration, by which water from soils passes through plants and then evaporates on the leaves and diffuses to the surrounding air through stomata. Technically, it is difficult to distinguish between these two different pathways, which are collectively defined as crop evapotranspiration (ETc). Although ETc can be directly measured by techniques such as lysimeter or eddy covariance, these techniques cannot be applied widely due to their high costs in observation and maintenance. Therefore, ETc is generally estimated alternatively via multiplying the reference crop evapotranspiration (ETo) by the crop coefficient (Kc). According to Allen et al. [1], ETo represents the evapotranspiration from an actively growing virtual vegetated surface that is 0.12 m tall, with a completely shaded ground and adequate water supply, and without pests or diseases. At present, ETo has become a fundamental parameter in the fields of regional water resource management and irrigation scheduling. The most widely accepted ETo estimating method is the FAO-56 Penman–Monteith method (P-M), which has been proposed as a standard method for validating other methods because of its global applicability and simplicity (no need for parameter calibration) [1]. The method is function of basic meteorological parameters only, including temperature, solar radiation, wind speed, and relative humidity. However, even basic meteorological information in reality is often lacking, especially in developing countries, because of limited observational capacity and costs. For example, the lack of solar radiation records is a common problem in most parts of the world. Although there are over 2000 weather stations in China, less than 130 of them have observational records of radiation [2–4]. To solve the problem of lacking field observational data, various alternative methods, including empirical models, statistical models, artificial intelligence algorithms, and the satellite remote sensing technology, have been proposed and studied worldwide. Although these methods have different application scopes and characteristics, they all play a great role in improving the management of water resources.Artificial intelligence (AI) algorithms are becoming popular recently because of their high flexibility and precision and have been applied widely in fields of hydrology, agriculture, and environment [5–8], such as drought forecasting [9], solar radiation modeling [2, 10], temperature estimation [11], precipitation modeling [12], soil moisture simulation [13], and reference evapotranspiration estimation [14–17]. In the early years, artificial neural networks (ANNs) were used to simulate ETo, generating more accurate results than the empirical model [18]. Kisi [19] compared least square support vector machine (LSSVM), multivariate adaptive regression splines (MARS), and M5 model tree (M5T) for estimating ETo in Antalya and Isparta of Iran and found that the MARS algorithm was most accurate when local input data were not sufficient. Sanikhani et al. [20] compared six AI algorithms, including multilayer perceptron (MLP), generalized regression neural network (GRNN), radial basis neural networks (RBNN), integrated ANFIS with grid partitioning and subtractive clustering (ANFIS-GP and ANFIS-SC), and GEP, with data from the Antalya and Isparta stations of Iran. The results tended to be station specific, of which GEP and GRNN produced better performance than other algorithms in the Antalya station, while RBNN and ANFIS-SC were superior to other algorithms in the other station. Saggi and Jain [21] used four AI algorithms (including deep learning-multilayer perceptron (DL), generalized linear model (GLM), random forest (RF), and gradient boosting decision tree (GBDT)) to determine daily ETo in Punjab, India. The DL algorithm presented high performance for modeling daily ETo, with NSE = 0.95–0.98, RMSE = 0.192–0.269 mm·d−1, and R2 = 0.95–0.99. Landeras et al. [22] proposed an alternative strategy, using the combined data of local and its neighboring stations as input for ETo estimation with GEP and ANN, results of which showed that GEP had very good performance in all stations. Feng et al. [23] evaluated extreme learning machine (ELM), GRNN, and the Hargreaves model for the estimation of daily ETo using temperature data. They found that the Hargreaves was worse than AI algorithms, while ELM and GRNN performed best in the local scenario and the cross station scenario, respectively. Fan et al. [24] evaluated four tree-based algorithms (including M5T, RF, GBDT, and XGBoost) for estimating ETo in eight stations of China and found that XGBoost could yield comfortable accuracy but with much less time cost when compared with other algorithms. Moreover, although the gaussian process regressive (GPR) algorithm with a clear mathematical basis has been applied widely, only a few studies to date used it for crop ETo estimation [25].Parameters of AI models usually need to be calibrated, which is generally time-consuming, and requires a solid professional knowledge, especially when the inputs are complicated (e.g., input data from multiple sites or as multiple combinations) [26]. To overcome this problem, optimization with metaheuristic algorithms has been widely adopted by scholars recently, such as using the genetic algorithm (GA) to optimize ANN for ETo estimation in Iran [27] and hybridizing the firefly algorithm with ANFIS for estimating ETo in Burkina Faso [28]. The effects of metaheuristic optimization tend to be algorithm specific, with both positive and negative effects reported in various studies. For instance, Petković et al. [29] found that the radial basis function neural network (RBFNN) coupled with the particle swarm optimization (PSO) algorithm was more accurate than traditional artificial neural network algorithms.In recent years, the bat algorithm (BA) has received a lot of attention because of its strong searching ability and robustness, and the application of this algorithm has been discussed in many fields [28, 30, 31]. In addition, as a tree-based assemble algorithm, XGBoost has been proved superior to other simple tree algorithms such as RF and M5T in terms of accuracy and efficiency (i.e., less computational cost) for ETo estimation [24]. However, XGBoost contains many parameters, which would have impacts on its stability and efficiency during processing. Therefore, using the bat algorithm to optimize XGBoost is expected to improve the accuracy and stability of the single XGBoost model. To our knowledge, studies using evolutionary algorithms to optimize model parameters of XGBoost in hydrology, agriculture, and environment fields are scarce to date, especially for ETo estimation. Thus, the aims of this study were (1) to assess the capability of XGBoost optimized by BA in estimating monthly ETo with limit meteorological data as inputs and (2) to evaluate the performance of the BAXGB model for ETo estimation in comparison with GPR and MARS models in the arid and semiarid regions of China.


## 2. Materials and Methods

### 2.1. Study Area and Data

Monthly meteorological parameters during 1991–2015 from three weather stations in North China were selected in this study for model training and testing, including mean air temperature (Tmean, °C), maximum air temperature (Tmax, °C), minimum air temperature (Tmin, °C), wind speed (U, m s−1), solar radiation (Rs, MJ m−2·d−1) or extraterrestrial radiation (Ra, MJ m−2·d−1), and relative humidity (RH, %). For the selected weather stations (Figure 1), Datong and Yinchuan locate at the arid region of China, with a mean annual precipitation of 183 mm and 364 mm, respectively, while Taiyuan belongs to the semiarid region of China and has an average of 432 mm annual precipitation.Figure 1
The geographic information of the three weather stations in this study.All meteorological data used in this study were obtained from the National Meteorological Information Center of China (NMICC) with quality control. In addition, a further quality control was applied to the data of solar radiation following guidelines in [3] before processing. The missing data of all datasets were no more than 1% in total. The geographic information and meteorological statistics during the training and testing stages are listed in Table 1. It is noted that, compared with the training stage, wind speed of Datong and Yinchuan decreased by 10% and 29.2% in the testing stage, respectively, while solar radiation of Taiyuan increased by 16% during the testing stage. Six input combinations that consisted of different meteorological parameters were used for model training and testing in this study (Table 2). The map road of this study is described in Figure 2.Table 1
Statistics of the meteorological data of the three stations during the training and testing stages in the study.
Station/Meteorological data
Average
Maximum
Minimum
Skew
Kurtosis
CV
Training
Testing
Training
Testing
Training
Testing
Training
Testing
Training
Testing
Training
Testing
Datong
T
max
14.52
14.54
32.38
29.05
−6.85
−5.13
−0.24
−0.31
−1.34
−1.41
0.79
0.82
T
min
1.44
0.37
19.6
18.01
−18.93
−18.06
−0.14
−0.14
−1.36
−1.26
7.97
30.97
T
mean
7.98
7.45
25.99
23.53
−12.89
−11.6
−0.19
−0.23
−1.35
−1.34
4.38
15.9
R
s
14.73
15.72
24.79
24.75
5.55
6.06
−0.05
−0.09
−1.28
−0.93
0.37
0.46
U
2.07
1.86
3.37
2.89
0.96
1.12
0.46
0.8
0.33
0.97
0.2
0.3
RH
51.34
54.32
76.94
78.47
26.47
29.53
−0.11
−0.25
−0.61
1.34
0.22
0.36
ETo
2.83
2.76
6.46
5.72
0.44
0.5
0.17
0.21
−1.33
−1.21
0.58
0.69
Yinchuan
T
max
16.55
17.03
32.06
31.43
−5.45
−2.27
−0.29
−0.34
−1.29
−1.35
0.65
0.7
T
min
3.98
4.94
20.09
19.5
−17.27
−13.41
−0.16
−0.27
−1.33
−1.17
2.66
2.17
T
mean
10.26
10.98
26.08
25.47
−11.36
−7.84
−0.23
−0.3
−1.31
−1.26
1.66
1.43
R
s
15.93
15.4
25.85
24.05
6.4
7.02
0.04
−0.03
−1.28
−0.5
0.33
0.43
U
1.78
1.26
2.77
1.78
0.98
0.84
0.27
2.06
−0.23
5.59
0.21
0.54
RH
54.1
48.7
78.67
77.4
23.93
24.19
−0.46
0.03
−0.58
0.91
0.21
0.38
ETo
2.98
2.85
6.17
5.37
0.46
0.48
0.13
0.09
−1.39
−1.33
0.57
0.66
Taiyuan
T
max
17.62
17.79
32.75
31.22
−1.89
−0.16
−0.25
−0.35
−1.34
−1.37
0.57
0.63
T
min
4.9
5.37
21.36
20.18
−13.84
−11.21
−0.08
−0.2
−1.37
−1.26
2.07
1.91
T
mean
11.26
11.58
27.06
25.7
−7.87
−5.69
−0.17
−0.27
−1.36
−1.32
1.32
1.27
R
s
13.29
15.46
22.3
27.81
4.41
6.08
0.03
0.12
−1.21
−0.66
0.35
0.47
U
1.44
1.39
2.6
2.49
0.57
0.82
0.6
0.89
0.23
0.48
0.27
0.35
RH
57.05
55.12
85.74
79.57
28.79
35.52
−0.09
0.03
−0.88
2.01
0.22
0.35
ETo
2.61
2.83
5.51
5.6
0.49
0.63
0.2
0.2
−1.29
−1.3
0.53
0.65Table 2
Input combinations of different meteorological parameters used in this study.
Combinations
Models
Input combinations
GPR
MARS
BAXGB
1
GPR1
MARS1
BAXGB1
T
mean
⋅
R
a
2
GPR2
MARS2
BAXGB2
T
max
⋅
T
min
⋅
R
a
3
GPR3
MARS3
BAXGB3
T
max
⋅
T
min
⋅
R
s
4
GPR4
MARS4
BAXGB4
T
max
⋅
T
min
⋅
U
⋅
R
a
5
GPR5
MARS5
BAXGB5
T
max
⋅
T
min
⋅
RH
⋅
R
a
6
GPR6
MARS6
BAXGB6
T
max
⋅
T
min
⋅
RH
⋅
U
⋅
R
aFigure 2
The technology roadmap of the three algorithms in this study.


### 2.2. Multivariate Adaptive Regression Splines (MARS)

MARS is a nonlinear nonparametric regression technique introduced by Friedman [32], which uses the divide-and-conquer strategy for data processing. Specifically, data for model training are divided into subsets, each of which has its own function with at least two regression lines. The endpoints of a line segment are defined as nodes. During processing, specific assumptions regarding the basic functional relationship between input variables and outputs are not needed, but using the knotting method instead. For each input variable, knotting nodes can be made on any data point within the range of the variable, and the strategy of knotting is to search all possible knotting positions based on the adaptive regression algorithm and form a piecewise curve consequently. A knotting node marks the end of one data region and the beginning of another. Between the two adjoining knotting nodes, the regression of data forms a basis function (BF). Therefore, MARS generates BFs through the step-by-step searching, during which the node positions are automatically selected by an adaptive regression algorithm. If there are multiple input variables, it will search for all possible univariate node positions and interactions between all variables. However, these processes generally build an over fit model. To reduce the degree of overfitting, pruning is incorporated into MARS as the second phase of the model development. In summary, the MARS model is actually a weighed sum of BF, which makes the model more flexible. More details about MARS can be found in [32].


### 2.3. Gaussian Process Regression (GPR)

Gauss process regression is a framework based on the Bayesian theory with wide application in many fields. A GPR model can be defined as a probability distribution over functions [33], which is denoted as(1)fx=GPmx,kx,x′,mx=Efx,kx,x′=Efx−mxfx′−mx′.Since the Gaussian process is a normal distribution in the domain of functions, it can be interpreted from the perspective of space. If there is a predictor groupx consisting of d variables, it can be written as D=xk,fxk,k=1,2,3,…,n, where xk∈Rn is the input vector and fxk∈Rn is the corresponding output vector. In general, the Bayesian regression function can be expressed as(2)fx=φxTω,where ω is the weight vector, whose prior distribution is ω∼N0,Σp. For instance, if Eω=0, then EωωT=Σp.Since the expectation and variance can be expressed as(3)Efx=φxTEω=0,Efxfx′=φxTEωωTφx′=φxT∑pφx′,then we can see that when testing data and training data obey the same distribution, for a given testing vector x∗, the corresponding prediction vector f∗ also belongs to the Gaussian distribution. Suppose they form a joint Gaussian distribution:(4)ff∗∼Nμμ∗,Kx,xKx,x∗Kx,x∗TKx∗,x∗,where(5)K=kx1,x1kx1,x2⋯kx1,xnkx2,x1kx2,x2⋯kx1,xn⋮⋮⋱⋮kxn,x1kxn,x2⋯kxn,xn,k∗=kx∗,x1kx∗,x2⋯kx∗,xn,k∗∗=kx∗,x∗,in which kx,x is the correlation between x and itself, while Kx,x∗=Kx∗,xT is the N ∗ 1 order covariance matrix between the input and the output, briefly recorded as K=KT, and Kx∗,x∗ is the correlation of the output itself, which can be written in short as K∗∗.The Bayesian posterior probability is as follows:(6)Pf∗x∗,f=Nμ∗,∑∗,where(7)μ∗=μx∗+K∗TK−1f−μx,∑∗=K∗∗−K∗TK−1K∗,where μ∗ is the posterior probability means of f∗ and Σ∗ is the variance of posterior probability.When there is noise in the regression model,(8)y=fx+ε,where x is the input vector, f is the function value vector, and y is the observations with noise. It is generally considered that noise obeys a normal distribution and has a mean value equal to 0, with the variation σn expressed as ε∼N0,σn2. The prior distribution of observed values can be obtained as(9)y∼Nμx,Kx,x+σn2In,where In is the n-dimensional identity matrix. Then, the joint distribution of the observation set x,y and the prediction set x∗,y∗ can be written as follows:(10)yy∗∼Nμμ∗,Kx,xσn2InKx,x∗Kx,x∗TKx∗,x∗.Therefore, the posterior probability can be calculated as follows:(11)Py∗x∗,y=Nμ∗,∑∗,where(12)μ∗=μx∗+K∗TK+σn2In−1y−μx,∑∗=K∗∗−K∗TK+σn2In−1K∗.In general, to make the symbol concise, data need to be preprocessed, and the mean value can be set asμx=μx∗=0. Therefore, formula (11) can be written as(13)Py∗x∗,x,y∼NK∗TK+σn2In−1y−μx,K∗∗−K∗TK+σn2In−1K∗.For regression problems, the mean value is used as predictive values (y^), which can be written as(14)y^=K∗TK+σn2In−1y.Further simplification is as follows:(15)K+σn2In−1y=K+σn2In−1y1,y2,…,ynT=αT,y^=K∗TαT=∑i=1nKx∗,xiαi.According to the Mercer theorem,Kx∗,xi can be changed to the Gaussian kernel method:(16)Kx∗,xi=γexpx∗−xi22σk2,where γ and σk2 can be determined automatically by the Bayesian method.


### 2.4. XGBoost, Bat Algorithm, and BAXGB

XGBoost, proposed by Chen and Guestrin [34], is a scalable artificial intelligence algorithm for tree boosting. The original model of XGBoost is gradient boosting decision trees, which combines many “weak” learning models into a “strong” learner through an iterative method. As shown in Figure 3, during each iteration of XGBoost, the previous predictors are corrected with residuals. The algorithm can independently determine the types of loss functions used for model evaluation. To reduce the risk of overfitting, an additional regularization term is added to the model. The mean score of each tree is used as the predictive value for regression. For the m-th decision tree, its calculation formula can be expressed as(17)y^i=∑i=1mfmxi,fm∈W,where m is the number of trees, fm is a function in the functional space W, and W is the space of all decision trees.Figure 3
The structure of the XGBoost algorithm.The objective function at thet-th iteration can be written as(18)Θt=Φt+Ωt,Θt=∑i=1nΦyi,y^i+∑k=1tΩfm,where n is the n-th prediction and y^it can be given as(19)y^it=∑m=1tfmxi=y^it−1+ftxi.The regularization termΩfk for a decision tree is defined by Chen and Guestrin [34] as follows:(20)Ωfk=γT+12λ∑j=1Twj2,where γ is the complexity of each leaf, T is the number of leaves in a tree, λ is a parameter to scale the penalty, and w is the vector of scores on the leaves. Then, the first-order along with the second-order Taylor expansions are taken to the loss function in XGBoost. For example, if the loss function is the mean square error (MSE), then the objective function can finally be derived as(21)Θt≈∑i=1ngiwqxi+12hiwqxi2+γT+12∑j=1Twj2,where q⋅ is a function that can assign data points to corresponding leaves and gi and hi are the first and second derivatives of MSE loss function, respectively. In formula (21), the loss function is determined by the sum of loss values for each data sample. Since each data sample corresponds to only one leaf node, the loss function can also be expressed as the sum of loss values for each leaf node, which is(22)Θt≈∑j=1T∑i∈Ijgiwj+12∑i∈Ijhi+λwj2+γT,where Ij represents all data samples in the leaf node j. Hence, optimization of the objective function is now transformed into finding the minimum value of a quadratic function. In other words, the change of model performance caused by a certain node split in the decision tree now can be evaluated based on the objective function. If the model performance of the decision tree is improved after the node split, then the change will be adopted; otherwise, the split will be stopped. The structure of the XGBoost algorithm can be seen in Figure 3. Besides, when optimizing the objective function, a predictive classifier can be trained against overfitting because of the regularization terms. In this study, three parameters were optimized, including the number of trees (nrounds), the ratio of subdatasets to all data in training the model (subset), and the minimum sum of instance weight (hessian) needed in a child (min_child_weight).The bat algorithm (BA) is a metaheuristic algorithm developed by Yang [35], which mirrors the echolocation behavior of microbats. Microbats emit sound pulses to search for prey, and each of them has its own echolocation characteristics (i.e., pulse emission rates and loudness), the feature of which can be used for developing various BA. The main processes of constructing a BA can be summarized as follows:(1)
A population of virtual bats is first generated for simulations. For each virtual bat, assume that it flies randomly with a velocityvi and a frequency fi at position xi.(2)
From the first iteration to the maximum iteration, positions and velocities of bats at timet are updated by(23)fi=fmin+fmax−fminβ,vit=vit−1+xit−x∗fi,xit=xit−1+vit,where β∈0,1 is a random vector from a normal distribution, xit and vit are the updated positions and velocities of bats at time t, respectively, and x∗ is the current best position (solution) a bat is located after comparing all the solutions among all the bats within the population.(3)
Once the current best solution is selected, the bats will update their positions based on a random walk:(24)xnew=xold+rand⋅At,where rand∈−1,1 is a random number and At is the mean loudness of all bats at time t.(4)
As the iteration proceed, the rate of pulse emissionri and the loudness Ai of each bat within the population will be updated according to the following rules:(25)Ait+1=αAit,rit+1=ri01−exp−γt,where α and γ are both constants. The emission rates and loudness of all bats will be updated until the best solutions are found, and the iteration will be terminated by then.In this study, a hybrid model integrating BA and XGBoost (i.e., BAXGB) was developed for ETo estimation. As stated previously, the performance of XGBoost models is largely dependent on the appropriate choice of parameters. In BAXGB, the parameter choosing for XGBoost was determined and optimized by BA. This paper optimizes three important parameters of the XGBoost algorithm, which are the number of trees (nrounds), the minimum sum of instance weight needed in a child (Min_child_weight), and the study rate (eta). The parameter values are shown in Table3.Table 3
Parameter values of the Bat-XGBoost algorithm.
Parameter
Value
Population
50
Max generation
200
A
0.98
R
0.98
nrounds
[20, 2000]
Min_child_weight
[1, 15]
Eta
[0.01, 0.3]


### 2.5. Model Performance Evaluation

Five statistical indicators were used for the performance evaluation of GPR, MARS, and BAXGB models in this study, which were as follows.(1)
Coefficient of determination:(26)R2=∑i=1nOi,m−O¯i,mOi,e−O¯i,e2∑i=1nOi,m−O¯i,m2∑i=1nOi,e−O¯i,e2.(2)
Root-mean-square error (RMSE):(27)RMSE=1n∑i=1nOi,m−Oi,e2.(3)
Mean absolute error (MAE):(28)MAE=1n∑i=1nOi,m−Oi,e.(4)
The maximum 10% RMSE (RMSEmax 10%):(29)RMSEmax10%=1z∑i=1zEmax10%2.(5)
The maximum 10% MAE (MAEmax 10%):(30)MAEmax10%=1z∑i=1zEmax10%,where n is the total number of test data, Oe and Om are the predicted ETo values by artificial intelligence methods and the FAO-56 P-M ETo values, respectively, O¯m is the mean P-M ETo value, and Emax10% is the maximum 10% absolute error of the test data and the FAO-56 P-M ETo.The value ofR2 is between 0 and 1, and the value closer to 1 indicates better regression fitting between predicted ETo and P-M ETo and hence the better model performance. On the contrary, if the value of RMSE or MAE is closer to 0, the model performs better. RMSE is generally useful when model errors follow a normal distribution, whereas MAE suits better for models with a uniform error distribution [24]. In addition, considering people’s different acceptability of the ETo error in different months, as well as the uncertainty caused by the water resource management process, it is necessary to evaluate not only the overall performance of the model but also whether the performance of the model meets the expectation in extreme cases. Therefore, two new statistical indicators using the maximum error point of 10% were adopted in this study to evaluate the model performance in extreme cases.


## 2.1. Study Area and Data

Monthly meteorological parameters during 1991–2015 from three weather stations in North China were selected in this study for model training and testing, including mean air temperature (Tmean, °C), maximum air temperature (Tmax, °C), minimum air temperature (Tmin, °C), wind speed (U, m s−1), solar radiation (Rs, MJ m−2·d−1) or extraterrestrial radiation (Ra, MJ m−2·d−1), and relative humidity (RH, %). For the selected weather stations (Figure 1), Datong and Yinchuan locate at the arid region of China, with a mean annual precipitation of 183 mm and 364 mm, respectively, while Taiyuan belongs to the semiarid region of China and has an average of 432 mm annual precipitation.Figure 1
The geographic information of the three weather stations in this study.All meteorological data used in this study were obtained from the National Meteorological Information Center of China (NMICC) with quality control. In addition, a further quality control was applied to the data of solar radiation following guidelines in [3] before processing. The missing data of all datasets were no more than 1% in total. The geographic information and meteorological statistics during the training and testing stages are listed in Table 1. It is noted that, compared with the training stage, wind speed of Datong and Yinchuan decreased by 10% and 29.2% in the testing stage, respectively, while solar radiation of Taiyuan increased by 16% during the testing stage. Six input combinations that consisted of different meteorological parameters were used for model training and testing in this study (Table 2). The map road of this study is described in Figure 2.Table 1
Statistics of the meteorological data of the three stations during the training and testing stages in the study.
Station/Meteorological data
Average
Maximum
Minimum
Skew
Kurtosis
CV
Training
Testing
Training
Testing
Training
Testing
Training
Testing
Training
Testing
Training
Testing
Datong
T
max
14.52
14.54
32.38
29.05
−6.85
−5.13
−0.24
−0.31
−1.34
−1.41
0.79
0.82
T
min
1.44
0.37
19.6
18.01
−18.93
−18.06
−0.14
−0.14
−1.36
−1.26
7.97
30.97
T
mean
7.98
7.45
25.99
23.53
−12.89
−11.6
−0.19
−0.23
−1.35
−1.34
4.38
15.9
R
s
14.73
15.72
24.79
24.75
5.55
6.06
−0.05
−0.09
−1.28
−0.93
0.37
0.46
U
2.07
1.86
3.37
2.89
0.96
1.12
0.46
0.8
0.33
0.97
0.2
0.3
RH
51.34
54.32
76.94
78.47
26.47
29.53
−0.11
−0.25
−0.61
1.34
0.22
0.36
ETo
2.83
2.76
6.46
5.72
0.44
0.5
0.17
0.21
−1.33
−1.21
0.58
0.69
Yinchuan
T
max
16.55
17.03
32.06
31.43
−5.45
−2.27
−0.29
−0.34
−1.29
−1.35
0.65
0.7
T
min
3.98
4.94
20.09
19.5
−17.27
−13.41
−0.16
−0.27
−1.33
−1.17
2.66
2.17
T
mean
10.26
10.98
26.08
25.47
−11.36
−7.84
−0.23
−0.3
−1.31
−1.26
1.66
1.43
R
s
15.93
15.4
25.85
24.05
6.4
7.02
0.04
−0.03
−1.28
−0.5
0.33
0.43
U
1.78
1.26
2.77
1.78
0.98
0.84
0.27
2.06
−0.23
5.59
0.21
0.54
RH
54.1
48.7
78.67
77.4
23.93
24.19
−0.46
0.03
−0.58
0.91
0.21
0.38
ETo
2.98
2.85
6.17
5.37
0.46
0.48
0.13
0.09
−1.39
−1.33
0.57
0.66
Taiyuan
T
max
17.62
17.79
32.75
31.22
−1.89
−0.16
−0.25
−0.35
−1.34
−1.37
0.57
0.63
T
min
4.9
5.37
21.36
20.18
−13.84
−11.21
−0.08
−0.2
−1.37
−1.26
2.07
1.91
T
mean
11.26
11.58
27.06
25.7
−7.87
−5.69
−0.17
−0.27
−1.36
−1.32
1.32
1.27
R
s
13.29
15.46
22.3
27.81
4.41
6.08
0.03
0.12
−1.21
−0.66
0.35
0.47
U
1.44
1.39
2.6
2.49
0.57
0.82
0.6
0.89
0.23
0.48
0.27
0.35
RH
57.05
55.12
85.74
79.57
28.79
35.52
−0.09
0.03
−0.88
2.01
0.22
0.35
ETo
2.61
2.83
5.51
5.6
0.49
0.63
0.2
0.2
−1.29
−1.3
0.53
0.65Table 2
Input combinations of different meteorological parameters used in this study.
Combinations
Models
Input combinations
GPR
MARS
BAXGB
1
GPR1
MARS1
BAXGB1
T
mean
⋅
R
a
2
GPR2
MARS2
BAXGB2
T
max
⋅
T
min
⋅
R
a
3
GPR3
MARS3
BAXGB3
T
max
⋅
T
min
⋅
R
s
4
GPR4
MARS4
BAXGB4
T
max
⋅
T
min
⋅
U
⋅
R
a
5
GPR5
MARS5
BAXGB5
T
max
⋅
T
min
⋅
RH
⋅
R
a
6
GPR6
MARS6
BAXGB6
T
max
⋅
T
min
⋅
RH
⋅
U
⋅
R
aFigure 2
The technology roadmap of the three algorithms in this study.


## 2.2. Multivariate Adaptive Regression Splines (MARS)

MARS is a nonlinear nonparametric regression technique introduced by Friedman [32], which uses the divide-and-conquer strategy for data processing. Specifically, data for model training are divided into subsets, each of which has its own function with at least two regression lines. The endpoints of a line segment are defined as nodes. During processing, specific assumptions regarding the basic functional relationship between input variables and outputs are not needed, but using the knotting method instead. For each input variable, knotting nodes can be made on any data point within the range of the variable, and the strategy of knotting is to search all possible knotting positions based on the adaptive regression algorithm and form a piecewise curve consequently. A knotting node marks the end of one data region and the beginning of another. Between the two adjoining knotting nodes, the regression of data forms a basis function (BF). Therefore, MARS generates BFs through the step-by-step searching, during which the node positions are automatically selected by an adaptive regression algorithm. If there are multiple input variables, it will search for all possible univariate node positions and interactions between all variables. However, these processes generally build an over fit model. To reduce the degree of overfitting, pruning is incorporated into MARS as the second phase of the model development. In summary, the MARS model is actually a weighed sum of BF, which makes the model more flexible. More details about MARS can be found in [32].


## 2.3. Gaussian Process Regression (GPR)

Gauss process regression is a framework based on the Bayesian theory with wide application in many fields. A GPR model can be defined as a probability distribution over functions [33], which is denoted as(1)fx=GPmx,kx,x′,mx=Efx,kx,x′=Efx−mxfx′−mx′.Since the Gaussian process is a normal distribution in the domain of functions, it can be interpreted from the perspective of space. If there is a predictor groupx consisting of d variables, it can be written as D=xk,fxk,k=1,2,3,…,n, where xk∈Rn is the input vector and fxk∈Rn is the corresponding output vector. In general, the Bayesian regression function can be expressed as(2)fx=φxTω,where ω is the weight vector, whose prior distribution is ω∼N0,Σp. For instance, if Eω=0, then EωωT=Σp.Since the expectation and variance can be expressed as(3)Efx=φxTEω=0,Efxfx′=φxTEωωTφx′=φxT∑pφx′,then we can see that when testing data and training data obey the same distribution, for a given testing vector x∗, the corresponding prediction vector f∗ also belongs to the Gaussian distribution. Suppose they form a joint Gaussian distribution:(4)ff∗∼Nμμ∗,Kx,xKx,x∗Kx,x∗TKx∗,x∗,where(5)K=kx1,x1kx1,x2⋯kx1,xnkx2,x1kx2,x2⋯kx1,xn⋮⋮⋱⋮kxn,x1kxn,x2⋯kxn,xn,k∗=kx∗,x1kx∗,x2⋯kx∗,xn,k∗∗=kx∗,x∗,in which kx,x is the correlation between x and itself, while Kx,x∗=Kx∗,xT is the N ∗ 1 order covariance matrix between the input and the output, briefly recorded as K=KT, and Kx∗,x∗ is the correlation of the output itself, which can be written in short as K∗∗.The Bayesian posterior probability is as follows:(6)Pf∗x∗,f=Nμ∗,∑∗,where(7)μ∗=μx∗+K∗TK−1f−μx,∑∗=K∗∗−K∗TK−1K∗,where μ∗ is the posterior probability means of f∗ and Σ∗ is the variance of posterior probability.When there is noise in the regression model,(8)y=fx+ε,where x is the input vector, f is the function value vector, and y is the observations with noise. It is generally considered that noise obeys a normal distribution and has a mean value equal to 0, with the variation σn expressed as ε∼N0,σn2. The prior distribution of observed values can be obtained as(9)y∼Nμx,Kx,x+σn2In,where In is the n-dimensional identity matrix. Then, the joint distribution of the observation set x,y and the prediction set x∗,y∗ can be written as follows:(10)yy∗∼Nμμ∗,Kx,xσn2InKx,x∗Kx,x∗TKx∗,x∗.Therefore, the posterior probability can be calculated as follows:(11)Py∗x∗,y=Nμ∗,∑∗,where(12)μ∗=μx∗+K∗TK+σn2In−1y−μx,∑∗=K∗∗−K∗TK+σn2In−1K∗.In general, to make the symbol concise, data need to be preprocessed, and the mean value can be set asμx=μx∗=0. Therefore, formula (11) can be written as(13)Py∗x∗,x,y∼NK∗TK+σn2In−1y−μx,K∗∗−K∗TK+σn2In−1K∗.For regression problems, the mean value is used as predictive values (y^), which can be written as(14)y^=K∗TK+σn2In−1y.Further simplification is as follows:(15)K+σn2In−1y=K+σn2In−1y1,y2,…,ynT=αT,y^=K∗TαT=∑i=1nKx∗,xiαi.According to the Mercer theorem,Kx∗,xi can be changed to the Gaussian kernel method:(16)Kx∗,xi=γexpx∗−xi22σk2,where γ and σk2 can be determined automatically by the Bayesian method.


## 2.4. XGBoost, Bat Algorithm, and BAXGB

XGBoost, proposed by Chen and Guestrin [34], is a scalable artificial intelligence algorithm for tree boosting. The original model of XGBoost is gradient boosting decision trees, which combines many “weak” learning models into a “strong” learner through an iterative method. As shown in Figure 3, during each iteration of XGBoost, the previous predictors are corrected with residuals. The algorithm can independently determine the types of loss functions used for model evaluation. To reduce the risk of overfitting, an additional regularization term is added to the model. The mean score of each tree is used as the predictive value for regression. For the m-th decision tree, its calculation formula can be expressed as(17)y^i=∑i=1mfmxi,fm∈W,where m is the number of trees, fm is a function in the functional space W, and W is the space of all decision trees.Figure 3
The structure of the XGBoost algorithm.The objective function at thet-th iteration can be written as(18)Θt=Φt+Ωt,Θt=∑i=1nΦyi,y^i+∑k=1tΩfm,where n is the n-th prediction and y^it can be given as(19)y^it=∑m=1tfmxi=y^it−1+ftxi.The regularization termΩfk for a decision tree is defined by Chen and Guestrin [34] as follows:(20)Ωfk=γT+12λ∑j=1Twj2,where γ is the complexity of each leaf, T is the number of leaves in a tree, λ is a parameter to scale the penalty, and w is the vector of scores on the leaves. Then, the first-order along with the second-order Taylor expansions are taken to the loss function in XGBoost. For example, if the loss function is the mean square error (MSE), then the objective function can finally be derived as(21)Θt≈∑i=1ngiwqxi+12hiwqxi2+γT+12∑j=1Twj2,where q⋅ is a function that can assign data points to corresponding leaves and gi and hi are the first and second derivatives of MSE loss function, respectively. In formula (21), the loss function is determined by the sum of loss values for each data sample. Since each data sample corresponds to only one leaf node, the loss function can also be expressed as the sum of loss values for each leaf node, which is(22)Θt≈∑j=1T∑i∈Ijgiwj+12∑i∈Ijhi+λwj2+γT,where Ij represents all data samples in the leaf node j. Hence, optimization of the objective function is now transformed into finding the minimum value of a quadratic function. In other words, the change of model performance caused by a certain node split in the decision tree now can be evaluated based on the objective function. If the model performance of the decision tree is improved after the node split, then the change will be adopted; otherwise, the split will be stopped. The structure of the XGBoost algorithm can be seen in Figure 3. Besides, when optimizing the objective function, a predictive classifier can be trained against overfitting because of the regularization terms. In this study, three parameters were optimized, including the number of trees (nrounds), the ratio of subdatasets to all data in training the model (subset), and the minimum sum of instance weight (hessian) needed in a child (min_child_weight).The bat algorithm (BA) is a metaheuristic algorithm developed by Yang [35], which mirrors the echolocation behavior of microbats. Microbats emit sound pulses to search for prey, and each of them has its own echolocation characteristics (i.e., pulse emission rates and loudness), the feature of which can be used for developing various BA. The main processes of constructing a BA can be summarized as follows:(1)
A population of virtual bats is first generated for simulations. For each virtual bat, assume that it flies randomly with a velocityvi and a frequency fi at position xi.(2)
From the first iteration to the maximum iteration, positions and velocities of bats at timet are updated by(23)fi=fmin+fmax−fminβ,vit=vit−1+xit−x∗fi,xit=xit−1+vit,where β∈0,1 is a random vector from a normal distribution, xit and vit are the updated positions and velocities of bats at time t, respectively, and x∗ is the current best position (solution) a bat is located after comparing all the solutions among all the bats within the population.(3)
Once the current best solution is selected, the bats will update their positions based on a random walk:(24)xnew=xold+rand⋅At,where rand∈−1,1 is a random number and At is the mean loudness of all bats at time t.(4)
As the iteration proceed, the rate of pulse emissionri and the loudness Ai of each bat within the population will be updated according to the following rules:(25)Ait+1=αAit,rit+1=ri01−exp−γt,where α and γ are both constants. The emission rates and loudness of all bats will be updated until the best solutions are found, and the iteration will be terminated by then.In this study, a hybrid model integrating BA and XGBoost (i.e., BAXGB) was developed for ETo estimation. As stated previously, the performance of XGBoost models is largely dependent on the appropriate choice of parameters. In BAXGB, the parameter choosing for XGBoost was determined and optimized by BA. This paper optimizes three important parameters of the XGBoost algorithm, which are the number of trees (nrounds), the minimum sum of instance weight needed in a child (Min_child_weight), and the study rate (eta). The parameter values are shown in Table3.Table 3
Parameter values of the Bat-XGBoost algorithm.
Parameter
Value
Population
50
Max generation
200
A
0.98
R
0.98
nrounds
[20, 2000]
Min_child_weight
[1, 15]
Eta
[0.01, 0.3]


## 2.5. Model Performance Evaluation

Five statistical indicators were used for the performance evaluation of GPR, MARS, and BAXGB models in this study, which were as follows.(1)
Coefficient of determination:(26)R2=∑i=1nOi,m−O¯i,mOi,e−O¯i,e2∑i=1nOi,m−O¯i,m2∑i=1nOi,e−O¯i,e2.(2)
Root-mean-square error (RMSE):(27)RMSE=1n∑i=1nOi,m−Oi,e2.(3)
Mean absolute error (MAE):(28)MAE=1n∑i=1nOi,m−Oi,e.(4)
The maximum 10% RMSE (RMSEmax 10%):(29)RMSEmax10%=1z∑i=1zEmax10%2.(5)
The maximum 10% MAE (MAEmax 10%):(30)MAEmax10%=1z∑i=1zEmax10%,where n is the total number of test data, Oe and Om are the predicted ETo values by artificial intelligence methods and the FAO-56 P-M ETo values, respectively, O¯m is the mean P-M ETo value, and Emax10% is the maximum 10% absolute error of the test data and the FAO-56 P-M ETo.The value ofR2 is between 0 and 1, and the value closer to 1 indicates better regression fitting between predicted ETo and P-M ETo and hence the better model performance. On the contrary, if the value of RMSE or MAE is closer to 0, the model performs better. RMSE is generally useful when model errors follow a normal distribution, whereas MAE suits better for models with a uniform error distribution [24]. In addition, considering people’s different acceptability of the ETo error in different months, as well as the uncertainty caused by the water resource management process, it is necessary to evaluate not only the overall performance of the model but also whether the performance of the model meets the expectation in extreme cases. Therefore, two new statistical indicators using the maximum error point of 10% were adopted in this study to evaluate the model performance in extreme cases.


## 3. Results and Discussion

### 3.1. Model Performance in the Datong Station

Performance of the three artificial intelligence models used in this study (i.e., GPR, MARS, and BAXGB) during the training and testing stages for the Datong station is listed in Table4. Both air temperature and extraterrestrial radiation were used for the first two input combinations. The difference was that the input combination 1 used the monthly average temperature, while the input combination 2 used the monthly maximum and minimum temperatures. Under these two input combinations, the accuracies of the three types of models varied. For GPR and BAXGB models, input combination 2 had higher accuracy than input combination 1. However, MARS models yielded an opposite result. Input combinations 3, 4, and 5 represented the effects of solar radiation, wind speed, and relative humidity or temperature data, respectively. Considering that the extraterrestrial radiation is only related to the latitude of the Earth and does not belong to the meteorological parameters to be observed, this information was added to the combinations 4 and 5. The importance of meteorological factors given by different models was different. MARS and BAXGB models showed that relative humidity was the most important, followed by wind speed and solar radiation. However, GPR model results showed that wind speed was the most critical, followed by relative humidity and solar radiation. Although there was difference to some degree in the importance between different meteorological parameters, it was still clear that relative humidity and wind speed were more important than solar radiation for the Datong station. Input combination 6 was a combination of temperature, wind speed, and relative humidity. For the three artificial intelligence models, results of input combination 6 were better than other combinations that with relative humidity or wind speed alone. Comparisons between the three artificial intelligence models under the same input combination indicated that BAXGB (RMSE = 0.116–0.438 mm·d−1 and MAE = 0.087–0.311 mm·d−1) always performed better than MARS (RMSE = 0.153–0.512 mm·d−1 and MAE = 0.221– 0.559 mm·d−1) and GPR (RMSE = 0.297–0.695 mm·d−1 and MAE = 0.120–0.379 mm·d−1) under each input combination. In addition, the MARS model generally performed better than the GPR model, except in combinations 2 and 3. Across algorithms and input combinations, values of RMSEmax10% and MAEmax10% during the testing stage were 2–2.5 times of their corresponding values of RMSE and MAE, suggesting that model performance in extreme cases would be worse than normal conditions. Although the overall performance of the BAXGB model was good, the maximum error point of 10% with input combination 3 (i.e., BAXGB3) exceeded 1.0 mm·d−1. This may be due to the difference in distribution of solar radiation data between the training and the testing stage.Table 4
Results of the statistical indicators for all models in the Datong station during the training and testing stages.
Model
Training
Testing
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
GPR1
0.841
0.66
0.519
1.392
1.373
0.832
0.695
0.559
1.333
1.328
GPR2
0.978
0.258
0.196
0.586
0.563
0.972
0.372
0.302
0.744
0.738
GPR3
0.98
0.242
0.181
0.565
0.547
0.956
0.444
0.35
0.898
0.891
GPR4
0.974
0.279
0.202
0.673
0.645
0.978
0.297
0.221
0.668
0.664
GPR5
0.988
0.197
0.138
0.483
0.46
0.951
0.398
0.251
1.089
1.019
GPR6
0.986
0.209
0.139
0.548
0.517
0.958
0.37
0.25
0.919
0.877
MARS1
0.964
0.315
0.236
0.723
0.71
0.97
0.292
0.227
0.649
0.631
MARS2
0.986
0.199
0.15
0.465
0.453
0.962
0.421
0.337
0.821
0.812
MARS3
0.984
0.209
0.157
0.481
0.468
0.953
0.512
0.379
1.158
1.148
MARS4
0.99
0.162
0.121
0.371
0.363
0.98
0.271
0.216
0.559
0.552
MARS5
0.992
0.154
0.11
0.377
0.362
0.992
0.153
0.12
0.321
0.32
MARS6
0.994
0.122
0.087
0.3
0.292
0.992
0.189
0.121
0.509
0.483
BAXGB1
0.978
0.238
0.179
0.541
0.531
0.972
0.291
0.23
0.592
0.587
BAXGB2
0.988
0.209
0.148
0.518
0.505
0.974
0.284
0.227
0.576
0.568
BAXGB3
0.998
0.076
0.057
0.175
0.168
0.955
0.438
0.311
1.067
1.057
BAXGB4
0.996
0.092
0.069
0.214
0.207
0.988
0.191
0.153
0.403
0.398
BAXGB5
0.996
0.108
0.076
0.266
0.255
0.992
0.146
0.112
0.337
0.328
BAXGB6
1
0.03
0.022
0.068
0.066
0.996
0.116
0.087
0.272
0.258
Note. Bold font indicates the best model.Figure4 presents scatter diagrams between the estimated ETo by the three artificial intelligence algorithms and the P-M ETo (mm·d−1) during the testing stage for the Datong station. In Figure 4, the R2 value ranged at 0.832–0.978, 0.953–0.992, and 0.955–0.996 for GPR, MARS, and BAXGB for all input combinations, respectively. The slopes of the GPR models were all less than 0.95 except with the input combination 3, which led to a problem of underestimation. These results clearly indicated that, for the Datong station, BAXGB was the most superior, while GPR was the worst among the AI algorithms, in terms of the performance in ETo estimation.Figure 4
Scatter plots between the estimated ETo by the three algorithms and the FAO-56 P-M ETo in the testing stage for the Datong station.
(a)
(b)
(c)
(d)
(e)
(f)Figure5 presents comparisons between the estimated ETo by AI algorithms and the P-M ETo (mm·d−1) during the testing stage for the Datong station. P-M ETo peaked in May with a value over 5 mm·d−1 while had the lowest value in December (0.6 mm·d−1). For the GPR model, the simulation effect of input combination 1 did not catch the trend well, resulting in a deviation more than 1 mm·d−1, presented as underestimation in May and overestimation in September, respectively. For the MARS model, input combinations 2 and 3 had a problem of overall overestimation, especially for months from August to October where the overestimation was more than 0.4 mm·d−1. For the BAXGB model, all input combination simulated an acceptable trend, except for the input combination 3 where an overestimation (more than 0.4 mm·d−1) was detected from July to October. Overall, among the AI methods in this study, the BAXGB model had the smallest deviation in all the 12 months, while GPR models performed worst.Figure 5
Relative errors by the different algorithms at each month for the Datong station.
(a)
(b)
(c)
(d)
(e)
(f)To further evaluate the performance of different AI models, Taylor diagrams were used to analyze the performance of models under different combinations, and the results of the Datong station can be seen in Figure6. It is not hard to see that combination 4 performs best in the GPR models, while combination 5 performs best in the MARS models. However, the combination 6 is superior to others in the BAXGB models. In addition, the point of the BAXGB model shows much closer than the other two models. This is consistent with the previous statistical indicators and scatter plot results.Figure 6
The Taylor diagram of (a) GPR, (b) MARS, and (c) BAXGB at the Datong station.
(a)
(b)
(c)


### 3.2. Model Performance in the Yinchuan Station

The performance of the three artificial intelligence models in the training and testing stages for the Yinchuan station is compared in Table5. In the first two input combinations, combination 2 had higher accuracy than combination 1 for all the three algorithms used in Yinchuan, which was different from Datong. However, the rank of importance of meteorological factors in different models showed similar trends with the Datong station. In addition, the combination including both wind speed and relative humidity (i.e., input combination 6) was still better than either factor alone for MARS and BAXGB models. It can be seen from the comparison of the three artificial intelligence models under the same input combination that BAXGB (RMSE = 0.175–0.261  mm·d−1 and MAE = 0.123–0.204 mm·d−1) performed better than MARS (RMSE = 0.172–0.340 mm·d−1 and MAE = 0.108– 266 mm·d−1) and GPR (RMSE = 0.289–0.663 mm·d−1 and MAE = 0.195–0.511 mm·d−1) under each combination. In general, the BAXGB model was better than the MARS model, and the GPR model was the worst, except for some individual combinations. Similar to the Datong station, RMSEmax10% and MAEmax10% during the testing stage were more than twice as high as the corresponding RMSE and MAE in the Yinchuan station. However, compared with Datong, values of RMSEmax10% and MAEmax10% were generally smaller in Yinchuan.Table 5
Results of the statistical indicators for all models in the Yinchuan station during the training and testing stages.
Model
Training
Testing
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
GPR1
0.856
0.653
0.514
1.371
1.351
0.87
0.663
0.511
1.424
1.401
GPR2
0.982
0.248
0.191
0.568
0.548
0.974
0.334
0.259
0.73
0.719
GPR3
0.984
0.23
0.176
0.53
0.518
0.974
0.289
0.197
0.751
0.696
GPR4
0.98
0.257
0.196
0.588
0.574
0.966
0.343
0.262
0.752
0.736
GPR5
0.984
0.23
0.169
0.536
0.519
0.972
0.461
0.378
0.941
0.907
GPR6
0.986
0.208
0.147
0.507
0.491
0.97
0.384
0.285
0.888
0.878
MARS1
0.968
0.3
0.229
0.677
0.666
0.974
0.34
0.262
0.728
0.723
MARS2
0.988
0.187
0.142
0.422
0.409
0.984
0.262
0.194
0.609
0.602
MARS3
0.986
0.203
0.153
0.472
0.466
0.984
0.227
0.184
0.461
0.453
MARS4
0.992
0.159
0.12
0.368
0.36
0.986
0.219
0.163
0.523
0.513
MARS5
0.992
0.155
0.117
0.348
0.344
0.99
0.323
0.266
0.687
0.676
MARS6
0.996
0.105
0.07
0.275
0.265
0.994
0.172
0.108
0.464
0.458
BAXGB1
0.978
0.261
0.189
0.614
0.594
0.978
0.261
0.204
0.58
0.576
BAXGB2
0.988
0.203
0.146
0.49
0.468
0.982
0.229
0.173
0.524
0.51
BAXGB3
0.99
0.176
0.13
0.412
0.402
0.986
0.19
0.142
0.416
0.409
BAXGB4
1
0.051
0.038
0.116
0.113
0.988
0.204
0.156
0.467
0.452
BAXGB5
0.992
0.205
0.145
0.502
0.487
0.988
0.222
0.158
0.534
0.532
BAXGB6
0.996
0.162
0.113
0.393
0.378
0.988
0.175
0.123
0.412
0.401
Note. Bold font indicates the best model.FigureS1 presents scatter diagrams between the estimated ETo and the P-M ETo (mm d−1) during the testing stage for the Yinchuan station. In Figure S1, the R2 value ranged at 0.871–0.974, 0.974–0.994, and 0.978–0.988 for GPR, MARS, and BAXGB among the input combinations, respectively. The slopes of the GPR models were less than 0.95 for input combinations 1, 3, 4, and 6, which led to a problem of underestimation. The slopes of MARS models were greater than 1 in all input combinations except the combination 3, suggesting a problem of overestimation. Overall, BAXGB and MARS models were superior to GPR models in the Yinchuan station.Comparisons between the estimated ETo and the P-M ETo during the testing stage for the three AI algorithms in Yinchuan are presented in FigureS2. P-M ETo had a peak in June (nearly 5 mm·d−1) but was generally small in winter (averaged at about 1 mm·d−1). For GPR models, the deviations of input combinations 1 and 5 were the most obvious, in which the combination 1 was seriously overestimated by 0.6 mm·d−1 or more from August to November, while the combination 5 was overvalued every month at a range of 0.2–0.6 mm·d−1. For MARS models, input combinations were performed similarly to the GPR model, but the extent of overestimation was less. In addition, MARS models in combinations 1 and 2 were significantly overestimated from May to July. For the BAXGB model, it had the smallest deviations (less than 0.4 mm·d−1) in each month among the three AI algorithms.


### 3.3. Model Performance in the Taiyuan Station

The performance of the three artificial intelligence algorithms in the training and testing stages for the Taiyuan station is compared in Table6. The statistical results in Taiyuan were similar to those in Datong because Datong is geographically closer to Taiyuan when compared with Yinchuan. The most significant difference between Taiyuan and Datong was that MARS models showed a difference in the importance rank of different meteorological parameters. However, the BAXGB model showed that solar radiation was the most important, followed by wind speed and relative humidity. Considering that each statistical index of the BAXGB model was better than the other two models, it is fair to suggest that solar radiation in this area was more important than other meteorological parameters.Table 6
Results of the statistical indicators for all models in the Taiyuan station during the training and testing stages.
Model
Training
Testing
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
GPR1
0.812
0.605
0.466
1.34
1.311
0.817
0.714
0.564
1.457
1.451
GPR2
0.972
0.248
0.19
0.561
0.55
0.96
0.418
0.311
0.949
0.942
GPR3
0.97
0.244
0.18
0.594
0.564
0.889
0.534
0.317
1.487
1.328
GPR4
0.972
0.243
0.176
0.582
0.564
0.98
0.403
0.31
0.835
0.827
GPR5
0.976
0.223
0.16
0.545
0.534
0.955
0.391
0.297
0.85
0.842
GPR6
0.978
0.22
0.152
0.559
0.546
0.976
0.397
0.301
0.858
0.851
MARS1
0.949
0.317
0.24
0.712
0.703
0.935
0.497
0.351
1.138
1.115
MARS2
0.97
0.241
0.187
0.542
0.534
0.976
0.377
0.275
0.865
0.853
MARS3
0.974
0.224
0.163
0.551
0.52
0.96
0.377
0.259
0.954
0.884
MARS4
0.984
0.178
0.132
0.422
0.406
0.986
0.364
0.284
0.811
0.803
MARS5
0.982
0.189
0.137
0.453
0.438
0.958
0.385
0.282
0.868
0.859
MARS6
0.988
0.151
0.101
0.39
0.361
0.978
0.33
0.22
0.817
0.809
BAXGB1
0.972
0.241
0.181
0.555
0.542
0.939
0.412
0.302
0.928
0.918
BAXGB2
0.994
0.115
0.087
0.266
0.256
0.97
0.314
0.24
0.712
0.699
BAXGB3
0.984
0.227
0.16
0.556
0.532
0.976
0.259
0.193
0.581
0.58
BAXGB4
0.996
0.097
0.071
0.232
0.221
0.986
0.297
0.235
0.643
0.632
BAXGB5
0.996
0.082
0.063
0.186
0.182
0.964
0.308
0.237
0.722
0.715
BAXGB6
1
0.036
0.026
0.083
0.08
0.988
0.253
0.182
0.594
0.578
Note. Bold font indicates the best model.FigureS3 presents scatter diagrams between the estimated ETo and the P-M ETo (mm d−1) during the testing stage for the Yinchuan station. In Figure S3, the R2 value ranged at 0.817–0.981, 0.935–0.986, and 0.940–0.987 for GPR, MARS, and BAXGB among the input combinations, respectively. The slopes of the three types of artificial intelligence models were all less than 0.93 for all combinations, suggesting that all combinations had been underestimated. Like results in other stations, the BAXGB and MARS models were superior to the GPR models in Taiyuan.FigureS4 presents comparison between the estimated ETo and the P-M ETo (mm d−1) during the testing stage for the three AI algorithms in the Taiyuan station. The monthly fluctuation of P-M ETo in Taiyuan was similar to that of Datong, with the peak appeared in May (nearly 5 mm d−1) as well. For the GPR model, the combination 1 was underestimated in spring but overestimated in autumn, the maximum deviation of which was over 1 mm d−1. For MARS models, the input combination 3 had a severe underestimation in April, while other combinations had overestimations at spring and summer. Although there was an overall trend of underestimation for the BAXGB models, combinations 1, 3, and 4 were underestimated by more than 0.4 mm d−1 in some months (in spring usually).


### 3.4. Comparison between Different Algorithms

Overall, the three AI algorithms presented different optimal combinations at different stations. For example, GPR models in combination 4 yielded better performance than in combination 6 across the three stations. On the contrary, BAXGB models identified that the combination 6 was the optimal input combination. This indicated that the BAXGB algorithm might have more advantages in identifying the contributions of multiple parameters compared with the other two algorithms, and therefore, its accuracy was the best.Given the fact that BAXGB overall performed better than GPR and MARS across the three stations and the six input combinations, statistical results of BAXGB can be used for analyzing the importance of each major meteorological parameter for estimating ETo. Based on the results of BAXGB models in this study, solar radiation is more important than the others in Yinchuan and Taiyuan stations. This is similar to other studies [24, 36], where they show that solar radiation is more important than wind speed and relative humidity in a semiarid environment of Qilian Mountain of Northwest China and Southeast China. However, relative humidity and wind speed seem more important than solar radiation for the Datong station with BAXGB models in this study, which may be related to the difference in the distribution of solar radiation data between the training and testing stages. Compared with the input combinations 4 and 5, statistical indicators of the input combination 3 are smaller during the training stage. These differences are much larger in the testing stage than the training stage, especially for values of RMSEmax10% and MAEmax10%.Tao et al. [28] developed a coupling model based on the firefly algorithm and ANFIS in Burkina Faso, a semiarid region of Africa, in which the RMSE in the testing period ranged at 0.24–1.10 mm×d−1 with different input combinations. Shamshirband [37] found that the cuckoo search algorithm optimized ANN and ANFIS had an RMSE value of 0.330 and 0.265 mm d−1 in the testing stage at Serbia, respectively. Another case study in Serbia showed that the genetic programming models coupling with GA generated RMSE at 0.649–1.188 and 0.657–1.193 mm d−1 in the testing stage and the validation, respectively [38]. In addition, a hybrid model (GA-SVM) conducted in the Pailugou Watershed, Northwest China, showed that the RMSE values ranged at 0.138–0.424 mm d−1 [21, 36]. Overall, RMSE values obtained in this study are similar to those achieved by other researchers using coupled artificial intelligence models in semiarid regions. Furthermore, in the case of limited observational data in this region, BAXGB models with only temperature and extraterrestrial radiation as inputs would offer an acceptable result and can be used as an alternative method for ETo estimation. On this basis, adding one or two meteorological factors would also be beneficial.Over the past decade, many highly efficient evolutionary algorithms have been proposed. The XGBoost is also a relatively new technology, and the study of combining evolutionary algorithms with XGBoost is very limit to date. This study confirms that the bat algorithm is an effective method for optimizing XGBoost. However, whether other evolutionary algorithms would be more efficient than BA for XGBoost optimization needs to be further studied.


## 3.1. Model Performance in the Datong Station

Performance of the three artificial intelligence models used in this study (i.e., GPR, MARS, and BAXGB) during the training and testing stages for the Datong station is listed in Table4. Both air temperature and extraterrestrial radiation were used for the first two input combinations. The difference was that the input combination 1 used the monthly average temperature, while the input combination 2 used the monthly maximum and minimum temperatures. Under these two input combinations, the accuracies of the three types of models varied. For GPR and BAXGB models, input combination 2 had higher accuracy than input combination 1. However, MARS models yielded an opposite result. Input combinations 3, 4, and 5 represented the effects of solar radiation, wind speed, and relative humidity or temperature data, respectively. Considering that the extraterrestrial radiation is only related to the latitude of the Earth and does not belong to the meteorological parameters to be observed, this information was added to the combinations 4 and 5. The importance of meteorological factors given by different models was different. MARS and BAXGB models showed that relative humidity was the most important, followed by wind speed and solar radiation. However, GPR model results showed that wind speed was the most critical, followed by relative humidity and solar radiation. Although there was difference to some degree in the importance between different meteorological parameters, it was still clear that relative humidity and wind speed were more important than solar radiation for the Datong station. Input combination 6 was a combination of temperature, wind speed, and relative humidity. For the three artificial intelligence models, results of input combination 6 were better than other combinations that with relative humidity or wind speed alone. Comparisons between the three artificial intelligence models under the same input combination indicated that BAXGB (RMSE = 0.116–0.438 mm·d−1 and MAE = 0.087–0.311 mm·d−1) always performed better than MARS (RMSE = 0.153–0.512 mm·d−1 and MAE = 0.221– 0.559 mm·d−1) and GPR (RMSE = 0.297–0.695 mm·d−1 and MAE = 0.120–0.379 mm·d−1) under each input combination. In addition, the MARS model generally performed better than the GPR model, except in combinations 2 and 3. Across algorithms and input combinations, values of RMSEmax10% and MAEmax10% during the testing stage were 2–2.5 times of their corresponding values of RMSE and MAE, suggesting that model performance in extreme cases would be worse than normal conditions. Although the overall performance of the BAXGB model was good, the maximum error point of 10% with input combination 3 (i.e., BAXGB3) exceeded 1.0 mm·d−1. This may be due to the difference in distribution of solar radiation data between the training and the testing stage.Table 4
Results of the statistical indicators for all models in the Datong station during the training and testing stages.
Model
Training
Testing
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
GPR1
0.841
0.66
0.519
1.392
1.373
0.832
0.695
0.559
1.333
1.328
GPR2
0.978
0.258
0.196
0.586
0.563
0.972
0.372
0.302
0.744
0.738
GPR3
0.98
0.242
0.181
0.565
0.547
0.956
0.444
0.35
0.898
0.891
GPR4
0.974
0.279
0.202
0.673
0.645
0.978
0.297
0.221
0.668
0.664
GPR5
0.988
0.197
0.138
0.483
0.46
0.951
0.398
0.251
1.089
1.019
GPR6
0.986
0.209
0.139
0.548
0.517
0.958
0.37
0.25
0.919
0.877
MARS1
0.964
0.315
0.236
0.723
0.71
0.97
0.292
0.227
0.649
0.631
MARS2
0.986
0.199
0.15
0.465
0.453
0.962
0.421
0.337
0.821
0.812
MARS3
0.984
0.209
0.157
0.481
0.468
0.953
0.512
0.379
1.158
1.148
MARS4
0.99
0.162
0.121
0.371
0.363
0.98
0.271
0.216
0.559
0.552
MARS5
0.992
0.154
0.11
0.377
0.362
0.992
0.153
0.12
0.321
0.32
MARS6
0.994
0.122
0.087
0.3
0.292
0.992
0.189
0.121
0.509
0.483
BAXGB1
0.978
0.238
0.179
0.541
0.531
0.972
0.291
0.23
0.592
0.587
BAXGB2
0.988
0.209
0.148
0.518
0.505
0.974
0.284
0.227
0.576
0.568
BAXGB3
0.998
0.076
0.057
0.175
0.168
0.955
0.438
0.311
1.067
1.057
BAXGB4
0.996
0.092
0.069
0.214
0.207
0.988
0.191
0.153
0.403
0.398
BAXGB5
0.996
0.108
0.076
0.266
0.255
0.992
0.146
0.112
0.337
0.328
BAXGB6
1
0.03
0.022
0.068
0.066
0.996
0.116
0.087
0.272
0.258
Note. Bold font indicates the best model.Figure4 presents scatter diagrams between the estimated ETo by the three artificial intelligence algorithms and the P-M ETo (mm·d−1) during the testing stage for the Datong station. In Figure 4, the R2 value ranged at 0.832–0.978, 0.953–0.992, and 0.955–0.996 for GPR, MARS, and BAXGB for all input combinations, respectively. The slopes of the GPR models were all less than 0.95 except with the input combination 3, which led to a problem of underestimation. These results clearly indicated that, for the Datong station, BAXGB was the most superior, while GPR was the worst among the AI algorithms, in terms of the performance in ETo estimation.Figure 4
Scatter plots between the estimated ETo by the three algorithms and the FAO-56 P-M ETo in the testing stage for the Datong station.
(a)
(b)
(c)
(d)
(e)
(f)Figure5 presents comparisons between the estimated ETo by AI algorithms and the P-M ETo (mm·d−1) during the testing stage for the Datong station. P-M ETo peaked in May with a value over 5 mm·d−1 while had the lowest value in December (0.6 mm·d−1). For the GPR model, the simulation effect of input combination 1 did not catch the trend well, resulting in a deviation more than 1 mm·d−1, presented as underestimation in May and overestimation in September, respectively. For the MARS model, input combinations 2 and 3 had a problem of overall overestimation, especially for months from August to October where the overestimation was more than 0.4 mm·d−1. For the BAXGB model, all input combination simulated an acceptable trend, except for the input combination 3 where an overestimation (more than 0.4 mm·d−1) was detected from July to October. Overall, among the AI methods in this study, the BAXGB model had the smallest deviation in all the 12 months, while GPR models performed worst.Figure 5
Relative errors by the different algorithms at each month for the Datong station.
(a)
(b)
(c)
(d)
(e)
(f)To further evaluate the performance of different AI models, Taylor diagrams were used to analyze the performance of models under different combinations, and the results of the Datong station can be seen in Figure6. It is not hard to see that combination 4 performs best in the GPR models, while combination 5 performs best in the MARS models. However, the combination 6 is superior to others in the BAXGB models. In addition, the point of the BAXGB model shows much closer than the other two models. This is consistent with the previous statistical indicators and scatter plot results.Figure 6
The Taylor diagram of (a) GPR, (b) MARS, and (c) BAXGB at the Datong station.
(a)
(b)
(c)


## 3.2. Model Performance in the Yinchuan Station

The performance of the three artificial intelligence models in the training and testing stages for the Yinchuan station is compared in Table5. In the first two input combinations, combination 2 had higher accuracy than combination 1 for all the three algorithms used in Yinchuan, which was different from Datong. However, the rank of importance of meteorological factors in different models showed similar trends with the Datong station. In addition, the combination including both wind speed and relative humidity (i.e., input combination 6) was still better than either factor alone for MARS and BAXGB models. It can be seen from the comparison of the three artificial intelligence models under the same input combination that BAXGB (RMSE = 0.175–0.261  mm·d−1 and MAE = 0.123–0.204 mm·d−1) performed better than MARS (RMSE = 0.172–0.340 mm·d−1 and MAE = 0.108– 266 mm·d−1) and GPR (RMSE = 0.289–0.663 mm·d−1 and MAE = 0.195–0.511 mm·d−1) under each combination. In general, the BAXGB model was better than the MARS model, and the GPR model was the worst, except for some individual combinations. Similar to the Datong station, RMSEmax10% and MAEmax10% during the testing stage were more than twice as high as the corresponding RMSE and MAE in the Yinchuan station. However, compared with Datong, values of RMSEmax10% and MAEmax10% were generally smaller in Yinchuan.Table 5
Results of the statistical indicators for all models in the Yinchuan station during the training and testing stages.
Model
Training
Testing
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
GPR1
0.856
0.653
0.514
1.371
1.351
0.87
0.663
0.511
1.424
1.401
GPR2
0.982
0.248
0.191
0.568
0.548
0.974
0.334
0.259
0.73
0.719
GPR3
0.984
0.23
0.176
0.53
0.518
0.974
0.289
0.197
0.751
0.696
GPR4
0.98
0.257
0.196
0.588
0.574
0.966
0.343
0.262
0.752
0.736
GPR5
0.984
0.23
0.169
0.536
0.519
0.972
0.461
0.378
0.941
0.907
GPR6
0.986
0.208
0.147
0.507
0.491
0.97
0.384
0.285
0.888
0.878
MARS1
0.968
0.3
0.229
0.677
0.666
0.974
0.34
0.262
0.728
0.723
MARS2
0.988
0.187
0.142
0.422
0.409
0.984
0.262
0.194
0.609
0.602
MARS3
0.986
0.203
0.153
0.472
0.466
0.984
0.227
0.184
0.461
0.453
MARS4
0.992
0.159
0.12
0.368
0.36
0.986
0.219
0.163
0.523
0.513
MARS5
0.992
0.155
0.117
0.348
0.344
0.99
0.323
0.266
0.687
0.676
MARS6
0.996
0.105
0.07
0.275
0.265
0.994
0.172
0.108
0.464
0.458
BAXGB1
0.978
0.261
0.189
0.614
0.594
0.978
0.261
0.204
0.58
0.576
BAXGB2
0.988
0.203
0.146
0.49
0.468
0.982
0.229
0.173
0.524
0.51
BAXGB3
0.99
0.176
0.13
0.412
0.402
0.986
0.19
0.142
0.416
0.409
BAXGB4
1
0.051
0.038
0.116
0.113
0.988
0.204
0.156
0.467
0.452
BAXGB5
0.992
0.205
0.145
0.502
0.487
0.988
0.222
0.158
0.534
0.532
BAXGB6
0.996
0.162
0.113
0.393
0.378
0.988
0.175
0.123
0.412
0.401
Note. Bold font indicates the best model.FigureS1 presents scatter diagrams between the estimated ETo and the P-M ETo (mm d−1) during the testing stage for the Yinchuan station. In Figure S1, the R2 value ranged at 0.871–0.974, 0.974–0.994, and 0.978–0.988 for GPR, MARS, and BAXGB among the input combinations, respectively. The slopes of the GPR models were less than 0.95 for input combinations 1, 3, 4, and 6, which led to a problem of underestimation. The slopes of MARS models were greater than 1 in all input combinations except the combination 3, suggesting a problem of overestimation. Overall, BAXGB and MARS models were superior to GPR models in the Yinchuan station.Comparisons between the estimated ETo and the P-M ETo during the testing stage for the three AI algorithms in Yinchuan are presented in FigureS2. P-M ETo had a peak in June (nearly 5 mm·d−1) but was generally small in winter (averaged at about 1 mm·d−1). For GPR models, the deviations of input combinations 1 and 5 were the most obvious, in which the combination 1 was seriously overestimated by 0.6 mm·d−1 or more from August to November, while the combination 5 was overvalued every month at a range of 0.2–0.6 mm·d−1. For MARS models, input combinations were performed similarly to the GPR model, but the extent of overestimation was less. In addition, MARS models in combinations 1 and 2 were significantly overestimated from May to July. For the BAXGB model, it had the smallest deviations (less than 0.4 mm·d−1) in each month among the three AI algorithms.


## 3.3. Model Performance in the Taiyuan Station

The performance of the three artificial intelligence algorithms in the training and testing stages for the Taiyuan station is compared in Table6. The statistical results in Taiyuan were similar to those in Datong because Datong is geographically closer to Taiyuan when compared with Yinchuan. The most significant difference between Taiyuan and Datong was that MARS models showed a difference in the importance rank of different meteorological parameters. However, the BAXGB model showed that solar radiation was the most important, followed by wind speed and relative humidity. Considering that each statistical index of the BAXGB model was better than the other two models, it is fair to suggest that solar radiation in this area was more important than other meteorological parameters.Table 6
Results of the statistical indicators for all models in the Taiyuan station during the training and testing stages.
Model
Training
Testing
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
R
2
RMSE
MAE
RMSEmax10%
MAEmax10%
GPR1
0.812
0.605
0.466
1.34
1.311
0.817
0.714
0.564
1.457
1.451
GPR2
0.972
0.248
0.19
0.561
0.55
0.96
0.418
0.311
0.949
0.942
GPR3
0.97
0.244
0.18
0.594
0.564
0.889
0.534
0.317
1.487
1.328
GPR4
0.972
0.243
0.176
0.582
0.564
0.98
0.403
0.31
0.835
0.827
GPR5
0.976
0.223
0.16
0.545
0.534
0.955
0.391
0.297
0.85
0.842
GPR6
0.978
0.22
0.152
0.559
0.546
0.976
0.397
0.301
0.858
0.851
MARS1
0.949
0.317
0.24
0.712
0.703
0.935
0.497
0.351
1.138
1.115
MARS2
0.97
0.241
0.187
0.542
0.534
0.976
0.377
0.275
0.865
0.853
MARS3
0.974
0.224
0.163
0.551
0.52
0.96
0.377
0.259
0.954
0.884
MARS4
0.984
0.178
0.132
0.422
0.406
0.986
0.364
0.284
0.811
0.803
MARS5
0.982
0.189
0.137
0.453
0.438
0.958
0.385
0.282
0.868
0.859
MARS6
0.988
0.151
0.101
0.39
0.361
0.978
0.33
0.22
0.817
0.809
BAXGB1
0.972
0.241
0.181
0.555
0.542
0.939
0.412
0.302
0.928
0.918
BAXGB2
0.994
0.115
0.087
0.266
0.256
0.97
0.314
0.24
0.712
0.699
BAXGB3
0.984
0.227
0.16
0.556
0.532
0.976
0.259
0.193
0.581
0.58
BAXGB4
0.996
0.097
0.071
0.232
0.221
0.986
0.297
0.235
0.643
0.632
BAXGB5
0.996
0.082
0.063
0.186
0.182
0.964
0.308
0.237
0.722
0.715
BAXGB6
1
0.036
0.026
0.083
0.08
0.988
0.253
0.182
0.594
0.578
Note. Bold font indicates the best model.FigureS3 presents scatter diagrams between the estimated ETo and the P-M ETo (mm d−1) during the testing stage for the Yinchuan station. In Figure S3, the R2 value ranged at 0.817–0.981, 0.935–0.986, and 0.940–0.987 for GPR, MARS, and BAXGB among the input combinations, respectively. The slopes of the three types of artificial intelligence models were all less than 0.93 for all combinations, suggesting that all combinations had been underestimated. Like results in other stations, the BAXGB and MARS models were superior to the GPR models in Taiyuan.FigureS4 presents comparison between the estimated ETo and the P-M ETo (mm d−1) during the testing stage for the three AI algorithms in the Taiyuan station. The monthly fluctuation of P-M ETo in Taiyuan was similar to that of Datong, with the peak appeared in May (nearly 5 mm d−1) as well. For the GPR model, the combination 1 was underestimated in spring but overestimated in autumn, the maximum deviation of which was over 1 mm d−1. For MARS models, the input combination 3 had a severe underestimation in April, while other combinations had overestimations at spring and summer. Although there was an overall trend of underestimation for the BAXGB models, combinations 1, 3, and 4 were underestimated by more than 0.4 mm d−1 in some months (in spring usually).


## 3.4. Comparison between Different Algorithms

Overall, the three AI algorithms presented different optimal combinations at different stations. For example, GPR models in combination 4 yielded better performance than in combination 6 across the three stations. On the contrary, BAXGB models identified that the combination 6 was the optimal input combination. This indicated that the BAXGB algorithm might have more advantages in identifying the contributions of multiple parameters compared with the other two algorithms, and therefore, its accuracy was the best.Given the fact that BAXGB overall performed better than GPR and MARS across the three stations and the six input combinations, statistical results of BAXGB can be used for analyzing the importance of each major meteorological parameter for estimating ETo. Based on the results of BAXGB models in this study, solar radiation is more important than the others in Yinchuan and Taiyuan stations. This is similar to other studies [24, 36], where they show that solar radiation is more important than wind speed and relative humidity in a semiarid environment of Qilian Mountain of Northwest China and Southeast China. However, relative humidity and wind speed seem more important than solar radiation for the Datong station with BAXGB models in this study, which may be related to the difference in the distribution of solar radiation data between the training and testing stages. Compared with the input combinations 4 and 5, statistical indicators of the input combination 3 are smaller during the training stage. These differences are much larger in the testing stage than the training stage, especially for values of RMSEmax10% and MAEmax10%.Tao et al. [28] developed a coupling model based on the firefly algorithm and ANFIS in Burkina Faso, a semiarid region of Africa, in which the RMSE in the testing period ranged at 0.24–1.10 mm×d−1 with different input combinations. Shamshirband [37] found that the cuckoo search algorithm optimized ANN and ANFIS had an RMSE value of 0.330 and 0.265 mm d−1 in the testing stage at Serbia, respectively. Another case study in Serbia showed that the genetic programming models coupling with GA generated RMSE at 0.649–1.188 and 0.657–1.193 mm d−1 in the testing stage and the validation, respectively [38]. In addition, a hybrid model (GA-SVM) conducted in the Pailugou Watershed, Northwest China, showed that the RMSE values ranged at 0.138–0.424 mm d−1 [21, 36]. Overall, RMSE values obtained in this study are similar to those achieved by other researchers using coupled artificial intelligence models in semiarid regions. Furthermore, in the case of limited observational data in this region, BAXGB models with only temperature and extraterrestrial radiation as inputs would offer an acceptable result and can be used as an alternative method for ETo estimation. On this basis, adding one or two meteorological factors would also be beneficial.Over the past decade, many highly efficient evolutionary algorithms have been proposed. The XGBoost is also a relatively new technology, and the study of combining evolutionary algorithms with XGBoost is very limit to date. This study confirms that the bat algorithm is an effective method for optimizing XGBoost. However, whether other evolutionary algorithms would be more efficient than BA for XGBoost optimization needs to be further studied.


## 4. Conclusions

Accurate and robust reference evapotranspiration (ETo) estimating models are of importance for agriculture and irrigation scheduling. In this study, a new model (BAXGB) that couples a bat algorithm with XGBoost was proposed for the estimation of monthly ETo. The developed BAXGB algorithm was studied by comparing with the FAO 56 Penman–Monteith model results in arid and semiarid environments of China. Data from three stations (Datong, Yinchuan, and Taiyuan) were used to build the BAXGB model as well as the multivariate adaptive regression splines (MARS) and the gaussian process regression (GPR) models. Six models for each algorithm were developed using different input combinations including mean air temperature, maximum air temperature, minimum air temperature, wind speed, relative humidity, and solar radiation or extraterrestrial radiation. The results showed that the BAXGB models exhibited a suitable accuracy for estimating ETo. Except for the input combination 6 (Tmax, Tmin, U, and RH), BAXGB models were superior to MARS models. The GPR model was generally worse than the MARS model. For the Datong station, the influence of solar radiation on model accuracy was less than that of wind speed and relative humidity. However, for the other two stations, there was an opposite trend. Findings of this research will be helpful for the agricultural and irrigation scheduling and may offer reference to other regions where accurate models for improving the local water management are needed.



---

*Source: 9575782-2019-10-17.xml*