Comparative Study on Spatial Digital Mapping Methods of Soil Nutrients Based on Different Geospatial Technologies

: Soil organic matter (SOM), total nitrogen (TN), available phosphorus (AP), and available potassium (AK) are important indicators of soil fertility when undertaking a quality evaluation. Obtaining a high-precision spatial distribution map of soil nutrients is of great signiﬁcance for the differentiated management of nutrient resources and reducing non-point source pollution. How-ever, the spatial heterogeneity of soil nutrients lead to uncertainty in the modeling process. To determine the best interpolation method, terrain, climate, and vegetation factors were used as auxiliary variables to participate in the investigation of soil nutrient spatial modeling in the present study. We used the mean error (ME), mean absolute error (MAE), root mean square error (RMSE), and accuracy (Acc) of a dataset to comprehensively compare the performance of four different geospatial techniques: ordinary kriging (OK), regression kriging (RK), geographically weighted regression kriging (GWRK), and multiscale geographically weighted regression kriging (MGWRK). The results showed that the hybrid methods (RK, GWRK, and MGWRK) could improve the prediction accuracy to a certain extent when the residuals were spatially correlated; however, this improvement was not signiﬁcant. The new MGWRK model has certain advantages in reducing the overall residual level, but it failed to achieve the desired accuracy. Considering the cost of modeling, the OK method still provides an interpolation method with a relatively simple analysis process and relatively reliable results. Therefore, it may be more beneﬁcial to design soil sampling rationally and obtain higher-quality auxiliary variable data than to seek complex statistical methods to improve spatial prediction accuracy. This research provides a reference for the spatial mapping of soil nutrients at the farmland scale.


Introduction
Soil nutrient indicators are key indicators for the evaluation of soil quality. Studying the spatial distribution of soil nutrients is the basis for understanding regional soil quality conditions, adjusting management measures and various material inputs, and obtaining maximum benefits [1][2][3]. However, due to the combined effect of structure and randomness [4,5], soil nutrients have a high degree of spatial heterogeneity and dependence [6]. Therefore, it is necessary but difficult to obtain more accurate spatial distribution information of soil nutrients in different regions.
In the past few decades, scholars have developed and applied many spatial interpolation methods, including deterministic interpolation methods (e.g., inverse distance weighted (IDW), radial basis function (RBF), global polynomial interpolation, and local polynomial interpolation), geostatistical methods (e.g., ordinary kriging (OK), simple kriging, and CoKriging (COK)), and hybrid techniques (e.g., regression kriging (RK) and geographically weighted regression kriging (GWRK)) [7]. Among these methods, the OK method is the most widely used geostatistical method. Through intensive field sampling, OK makes full use of the data information of the sample point, and can provide an estimation error for the interpolation results [8]. However, OK would consume large amounts of labor, material, and financial resources [9]. The sampling density and method have a considerable influence on the interpolation accuracy of OK [10]. Moreover, this method does not incorporate environmental factors into the model, which are closely related to soil nutrients. These limitations lead to the simulation of the spatial distribution of soil properties in complex landscapes being greatly restricted [11,12].
With the rapid development of geographic information technology and remote sensing technology, other environmental factors that have the advantages of high accuracy and easy access have gradually been incorporated into spatial modeling analyses. Many studies have shown that methods that use environmental covariates (e.g., RK) usually produce more accurate simulations than OK. Watt and Palmer [13] used RK to draw a surface map of New Zealand's carbon/nitrogen ratio. Bangroo et al. [14] used RK to study the influence of topographic factors on the spatial distribution of soil organic carbon (SOC) and total soil nitrogen (TSN); the authors found that the RK method was significantly better than the OK model in predicting SOC and TSN in the case of residuals with a moderate degree of spatial autocorrelation. Mondal et al. [15] selected eight variables to predict SOC; the RK method also provided satisfactory results from the accuracy point of view due to the introduction of additional auxiliary information. The RK method combines the trend item and kriging residual item generated by global regression fitting [16,17]. Thus, it has the advantages of easy implementation, detailed results, and high prediction accuracy [18].
However, the relationships between soil nutrients and environmental covariates are always non-stationary over space, but the RK method assumes that the relationship between soil properties and environmental covariates is constant globally [12,19]. As a result, RK cannot capture the local characteristics of the spatial variation of soil nutrients, and the improvement of interpolation accuracy is again restricted. To solve this problem, the local model GWRK has been applied. As a combination of geographically weighted regression (GWR) and kriging, this method embeds the spatial position into the regression parameters through the idea of local weighted least squares [20,21]. It not only considers the non-stationarity of the spatial parameters and the local relationship between the target variable and the explanatory variable, but also the spatial autocorrelation of the regression residuals [22]. The GWRK method provides good prediction results in most cases [23][24][25].
However, GWR also has certain limitations. The method uses a single bandwidth for all variables, ignoring the different scales of the relationship between the independent variable and the dependent variable. This may cause a lot of noise and bias in the regression coefficient, which may lead to unstable regression coefficients [26]. Nevertheless, scale is an important term in geographic research [27,28], and spatial phenomenon are inherently affected by scale effects [29]. Therefore, the change in scale has a considerable impact on the model performance [30,31]. Murakami et al. [32] also emphasized the importance of scale in local regression modeling. The first model to consider multi-scale effects was the hybrid GWR model. The regression coefficients of some parameters in this model are globally fixed, while the remaining parameters are spatially variable [33][34][35][36]. However, the hybrid GWR is limited; it can only distinguish the scale influence of different variables into global and local scales, but the local coefficients of the variables still change within the same scale [37]. Fotheringham et al. [38] proposed the multi-scale geographic weighted regression (MGWR) method in 2017 to allow multi-scale modeling, which allows the model to run at different spatial scales by introducing a specific bandwidth for each variable. As an extension of MGWR, MGWRK [39] is similar to other hybrid methods. It not only retains the characteristics of MGWR, but also accounts for the relevance of the regression residuals; thus, it appears to be a promising method.
In summary, scholars have carried out a lot of research on the digital mapping of soil nutrients. Although the emergence of different methods is necessary, the portability of these methods is generally poor. The selected explanatory variables differ between research areas and for various scales, the best interpolation method may also differ between target variables. At the same time, as an emerging model, few studies have reported on the application of MGWRK at the farmland scale. Therefore, the main aims of the present study are to select environmental factors such as, vegetation, terrain, and climate factors for modeling, and to compare the performance of different geospatial technologies (OK, RK, GWRK, and MGWRK) in soil nutrient mapping.

Study Area
Gaoping City is located at the foothills of the Taihang Mountains in the southeast area of Shanxi Province in central China (35 • 40 -36 • 0 N, 112 • 40 -113 • 10 E; Figure 1). The total land area of the city is 946 km 2 , with an altitude ranging from 780 m to 1391 m. Gaoping is surrounded by mountains in the east, west, and north, with flat land in the middle. There are many types of land use and large areas of arable land. The main soil types are cinnamon soil and fluvo-aquic soil. The study area has a warm temperate continental monsoon climate. The annual average temperature is approximately 10.5 • C, the annual average sunshine duration is 2532.5 h, the frost-free period is 180-200 days, and the annual average rainfall is about 600 mm (mainly from July to September). Shanxi Province, Gaoping, where the research was carried out is characterized by good water and air quality and as such offers optimum research possibilities.

Soil Sampling and Determination
A total of 2583 samples were collected before crops were planted in the spring of 2010, in accordance with the "Rules for soil quality survey and assessment" (NY/T 1634-2008). We used the "S" method to layout points in each field by employing stainless steel soil drills and other tools. The spatial distribution of the data points is shown in Figure 1. After pretreatment of the soil samples, the contents of soil organic matter (SOM), total nitrogen (TN), available phosphorus (AP), and available potassium (AK) of each sample were determined. Specific determination methods are described elsewhere [40].

Environmental Data
Numerous studies have demonstrated that soil nutrient content is comprehensively affected by factors such as climate, vegetation, and terrain [3,14,41,42]. Therefore, this paper selected these factors as auxiliary variables to participate in the investigation of soil nutrient spatial modeling.
The climate factors included mean annual precipitation (MAP) and mean annual temperature (MAT), which were obtained from National Meteorological Information Center [43]. We downloaded the precipitation and temperature data of local meteorological stations from 1980 to 2010. By calculating the average multi-year precipitation and average multi-year temperature, and by performing kriging interpolation, meteorological data at each soil sampling point in Gaoping City were extracted.
The topographic factors were derived from the DEM (30 m spatial resolution) obtained from Geospatial Data Cloud [44]. Based on the DEM data, we obtained various composite terrain factor indicators in the SAGA-GIS software.
The remote-sensing data include the NDVI and vegetation coverage (VC), which were also derived from Geospatial Data Cloud [44]. We downloaded Landsat 8 remote sensing images that were consistent with the sampling time. After preprocessing the image, we used the band calculator to calculate the NDVI and VC to reflect the growth status of farmland vegetation in that year. The data of various environmental covariates are shown in Table 1. The OK model is based on spatial autocorrelation and second-order stationary assumptions. It uses the semivariogram theory as a tool to analyze its spatial structure, and realizes the optimal unbiased linear estimation of the unknown area based on the existing sampling point data in a limited area [11]. The semivariance function equation is described as Equation (1).
where γ(h) is the semivariogram value; h is the lag; N(h) is the total number of point pairs with an interval of h; Z(x i ) and Z(x i + h) are the measured values of regionalized variables at positions x i and x i + h, respectively.

Regression Kriging (RK)
Regression kriging combines multiple linear regression (MLR) and kriging interpolation. The steps for establishing this method are as follows: a stepwise MLR analysis is performed between the target variable and its highly correlated auxiliary variables to obtain the trend item and residual item, which represent certainty and randomness, respectively. Then, OK is applied to interpolate the residual item of the regression model before finally summing them up [16,18]. The fundamental equation of the RK model is described as Equation (2): where y RK (x i ) denotes the predicted value of the RK model at the point is the residual item estimated by the OK method at the point x i .

Geographically Weighted Regression Kriging (GWRK)
As a hybrid method, the modeling process of GWRK [25] is similar to that of RK, with the exception of the fact that the global fitting in RK is replaced with the local fitting in GWR. The GWR approach is used to deal with spatial non-stationarity by introducing a spatial location into regression coefficients [21].
Assuming that there is a total of n observation points, the location of each observation point is (u i , v i ), and there is a total of m independent variables involved in modeling. The expression of the GWR model is shown in Equation (3): where x ij is the jth independent variable at the observation point; Two types of spatial weight function methods, the Gaussian function and Bi-square function, are usually used to determine the spatial weight in GWR analysis. Different types of kernel functions can be used during operation, including fixed kernel functions and adaptive kernel functions. Research has shown that an adaptive kernel is arguably more favorable when dealing with non-uniform spatial distributions [21,24]. Because the sampling points are not uniformly distributed, this study selects an adaptive Gaussian kernel function for modeling. In addition, geographic weighted regression analysis is very sensitive to the bandwidth selection of a specific weight function [45]. This study uses the corrected Akaike information criterion (AICc) to determine the model bandwidth, where the principle is to minimize the AICc. This approach helps to evaluate whether the GWR model simulates the data better than the MLR model.

Multiscale Geographically Weighted Regression Kriging (MGWRK)
As mentioned, MGWRK [39] is a combination of MGWR and kriging interpolation. The MGWR model can be expressed as Equation (4): where β bwj is the regression coefficient corrected by the effective bandwidth of the jth independent variable. In our study, the kernel functions and kernel types of MGWR and GWR models are consistent. The MGWR model accounts for the spatial scale effect of the influence of different independent variables on the dependent variable, and provides different independent variables with different bandwidths. This is the main difference between the MGWR model and the GWR model. All independent variables in the GWR model share a bandwidth, which can be regarded as a weighted average of different levels of spatial heterogeneity [38,46].

Model Validation and Evaluation
Ninety percent of sampling points (2284) were randomly selected from the original data as the prediction set to build the model, and the remaining 10% samples (254) were used as the verification set to test the model's prediction accuracy. We evaluated the performance of different models by calculating the mean error (ME), mean absolute error (MAE), root mean square error (RMSE), and prediction accuracy (Acc) of the validation set [33,47]. These four indices can be written as Equations (5)-(8): whereẐ(x i ) is the predicted value at position x i ; Z(x i ) is the measured value at position x i ; n is the number of samples.
In general, ME, MAE, and RMSE values closer to zero and a higher Acc give a better prediction.

Descriptive Statistics
In this study, 11 indicators (elevation, slope, aspect, PC, TR, SR, TWI, MAP, MAT, NDVI, and VC) were initially selected as environmental variables ( Figure 2).  Table 2 shows the descriptive statistics of the soil nutrients. The SOM content ranged from 3.10 g/kg to 45 g/kg. The TN, AP, and AK contents of soil ranged from 0.11 g/kg to 2.51 g/kg, 1.40 mg/kg to 39.90 mg/kg, and 38 mg/kg to 388 mg/kg, respectively. The coefficient of variation (CV) is defined as the ratio of the standard deviation to the mean, and is a normalized measure of the dispersion of a probability distribution or frequency distribu-tion [20]. The CV values of the SOM, TN, AP, and AK contents were 26.92%, 21.65%, 48.48%, and 28.95%, respectively. Thus, the CV of the AP content was the largest, which is consistent with previous research results [48]. In addition, Table 3 shows the descriptive statistics of the environmental covariates. The CV of different environmental covariates in the study area also varied significantly.

Correlation Analysis of Soil Nutrients and Environmental Factors
The results of a Pearson's correlation analysis of soil nutrients and environmental factors in the study area are presented in Table 4. There were significant correlations between the SOM content and (i) topographic factors (elevation, slope, SR, TR), (ii) climate factors (MAP), and (iii) vegetation factors (NDVI and VC). The significant influencing variables of the TN content of soil were elevation, TR, and the NDVI. Both MAP and MAT had significant effects on the AP and AK contents. The AP content was also significantly negatively correlated with the TWI.
Terrain conditions affect soil nutrients by regulating the distribution of soil moisture as well as heat and air permeability conditions in a landscape [49,50]. In this study, topographic factors and soil nutrient content were significantly negatively correlated, thus indicating that environments with a high elevation, steep slopes, a fragmented ground surface, or low TWI usually have a low soil nutrient content.
Climate is the main driving force affecting plant type, crop growth, and litter decomposition [51], but soil nutrients have different response mechanisms to the different synergy of temperature and precipitation. Generally speaking, an environment with high precipitation and low temperature is conducive to the accumulation of soil nutrients [49]; however, an increase in precipitation may also increase soil erosion, which leads to the loss of effective soil elements [52]. Studies have also shown that temperature changes will affect the decomposition of microorganisms in the soil, thus complicating the relationship between temperature and nutrients [53]. In our study, the MAP had a positive effect on the SOM content, whereas it had a negative effect on AP and AK. Temperature exhibited a significant positive relationship with AP and AK, while the relationships between SOM, TN, and temperature were not significant.
Vegetation is also an important factor affecting the soil nutrient content. The influence of vegetation is usually mediated through terrain and climate variables [54]. Significant negative correlations were observed between the various vegetation factors and the soil nutrient content in Gaoping City, which is inconsistent with the results of Bangroo et al. [14]. It may be that the part of the plant residues returned to the soil by complex vegetation types hindered the accumulation of soil nutrients.

Selection of Environmental Covariates and MLR Modeling
Using a stepwise regression method to select variables can ensure that only independent variable factors that are significant to the soil nutrient content enter the regression model, and can also check the multicollinearity of the model [55]. The probability level of variable entry in stepwise regression is 0.05, and the probability level of elimination is 0.1 The final fitting model of the MLR of each dependent variable is shown in Table 5. Notes: SOM = soil organic matter, TN = total nitrogen, AP = available phosphorus, AK = available potassium; MAP = mean average precipitation, MAT = mean average temperature, SR = surface roughness, TR = topographic relief, PC = plane curvature, TWI = topographic wetness index, NDVI = normalized differential vegetation index, VC = vegetation coverage.
To further explore the influence of environmental covariates on the soil nutrient content in the study area at different locations, and to compare the differences in spatial mapping of different methods, the factors of the MLR model were also applied for GWR and MGWR modeling.

Spatial Variability Characteristics of Soil Nutrients and Residuals of MLR, GWR, and MGWR Models
The residuals of the regression models were used to fit the semivariogram model. We selected the optimal semivariogram model according to the maximum determination coefficient (R 2 ). Table 6 showed that the R 2 value of each model was >80%, which indicated that a good model simulation effect was achieved in each case. Semivariograms of soil nutrients and the residuals of MLR, GWR, and MGWR are shown in Figure 3.  The nugget/sill ratio of the optimal residual model was between 25% and 75%, which indicates that the model residual had a moderate degree of spatial autocorrelation. Therefore, it was necessary to use hybrid technologies (RK, GWRK, and MGWRK) for modeling [33]. The obtained semivariogram and various parameters were used to perform the kriging interpolation on the residuals. The results are presented in Figure 4, which shows that the three models produced very different spatial distributions of the residuals: (i) The extremum of the MLR residuals were the largest compared with other models, and the high and low values were clearly clustered, (ii) the spatial distribution of the GWR residuals presented a "low-high-staggered" phenomenon. (iii) the extremum of the MGWR residuals were the lowest for some indices, but the phenomenon of "lowhigh-staggered" distribution was much more serious than that of the GWR residuals. This may be because different bandwidths for various independent variables were used in the MGWR model, while all independent variables in the GWR model shared one mean bandwidth. Hence, the "fragmentation" phenomenon of the spatial distribution of GWR residuals was much less intense than that of the MGWR residuals.

Results of Soil Nutrient Mapping Based on Different Models
Based on the above analysis, we drew spatial distribution maps of the soil nutrient content in the study area based on different geospatial technologies (OK, RK, GWRK, MG-WRK) ( Figure 5). A comparison of the results revealed that the overall spatial distribution trends of nutrient indicators obtained by these different models were consistent. The spatial prediction effect of the OK method was smoother. The mapping effects of RK and GWRK were very similar, while the "fragmentation" phenomenon of the MGWRK mapping effect was more prominent, which was basically consistent with the performance characteristics of the residual spatial distribution map. The high-value areas of the SOM content in the study area were concentrated in the east, whereas the low-value areas were concentrated in the west. The TN content of soil decreased from the center of the study area toward the surrounding regions, which was contrary to the distribution patterns of various topographical factors. This implies that the TN content of soil in the study area is greatly affected by topographical factors. The AP content of soil in the western region of the study area was higher than that in the east, which was contrary to the characteristics of the spatial distribution of SOM. This may be because the regional differences in the amount of fertilizer applied, along with the fixed and difficult-to-move nature of phosphate fertilizers complicate the distribution of soil phosphorus over a small area [56,57]. The AK content was obviously higher in the southeast region of the study area compared to the northwest. The study also found that the soil nutrient content in the river valley is relatively high. This may be due to the flat terrain and convenient irrigation conditions, which is beneficial to the accumulation of various soil nutrient contents while soil nutrients are more likely to be lost in mountainous and hilly areas. At the same time, human activities such as different main crops and different management measures in various towns and small farm families may also cause differences in the spatial pattern of soil nutrient content.

Evaluation of Model Accuracy
The prediction accuracy of the OK, RK, GWRK, and MGWRK models for soil nutrients was comprehensively evaluated by comparing each model's performance based on the ME, MAE, RMSE, and Acc of 254 samples (Table 7). MLR, GWR, and MGWR participated in model accuracy evaluation simultaneously. The results showed that, for different soil properties, the best interpolation methods also differed depending on the selected auxiliary factors. We concluded that the best prediction model for (i) SOM was RK, (ii) TN and AP was RK or OK, and (iii) AK was GWRK.

Spatial Non-Stationary Relationship between Soil Nutrients and Environmental Variables
In actual situations, the relationships between environmental variables and soil nutrients may vary geographically, i.e., spatial non-stationarity. Liu et al. [20] demonstrated in their study that nearest distance to residential area (TRA), Normalized Difference Moisture Index (NDMI), slope, Normalized Difference Vegetation Index (NDVI), nearest distance to road (TRD), and elevation were the dominant variables affecting the spatial distribution of SOC stocks; while the other three factors (i.e., nearest distance to river (TRR), aspect, and the land cover degree comprehensive index (LDCI)) did not play a dominant role in any part of the study region.
SOM is a critical indicator for the evaluation of soil quality. Thus, this study took SOM as an example to compare and analyze the regression effects of MLR, GWR and MGWR. Table 8 showed that the GWR and MGWR models based on local regression idea could better explore and explain the spatial non-stationarity characteristics of variables than the global model MLR. It could also be seen from Table 7 that the MGWR model performed better than the GWR model in various diagnostic indicators. For example, the AICc value is reduced from 6947.297 (GWR) to 6823.630 (MGWR), and the R 2 value is increased from 0.165 (GWR) to 0.201 (MGWR), indicating that the MGWR model achieved a higher fitting accuracy. The maps of standardized local coefficients for the different explanatory variables by the MGWR model were provided in Figure 5, which could be a reference for the contribution of each influencing factor to SOM in various locations. The positive regression coefficients represent positive correlation, the negative regression coefficients correspond to negative correlation, and the absolute coefficient values reflect the degree of the influence of the factor on SOM [58]. Specifically, Figure 6b showed the standardized local coefficients of MAP (−2.39 to 3.86). MAP was negatively correlated with SOM in the northern region of the study area, and with the increasing MAP, the negative effect between MAP and SOM also increased. Whereas positive correlation was presented in the south, which indicated that SOM content increased with rising MAP over this area. Figure 6c showed the negative correlation between elevation and SOM in the western region, and the positive correlation in the eastern region. The local coefficients of slope ( Figure 6d) and SR (Figure 6e) showed similar trends across the space. Both of them varied from positive to strongly negative from northwest to southeast, which demonstrated that with the increase of slope and MAP in the northwest of the study area, higher SOM could be observed, while the southeast is the opposite. In terms of TR (Figure 6f), negative correlation between SOM and TR was found in the north and southwest, while positive correlation was seen in the southeast. The absolute coefficient value was the lowest in the central area, which indicated that the central area was less effected by TR. Such spatial patterns could be due to the distribution of TR (Figure 2g); the terrain of the central region is very flat, therefore the weak influence from TR. In Figure 6g,h, it could be seen that the two explanatory variables (PC and TWI) had a strong negative impact on SOM in the eastern part of the study area and positive impact in the west. The negative effect in the eastern fringe was the greatest, which indicated that PC and TWI play a critical role in the distribution of SOM here. The local estimates of VC coefficients (Figure 6i) varied from strongly negative to strongly positive from northwest to southeast. We also found that the impact of VC on SOM is very weak in the central area, where the main land use type was construction land with low vegetation coverage. Understanding the characteristics of spatial non-stationarity is helpful to further explore the spatial relationships between environmental variables and soil nutrients.

Comparison of Model Performance in Soil Nutrient Mapping
Accuracy indicators were used to prepare radar charts to visualize the spatial prediction accuracy of the different models (Figure 7). The prediction accuracy of the MLR was not very good, and this approach provided the worst predictions of SOM and TN. Because MLR is less flexible, it can only estimate soil nutrient content by fitting a global function. If the relationship between soil properties and covariates is spatially variable rather than globally fixed, MLR usually cannot simulate the distribution of soil nutrients accurately. Although the regression equations and regression coefficients of the four nutrient indicators passed the significance test, the low R 2 values proved that each regression equation could only explain a small part of the spatial variation of the corresponding data [59]. This may be due to the non-linear relationship between soil nutrient indicators and environmental parameters in Gaoping City; thus, traditional linear regression models are restricted for the study area. However, by performing OK interpolation on the regression residuals, the prediction accuracy of the RK model was significantly improved, and similar conclusions have been drawn in previous studies [60][61][62][63].
The local regression model GWRK achieved a good prediction accuracy in most cases because it considered the non-stationarity of the spatial parameters, the relationship between each target variable and the explanatory variable, and the spatial distribution characteristics and trends of the residuals [12]. Brian et al. [23] used the GWRK method to predict the spatial pattern of soil pH at a depth of 15 cm. The results showed that, compared with the global model based on MLR, the local model based on GWR significantly improved the prediction accuracy of soil pH. Kumar et al. [24] estimated the SOC density in the Midwest of the United States based on the GWRK method. The authors proved that GWRK provides a useful visualization tool for SOC estimation on a regional scale when good correlations exist between SOC and environmental variables, and when the residual has spatial autocorrelation. These studies demonstrated that GWRK is a promising method for spatial predictions of soil attributes. As a newly proposed method, MGWRK has certain advantages in reducing the overall residual level. The R 2 of MGWR is also higher than that of GWR, but its prediction accuracy is not very high compared with other models. Only in the prediction of TN, MGWRK ranked in the top three. The MGWRK model accounts for the spatial scale effect of different environmental covariates, which is a useful exploration. However, different bandwidths also have an impact on the spatial distribution of the residuals, which may lead to the failure of the MGWRK method to achieve the desired ideal effect. This was also reported by Qiao Lei et al. [39], who selected 13 indices (i.e., aspect, slope, MAP, MAT, altitude, annual cross primary productivity of vegetation, annual evapotranspiration, topographic wetness index, plane curvature, stream power index, terrain position index, terrain ruggedness index, and the annual average NDVI) as the environmental covariates, and compared the capabilities of OK, RK, GWRK, and MGWRK on SOM prediction; their results showed that the spatial prediction accuracy of the MGWRK method reached 69.0% of the RK method, 71.7% of the OK method, and 71.2% of the GWRK method. Therefore, MGWRK can provide a method reference for the digital mapping of the soil nutrient content in the study region, but its characteristics and performance require further investigation.
In summary, the use of hybrid technologies (RK, GWRK, and MGWRK) could improve the spatial prediction accuracy of soil properties than pure regression models (MLR, GWR, and MGWR). Hybrid methods also had advantages compared with OK, thus indicating that taking heterogeneity and dependence into consideration is necessary when model residuals showed moderate or strong spatial autocorrelation. Other studies reported similar results [23,64,65]. Yang et al. [33] compared the performance of MLR, RK, GWR, GWRK, MGWR, and MGWRK in mapping topsoil electrical conductivity (EC). The authors also found that hybrid approaches (i.e., MGWRK, RK, and GWRK) occupied the first three positions in the digital mapping of soil EC. However, the improvement of hybrid methods is very limited in our study, which may be due to the fuzzy correlations between soil properties and auxiliary variables [66]. The OK method, which makes full use of the spatial structure information of the data, has also achieved a satisfactory accuracy at a high sampling density and with a strong spatial correlation of the sampled data. Therefore, it may be more beneficial to design soil sampling rationally and obtain higher-quality data of auxiliary variables than to look for complex statistical methods to improve spatial prediction accuracy [10,67].

Conclusions
To summarize, this study compared the performance of OK, MLR, RK, GWR, GWRK, MGWR, and MGWRK for the digital mapping of soil nutrients in Gaoping City, China. The results showed that these models performed differently when predicting various indicators. The hybrid methods (RK, GWRK, and MGWRK) improved model prediction accuracy to a certain extent when the residuals were spatially correlated; however, this improvement was not significant. The new method MGWRK has certain advantages in reducing the overall residual level, but it failed to achieve the desired accuracy. Further verification of the applicable conditions and deviations of the MGWRK approach is required. Considering the cost of modeling, the OK method still provided an interpolation method with a relatively simple analysis process and relatively reliable results. This does not mean that more complex models would lead to an improved model performance.
In addition, the prediction accuracy of a spatial regression model largely depends on the correct selection of auxiliary variables. Although all the variables selected in this study were natural factors, soil nutrients are simultaneously affected by five major soilforming factors and human factors. Human activities in low-altitude areas increase the spatial variability of soil nutrients. Future research should select appropriate methods to quantify management measures (e.g., planting systems, farming measurements), and then incorporate them into models to further improve prediction accuracy. Data Availability Statement: Publicly available datasets were analyzed in this study. These data can be found here: http://data.cma.cnr (acceded on 15 October 2020); http://www.gscloud.cn (acceded on 15 October 2020).