How Can We Realize Sustainable Development Goals in Rocky Desertiﬁed Regions by Enhancing Crop Yield with Reduction of Environmental Risks?

: To meet the sustainable development goals in rocky desertiﬁed regions like Guizhou Province in China, we should maximize the crop yield with minimal environmental costs. In this study, we ﬁrst calculated the yield gap for 6 main crop species in Guizhou Province and evaluated the quantitative relationships between crop yield and inﬂuencing variables utilizing ensembled artiﬁcial neural networks. We also tested the inﬂuence of adjusting the quantity of local fertilization and irrigation on crop production in Guizhou Province. Results showed that the total yield of the selected crops had, on average, reached over 72.5% of the theoretical maximum yield. Increasing irrigation tended to be more consistently effective at increasing crop yield than additional fertilization. Conversely, appropriate reduction of fertilization may even beneﬁt crop yield in some regions, simultaneously resulting in signiﬁcantly higher fertilization efﬁciency with lower residuals in the environment. The total positive impact of continuous intensiﬁcation of irrigation and fertilization on most crop species was limited. Therefore, local stakeholders are advised to consider other agricultural management measures to improve crop yield in this region.


Introduction
In 2015, the UN Sustainable Development Goals (SDGs) provided a shared blueprint for peace and prosperity for people and the planet [1]. The second, sixth and fifteenth themes of this blueprint set specific goals for food security, clean water supply and halting land degradation [2]; requiring balancing of agricultural development and the environment [3]. Over the past decades, the intensification of irrigation and fertilization has greatly improved food production globally, but at the cost of environmental pollutions like soil salination and eutrophicated water resources [4,5]. The difficulty of achieving these SDGs is spatially unbalanced. This is especially evident for some underdeveloped areas due to their environmental and economical limitations [6]. For example, in karst regions, rocky desertification, a process of land degradation involving serious soil erosion, results in extensive exposure of basement rocks, a drastic decrease in soil productivity, and the appearance of a desert-like landscape [7][8][9]. The poor fertility of the edaphic environment, along withwater and soil losses, induces low crop productivity [10,11]. Nitrogen and other nutrition loss from the karst system also lead to lower nutrient use efficiency of vegetation and contaminates water quality, further impacting the long-term soil security and sustainability in karst regions [12].
Guizhou province, located in southwest China, contains large karst regions with their highly sensitive and vulnerable ecosystem and is one of the world's largest exposed relationship of different factors and crop yield with nonlinear regressions [45]. Based on ANN training, higher accuracy in the simulation of crop yield could be achieved [46,47].
In this paper, we identified the current yield gap in Guizhou Province and investigate strategies to close this gap. We first calculated the yield gap for 6 dominant crop species grown in Guizhou Province. We subsequently employed ensembled back propagation (BP) networks to simulate the crop yield and fertilizer balance (defined as the difference of fertilization investment and demand for crop growth) in the study region by inputting different influencing variables inside the networks. Lastly, we forecast the variation in crop yield and fertilization efficiency under different agricultural management scenarios. The results of this work are useful as a DST for policymakers and other stakeholders to optimize food demand and environmental impact in agriculture development.

Study Region
Guizhou Province is located in southwest China, at the heart of the East Asia Karst, one of the three largest areas of almost unbroken karst globally [48] (Figure 1). There are 9 prefectures spatially distributed inside Guizhou. However, due to historical reasons and its geological environment, which causes difficulties with transportation and communication with outsiders, Guizhou is still one of the less developed provinces in China [49][50][51]. In this study, we chose 6 crops grown extensively in the 9 prefectures of Guizhou Province for the simulation of yield, which are maize, potato, rice, soybean and wheat (food crops), along with groundnut (cash crop).
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 24 learning methods, such as artificial neural networks (ANN), provides a new solution to identify the relationship of different factors and crop yield with nonlinear regressions [45]. Based on ANN training, higher accuracy in the simulation of crop yield could be achieved [46,47].
In this paper, we identified the current yield gap in Guizhou Province and investigate strategies to close this gap. We first calculated the yield gap for 6 dominant crop species grown in Guizhou Province. We subsequently employed ensembled back propagation (BP) networks to simulate the crop yield and fertilizer balance (defined as the difference of fertilization investment and demand for crop growth) in the study region by inputting different influencing variables inside the networks. Lastly, we forecast the variation in crop yield and fertilization efficiency under different agricultural management scenarios. The results of this work are useful as a DST for policymakers and other stakeholders to optimize food demand and environmental impact in agriculture development.

Study Region
Guizhou Province is located in southwest China, at the heart of the East Asia Karst, one of the three largest areas of almost unbroken karst globally [48] (Figure 1). There are 9 prefectures spatially distributed inside Guizhou. However, due to historical reasons and its geological environment, which causes difficulties with transportation and communication with outsiders, Guizhou is still one of the less developed provinces in China [49][50][51]. In this study, we chose 6 crops grown extensively in the 9 prefectures of Guizhou Province for the simulation of yield, which are maize, potato, rice, soybean and wheat (food crops), along with groundnut (cash crop).

Dataset
To simulate crop yield (per unit area), 14 variables were selected as input factors (Table 1, excluding crop yield, crop area, N and P balance), of which annual total rainfall, soil bulk density, irrigation, nitrogen (N) fertilizer rate, phosphorus (P) fertilizer rate and slope were chosen to simulate the Nor P balance. Specifically, we extracted all the data from (circa) year 2000 and unified their spatial resolution into 5′ by aggregation (overlaying meteorological data with a spatial resolution of 5° on 5′ data and extracted meteorological

Dataset
To simulate crop yield (per unit area), 14 variables were selected as input factors (Table 1, excluding crop yield, crop area, N and P balance), of which annual total rainfall, soil bulk density, irrigation, nitrogen (N) fertilizer rate, phosphorus (P) fertilizer rate and slope were chosen to simulate the Nor P balance. Specifically, we extracted all the data from (circa) year 2000 and unified their spatial resolution into 5 by aggregation (overlaying meteorological data with a spatial resolution of 5 • on 5 data and extracted meteorological data values for each pixel of 5 ). The dataset of crop yield was created by merging a detailed library of global census data (from the national, state and county level for over 20,000 political units) with three different, detailed satellite datasets of global landcover conditions, and this is by far the most complete dataset of agriculture [52].

Calculation of Yield Gap
To calculate the yield gap for the selected species, we referred to the approaches raised in the papers of Foley's group [66]. First, we conducted an unsupervised spatial classification method (self-organizing feature mapping, or SOFM) to define different zones of "natural environment" (natural zones hereinafter) based on temperature, rainfall, shortwave radiation and slope, which human behaviors cannot change over a short period. This idea allowed us to control the same basic environmental conditions for crop growth. Inside each natural zone, we assumed that there is the same theoretical value of maximum crop yield, which is calculated based on the true yield distribution. Lastly, by comparing the true (observed) value with the theoretically maximum value, we finally calculated the yield gap in each prefecture.

SOFM
The SOFM network ( Figure 2) consists of a fully interconnected array of neurons with a topology input layer and the competition layer [67]. The inputs are connected to each node on the network grid, and each grid node is only connected to adjacent nodes. They use a neighborhood function to preserve the input space's topological properties and realize unsupervised classification by different data resources. The training process is completed to adjust the weight between each node in the output layer until it meets the fixed terminal condition. The whole process of this work is to reflect the input variables into lower dimensional space to realize the classification. The specific steps are as follows [68]: (1) Initializing each weight with a group of random numbers at the beginning; (2) Selecting the best matching unit. This step is also regarded as a competing process.
The weight vector, which has the minimum Euclidean distance (calculated as below), with sample x (randomly selected), is the winning unit: where c represents the winning unit, and i is the sequence number for x and the weight vector; (3) Updating the weight until it meets the terminal conditions: where p i and p c indicate the positions of the output units of i and c, respectively, while σ is the width of the neighborhood function.
where c represents the winning unit, and i is the sequence number for x and the weight vector; (3) Updating the weight until it meets the terminal conditions:    (2) where pi and pc indicate the positions of the output units of i and c, respectively, while σ is the width of the neighborhood function.

Calculation of Yield Gap
The yield gap was calculated using three steps: (1) For each natural zone, the grid cells were sorted by their yield value, ranked from lowest to highest; (2) The grid cells' respective harvested areas were accumulated into a histogram, then the percentile of 5th and 95th were extracted ( Figure 3). The corresponding yield data from these two points were defined as theoretical minimum and maximum yield, indicating the potential yield for crops growing under purely natural conditions and under optimal agricultural management inside each zone. Extraction from the histogram of harvested area, instead of yield, can exclude the grid cells with too large or small cultivated areas that could skew the distribution of yield histogram; (3) For each grid cell, the crop yield gap was the difference between actual yield per area (AYPA) and potential yield per area (PYPA-referring to the theoretical yield under optimal growth conditions for each crop per area). Additionally, the yield gap of total production (AYPA/PYPA multiplied by corresponding cultivated area) was also calculated (shortened as AY and PA). To provide a better reference for policymakers and other stakeholders, we integrated each cell's gap into the yield gap for each prefecture

Calculation of Yield Gap
The yield gap was calculated using three steps: (1) For each natural zone, the grid cells were sorted by their yield value, ranked from lowest to highest; (2) The grid cells' respective harvested areas were accumulated into a histogram, then the percentile of 5th and 95th were extracted ( Figure 3). The corresponding yield data from these two points were defined as theoretical minimum and maximum yield, indicating the potential yield for crops growing under purely natural conditions and under optimal agricultural management inside each zone. Extraction from the histogram of harvested area, instead of yield, can exclude the grid cells with too large or small cultivated areas that could skew the distribution of yield histogram; (3) For each grid cell, the crop yield gap was the difference between actual yield per area (AYPA) and potential yield per area (PYPA-referring to the theoretical yield under optimal growth conditions for each crop per area). Additionally, the yield gap of total production (AYPA/PYPA multiplied by corresponding cultivated area) was also calculated (shortened as AY and PA). To provide a better reference for policymakers and other stakeholders, we integrated each cell's gap into the yield gap for each prefecture in Guizhou Province. in Guizhou Province.
The theoretical maximum and minimum value also acted as the limiting threshold for further simulation of crop yield under different management scenarios.

Single BP Network
BP network is widely applied in classification, simulation and forecast across different research disciplines, and in our previous work it offered the best balance between accuracy in predicting crop yield in Guizhou Province and time cost [69,70]. This supervised learning is conducted using the mean square error and gradient descent algorithm to achieve the correction of the network connection weights, to minimize the difference between the mean square error of the actual output to achieve prediction accuracy. The common BP network structure consists of an input layer, (single or multiple) hidden layer (s) and an output layer. The BP process integrates the value of each node (or neuron) by weights, bias and activation function, then transfers this into the output layer, while the error signal is transmitted back along the original connection path, through the network, to modify the weights of neurons until the desired target is reached [71,72]. The training process can be characterized by numbers of iteration, and each run of them is called an epoch. Before commencing a BP network, several hyper-parameters, which may indirectly influence the accuracy and training efficiency of the model, need to be set, such as the numbers of hidden layers and nodes, as well as learning rate, etc. In this study, we imported one hidden layer with 10 nodes and used a tan-sigmoid (Equation (3)) function and linear function to transfer values between the different layers. Meanwhile, after testing different groups of hyper-parameters, the learning rate of the network and maximum epoch number was ultimately assigned as 0.01 and 5000, respectively. The theoretical maximum and minimum value also acted as the limiting threshold for further simulation of crop yield under different management scenarios.

Ensembled BP (Back Propagation) Artificial Neural Networks 2.4.1. Single BP Network
BP network is widely applied in classification, simulation and forecast across different research disciplines, and in our previous work it offered the best balance between accuracy in predicting crop yield in Guizhou Province and time cost [69,70]. This supervised learning is conducted using the mean square error and gradient descent algorithm to achieve the correction of the network connection weights, to minimize the difference between the mean square error of the actual output to achieve prediction accuracy. The common BP network structure consists of an input layer, (single or multiple) hidden layer (s) and an output layer. The BP process integrates the value of each node (or neuron) by weights, bias and activation function, then transfers this into the output layer, while the error signal is transmitted back along the original connection path, through the network, to modify the weights of neurons until the desired target is reached [71,72]. The training process can be characterized by numbers of iteration, and each run of them is called an epoch. Before commencing a BP network, several hyper-parameters, which may indirectly influence the accuracy and training efficiency of the model, need to be set, such as the numbers of hidden layers and nodes, as well as learning rate, etc. In this study, we imported one hidden layer with 10 nodes and used a tan-sigmoid (Equation (3)) function and linear function to transfer values between the different layers. Meanwhile, after testing different groups of hyper-parameters, the learning rate of the network and maximum epoch number was ultimately assigned as 0.01 and 5000, respectively.

Evaluating Accuracy of Simulation
The dataset was first randomly grouped into training data (75% of all data) and validation data (the remaining 25%). Three indices of validation data were selected to evaluate the accuracy of training performance of the BP network, including the correlation coefficient (between the true (observed) value and forecasted (simulated) value of the validation group), root-mean-square error (RMSE) and relative mean error (RME), the latter Remote Sens. 2021, 13, 1614 7 of 23 two, of which indicate absolute and relative deviation of the simulation results, respectively (defined in Equations (4) and (5)).
where T and F represent true value and forecasted value, and n indicates the total number of validation samples.

Ensemble Approach
The ensemble approach for ANN is to integrate a finite number of ANNs that are trained for the same purpose whose predictions are combined to generate a unique output ( Figure 4). The relationship between a group of ensembled ANNs and a single ANN is similar to a random forest and each tree inside. These ensembles have several advantages over a single ANN as they have the potential for improved generalization and increased stability, although at the cost of a longer training time, enabling better convergence and accuracy [73]. The commonly used idea of an ensemble approach to gain a final result from a group of ANNs can be categorized into "bagging" and "boosting". They both build a robust classification or regression network based on a single ANN but with different specific algorithms [74]. Herein, for each process of simulation, we utilized the idea of "bagging". Namely, we created several networks (the number of which was set following performance testing of the simulation-detailed in Section 2.4.4) and subsequently calculated the average value for each as the final result of the simulation, and this is the commonly used nonlinear regression model ensemble algorithm.

Evaluating Accuracy of Simulation
The dataset was first randomly grouped into training data (75% of all data) and validation data (the remaining 25%). Three indices of validation data were selected to evaluate the accuracy of training performance of the BP network, including the correlation coefficient (between the true (observed) value and forecasted (simulated) value of the validation group), root-mean-square error (RMSE) and relative mean error (RME), the latter two, of which indicate absolute and relative deviation of the simulation results, respectively (defined in Equations (4) and (5)).
where T and F represent true value and forecasted value, and n indicates the total number of validation samples.

Ensemble Approach
The ensemble approach for ANN is to integrate a finite number of ANNs that are trained for the same purpose whose predictions are combined to generate a unique output ( Figure 4). The relationship between a group of ensembled ANNs and a single ANN is similar to a random forest and each tree inside. These ensembles have several advantages over a single ANN as they have the potential for improved generalization and increased stability, although at the cost of a longer training time, enabling better convergence and accuracy [73]. The commonly used idea of an ensemble approach to gain a final result from a group of ANNs can be categorized into "bagging" and "boosting". They both build a robust classification or regression network based on a single ANN but with different specific algorithms [74]. Herein, for each process of simulation, we utilized the idea of "bagging". Namely, we created several networks (the number of which was set following performance testing of the simulation-detailed in Section 2.4.4) and subsequently calculated the average value for each as the final result of the simulation, and this is the commonly used nonlinear regression model ensemble algorithm.

Scenario Settings and Convergence Test of Ensembled ANNs
To evaluate different future scenarios of agricultural management, we changed the input of N&P fertilizer (by ±10%) and irrigation rates (only by +10%, as the study region was proven to be lacking water resources for crops in previous studies) to explore the variation in crop yield in the 9 different prefectures of Guizhou Province. Moreover, we calculated Remote Sens. 2021, 13, 1614 8 of 23 the fertilization efficiency (FE) defined in Equation (6) to evaluate the comprehensive impact on both crop growth and the environment.
Here, Y is the yield per area for each selected crop species (6 in total), A represents the corresponding cultivated area. In a specific region, FE was calculated as total crop production divided by total N or P balance. A higher FE indicates the strategy of fertilization is more effective and reasonable, with the significant promotion of crop growth and low environmental cost.
The convergence test determines the final number (n) of BP networks brought into the ensemble method. The indicator for testing the convergence was set as the percentage of pixels with a 5-10% increase of yield after increasing N&P fertilizer and irrigation rates. The specific process of this method can be divided into four steps: (1) The number (n) of BP ANN was set from 1 to 40, respectively; (2) For each number of n, we built ensembled ANNs and ran the simulation 10 times; (3) After 10 simulations, we evaluated the performance of the ensemble by calculating the value of the defined indicator, correlation coefficient and time cost; (4) We compared the performance of ensembled ANNs to determine the final n for further simulations.
The optimization of n in the ensembled group can help us balance simulation stability with the time cost.

Simulation of Single BP Network
To build the ensemble model, we first examined the simulation of crop yield per area and N, P balance by a single BP network with the yield influencing parameters detailed in Methods. The simulation for N, P balance is better than that for crop yield, with correlation coefficients (between the true value and fitted value) larger than 0.99 ( Figure A1 in Appendix A). In contrast, the crop yield simulation accuracy varied among different species. For example, groundnut simulation had the best performance, indicated by the largest R (0.87) and least RME (11.67%), compared with the worst performance for potato (with R of 0.68 and RME of 15.38%). In terms of the time cost, the process of simulation of N or P balance is less than that of crop yield, about 2-3 (s) for training each network, while the time cost for simulation of crop yield ranges from 11.1 (s) to 19.7 (s) ( Figure A1).

Convergence Test for Ensembled Networks
As introduced in Section 2.4.4, we used the change of crop yield (with an interval of increase of 5-10%) after changing the investment of N&P fertilizer and irrigation as an indicator to test the stability and convergence of ensembled BP networks. The number of networks in each ensembled group ranged from 1 to 40, and each group was conducted 10 times to analyze the distribution of the simulation results (detailed in Methods). Considering the efficiency of the process, we set the epoch number to 5000 for each time. In Figures A2 and A3 we can see that when the number of networks is less than 10, the simulation results fluctuated greatly in terms of the indicator percentage.

Yield Gap in Guizhou Province
The extent of yield gap (or potential yield) varied among the different pairs of crop and prefecture, with some crop species having relatively more potential for yield improvement in some of the prefectures compared with others. For example, the total yield gap of potato in Bijie and Guiyang is 365,666.0 (t) and 17,665.9 (t), respectively ( Figure 5). Meanwhile, the ratio of actual yield and theoretical yield represents the proportion of local crop production that has achieved its maximum value. It can be seen that rice had the greatest ratio among the different crops, achieving more than 60% of the theoretical maximum yield per area and 72% for total yield in all the prefectures. By contrast, the yield per area of groundnut in Liupanshui and Qiandongnan only reached 47% and 52% of the theoretical maximum value. Lastly, the spatial difference of the ratio of maize is large, shown by the value of 88% in Bijie, but only 53% in Qiandongnan for the total production. This average percentage of achieved total production is all over 60% for each crop species in Guizhou, with a total mean value of 72.5% for all the selected crops.

Crop Yield and Fertilization Efficiency under Different Scenarios
After training the ensembled networks, we first increased the input variables of N&P and irrigation by 10% and predicted the variation of crop yield and fertilizer balance. On

Crop Yield and Fertilization Efficiency under Different Scenarios
After training the ensembled networks, we first increased the input variables of N&P and irrigation by 10% and predicted the variation of crop yield and fertilizer balance. On the whole, the change in yield was both crop-and prefecture-specific. Increasing fertilizer application benefited the yield (per area) of potato in most prefectures, except for Anshun, where potato yield decreased (Figure 6a). In particular, potato yield in Zunyi increased the greatest by 2.59 (t/ha). By contrast, other crops like maize in Zunyi, Bijie and Tongren decreased under the same scenario, although the degree of change was less than that of the increase in potato yield. The intensification of fertilization also caused the increase of N, P balance in all prefectures, with an average total increase in N of 8045.8 (kg/grid) and P of 1614.6 (kg/grid) in Guizhou. The increase in N and P balance was highly correlated in all prefectures ( Figure A7). Moreover, even though some crop yields increased in this scenario, the total FE in all prefectures decreased, except for Bijie (Figure 6c). On the other hand, implementing an increase in irrigation by 10% could favor the growth of crops in most cases (Figure 6b). The average rise of yield per area is 0.54 (t/ha) and 0.47 (t/ha) for potato and rice. Spatially, however, the impacts are highly heterogeneous, as shown by the significant yield increase in Qianxinan, Qiandongnan and Qiangnan, with an average increase of 0.064 (t/ha), 0.031 (t/ha) and 0.027 (t/ha) for all crops, whereas the yield increase in Guiyang, Anshun and Tongren is negligible (detailed distribution of yield change in each prefecture illustrated in Figures A4 and A5).  In the scenario of decreasing N&P fertilizer, the change in crop yield varied both spatially and for the species distribution (Figure 7a). Similarly to the former scenario, potato yield was the most impacted by changing the quantity of fertilization. Except for Qianxinan and Qiannan, the absolute yield decrease ranged from 1.94 t/ha (Tongren) to 0.28 t/ha (Zunyi) across Guizhou Province. By contrast, maize yield even showed a small increase following decreasing fertilizer investment, which is significant in mid-eastern prefectures like Zunyi, Tongren and Qiandongnan. Spatially, the greatest reduction in crop In the scenario of decreasing N&P fertilizer, the change in crop yield varied both spatially and for the species distribution (Figure 7a). Similarly to the former scenario, potato yield was the most impacted by changing the quantity of fertilization. Except for Qianxinan and Qiannan, the absolute yield decrease ranged from 1.94 t/ha (Tongren) to 0.28 t/ha (Zunyi) across Guizhou Province. By contrast, maize yield even showed a small increase following decreasing fertilizer investment, which is significant in mid-eastern prefectures like Zunyi, Tongren and Qiandongnan. Spatially, the greatest reduction in crop yield following fertilization reduction was shown in some southwestern prefectures (e.g., Anshun and Liupanshui). The decrease of all crop yields in Anshun ranged from 0.10 t/ha (groundnut) to 0.34 t/ha (potato). Moreover, the N, P balance also showed a synchronous shift after decreasing fertilization. The decline of N balance ranged from 13,605.0 kg/grid (Guiyang) to 4736.7 kg/grid (Qiandongnan), while P balance decreased much less (from 2894.3 kg/grid in Guiyang to 948.8 kg/grid in Qiandongnan). This phenomenon was following the former scenario. Moreover, it is worth noting that there is an evident spatial coherency for the change in N, P balance in Guizhou Province (R > 0.99), which implied that these shifts were mainly determined by the local environmental conditions, under the same change of fertilization. Lastly, this strategy demonstrated strong optimization on FE, shown by an increased value for all prefectures, except for Liupanshui, which had no significant variation in this scenario (Figure 7b; detailed distribution of yield change in each prefecture was illustrated in Figure A6).

Yield Gap Closing
Finally we calculated the ratio of (total) crop yield variation and (total) yield gap obtained from the ensemble networks and investigated the influence of intensifying management strategies (N&P fertilizer and irrigation combined) on closing the existing local yield gap. The impacts on the yield gap were both species and spatially heterogeneous (Figure 8). From a species perspective, potato and soybean gained the most crop yield, reducing the yield gap by an average of 33.9% and 16.6%, respectively. However, the yields of groundnut and maize decreased following additional agricultural investment, with an average reduction of −7.6% and −7.4%. From a spatial perspective, most crops in southern parts of Guizhou Province in Qiangdongnan, Qianxinan and Qiannan benefited from closing the yield gap by increasing N&P fertilizer and irrigation, where the averaged ratio value is 17.2%, 15.0% and 4.0%, respectively. At the same time, other prefectures had more heterogeneity in their yield gap variation.

Yield Gap Closing
Finally we calculated the ratio of (total) crop yield variation and (total) yield gap obtained from the ensemble networks and investigated the influence of intensifying management strategies (N&P fertilizer and irrigation combined) on closing the existing local yield gap. The impacts on the yield gap were both species and spatially heterogeneous (Figure 8). From a species perspective, potato and soybean gained the most crop yield, reducing the yield gap by an average of 33.9% and 16.6%, respectively. However, the yields of groundnut and maize decreased following additional agricultural investment, with an average reduction of −7.6% and −7.4%. From a spatial perspective, most crops in southern parts of Guizhou Province in Qiangdongnan, Qianxinan and Qiannan benefited from closing the yield gap by increasing N&P fertilizer and irrigation, where the averaged ratio value is 17.2%, 15.0% and 4.0%, respectively. At the same time, other prefectures had more heterogeneity in their yield gap variation. reducing the yield gap by an average of 33.9% and 16.6%, respectively. However, the yields of groundnut and maize decreased following additional agricultural investment, with an average reduction of −7.6% and −7.4%. From a spatial perspective, most crops in southern parts of Guizhou Province in Qiangdongnan, Qianxinan and Qiannan benefited from closing the yield gap by increasing N&P fertilizer and irrigation, where the averaged ratio value is 17.2%, 15.0% and 4.0%, respectively. At the same time, other prefectures had more heterogeneity in their yield gap variation. Figure 8. Percent of yield gap closing by increasing N&P fertilizer and irrigation by 10%. Figure 8. Percent of yield gap closing by increasing N&P fertilizer and irrigation by 10%.

Discussion
Establishing the quantitative relationship between crop yield and the environment is essential for understanding each factor's interaction and predicting future crop yield under different climate or agricultural scenarios. Previously, these demands were mainly met by conducting different regression models or process-based models [37,38,75]. Nevertheless, these approaches have three issues. First, the performance of traditional models relied on the characteristics of the data distribution. For example, higher accuracy can be obtained when the crop yield variation is large [76]. Second, the large differences occurred between existing model outputs for simulating crop yield, and the relative error in some cases surpassed 50% [35]. Third, the reported accuracy for traditional regression and process-based models could only explain 51% and 42% of observed yield variability, respectively [77]. The technology of machine learning, including ANN, can counter these deficiencies to some extent. However, in previous studies, researchers have been cautious about predictions of crop yield under different scenarios by trained networks [31,37,45]. Most of these studies were constrained to fitting the existing crop yield; thus the work of prediction was mainly finished by existing traditional models [78][79][80]. Utilizing the characteristics of convergence of ensemble networks, we found that the usage of a single network may generate issues in terms of the stability of the result, revealed by the fluctuation of correlation coefficient of actual yield and fitted value. Namely, the correlation coefficient for each simulation can be high, but the distribution of the fitted value may vary dramatically, which was demonstrated by the indicator value fluctuating (Figures A2 and A3). This issue could be prominent, especially when we conducted the prediction by a single network. Clustering of multiple networks (forming ensembled networks) can resolve these issues, providing a more reliable forecast of crop yield under future scenarios, but the time cost should also be fully considered.
To date, two main approaches are commonly used for estimating crop yield gaps in specific locations. The first is to apply crop models to set optimal environmental parameters to calculate the potential yield and gain the difference between the result with actual yield [81,82]. This method requires comprehensive environmental parameters throughout the growing season for each crop (often with daily temporal resolution), relying on the availability of historical data records. In contrast, the second method, which we employed in this study, calculates the yield gap by statistical data. This method is simple to conduct and fully considers the actual yield to obtain the locally exploitable yield gap. However, previous studies utilizing this method only imported variables of temperature elements to classify natural zones [52]. To strengthen our analysis, we added shortwave radiation, rainfall and slope as well in this study. Particularly, slope (representing topography and, therefore, areas of severe soil degradation in the region) has been proven to be the main factor for limiting the absorption of water and nutrient resources belowground for vegetation (including crop) growth in the karst region [9,83]. Therefore, this approach enables our subsequent analyses and conclusions for agricultural management to be more reliable.
In the past, large-scale yield gap studies have demonstrated a big difference between the yield gap in developed and developing countries, which means that some countries like European countries have achieved a high percentage of maximum yield in the local environment, while others have more scope to further increase local crop production [34]. This phenomenon was mainly caused by the imbalance of economic development and the spread of new agricultural techniques among different countries. For example, in Nebraska, USA, the production of maize has reached nearly 90% of the theoretical maximum, while in some African countries, this value is even less than 50% [84,85]. By contrast, the yield (per area) gap of maize in Guizhou was roughly ranked in between these two groups, ranging from 46% in Qiandongnan to 89% in Bijie. Meanwhile, the growth in crop yield in some areas has slowed in recent years, like maize in African countries and rice in the southern part of China [29]. Therefore, closing the yield gap by different approaches in these regions requires further investigation, especially in less developed areas [85].
For most crop species (excepting potato) in Guizhou Province, we found that the effect of fertilization has reached a tipping point, whereby increasing inputs of fertilizer alone will not achieve a significant rise in yield. Indeed, for some crops, additional fertilizer investment could even cause a reduction in crop yield (Figure 6a). Three reasons can explain this. First, the absorption of nutrient elements for crops like wheat has reached "saturation" in the soil, resulting in a low additional increase of absorption being detected after the fertilization tipping point is reached [86]. Second, the facilitation of crop growth via fertilization requires reasonable consideration of the ratios of different fertilizers (N, P and K) specific to the local area's nutrient limitations and the crop species physiology, rather than a prescribed addition of each [87]. Third, the overuse of inorganic fertilizers can lead to negative environmental consequences, indicated by hardening of the soil, decreasing fertility, strengthening pesticides, polluting air and water resources, which further influences crop growth and production [88]. Therefore, more fertilization on crops in Guizhou Province could only cause more fertilizer balance left in the environment and less fertilization efficiency, further resulting in more environmental problems, which is illustrated by Figure 6c. According to statistics, average chemical fertilizer (N, P, and K) and pesticide use per hectare of cropland in China are between double and three times that of other countries/regions in the world [89]. Field surplus nitrogen (N) and farm disposal N are major sources of water pollution in farming systems, which contaminates the water supply for livelihood [90]. Hence, it is very important to manage the investment of different agricultural chemicals more precisely and more rationally. Our results also demonstrated that an appropriate reduction in fertilizer use might benefit the crop yield in some regions, as well as causing an evident increase in fertilization efficiency. For example, higher maize production in some mid-eastern prefectures existed with a distinct decrease in the N, P balance in Guizhou Province. However, the implementation of this strategy should fully consider the spatial heterogeneity. Among all the prefectures, Qianxinan was highly recommended to reduce fertilization, resulting in slightly higher crop yield, reduced fertilizer balance and agricultural management cost savings.
Furthermore, an appropriate increase of irrigation is able to favor crop yield in Guizhou Province, aside from the economic and natural resources costs involved (Figure 6b). The lack of water resources for crops in Guizhou karst farmland was also shown by the high crop water stress index for irrigation due to seasonal rainfall shortage and aggravated by the thin soil layer and poor water holding capacity in this area [91]. Although the impact of changing irrigation was less significant than that of fertilization (e.g., the increase of yield for potato by increasing irrigation is less than that from increasing fertilization), it has a more stable and positive consequence for most species in most regions. Compared with irrigation rates in the North China Plain and other northern farming areas, rainfed farmland accounted for a high percentage in southwest China. This region may also encounter additional water scarcity under global warming [92]. Therefore, considering the local topography, fragmentation of cropland, and other environmental characteristics, the combination of both "soft-path" (capturing water resources in small and check dams) and "hard-path" (with large, centralized, capital intensive irrigation projects and water storage infrastructure) may be necessary to enable sufficient irrigation water to meet the demand for crop growth in the future [93]. Moreover, it is notable that the implementation of these development patterns is often limited by "economic water scarcity" resulting from the cost of building irrigation infrastructure under local financial conditions [94].
Together, the continuous increase of fertilization and irrigation have a relatively limited effect on improving crop yield in Guizhou Province, except for potato and soybean. These indications can be used as a DST to target agricultural management to the areas where it will be most effective, and the stakeholders should consider other strategies to meet the local food demand. For example, optimizing plant developmental features, crossbreeding and hybridization, new genetic techniques and ecological engineering strategies, such as increasing biodiversity and soil fertility in farmland to harmonize agricultural development and the environment [95][96][97].
Other than the scale effect and errors brought by different spatiotemporal remote sensing datasets, this study's uncertainty was also caused by the "extrapolation" effect of ANNs. Indeed, just like traditional regression models, the simulation accuracy of a neural network also relies on the distribution of training and input data for simulation. Where there is limited training data (e.g., the input data for prediction is beyond the range of training), the corresponding trained network tends to have an arbitrary prediction or classification [98], which can reduce the accuracy of the predictions. Different from the traditional classification of computer graphics or remotely sensed images [99,100], we often set future scenarios for environmental parameters (like increasing fertilization and irrigation in this study or higher global temperature) beyond that found in historical records, which may put the data out of the range of the available training data in the region and further impact the accuracy of the predictions. This error may also be expanded if there are large differences between the data utilized for prediction and the available training data [101]. Therefore, we only set the absolute rate of change at 10% for fertilization and irrigation.

Conclusions
In this paper, we calculated the yield gap for 6 crop species grown in Guizhou Province and used ensembled ANNs to predict the change in crop yield in Guizhou under different management scenarios of increasing/decreasing fertilization and irrigation. The conclusions are as follows: (1) The ensemble of ANNs can improve the robustness of simulation of crop yield by multiple input factors, which is especially beneficial for predictions under different future scenarios; (2) The selected crops' total yield has already realized over 60% for each of the theoretical maximum, with an average value of 72.5% in all cases. The yield gap in the study area is still larger than that in most developed countries; (3) An appropriate increase in irrigation can benefit the crop yield for most Guizhou species, compared with increasing fertilization. Moderate reduction of fertilizer in this region may increase some crop production and enhance local fertilization efficiency.
Combing the management strategies of increasing fertilization and irrigation was beneficial for increasing the yield of potato and soybean, but this approach has a limited effect on the existing yield gap for most other crops in Guizhou Province.

Funding:
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was funded by the National Natural Science Foundation of China under grant no. 41571130044, and the scholarship of the China Scholarship Council. The UK team was supported by the Natural Environmental Research Council MIDST-CZO project (-NE/S009175/1).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data utilized in the study is available from the authors, on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.         Figure A6. Distribution of change of yield (t/ha) under the scenario of decreasing fertilizer.