A machine learning and geostatistical hybrid method to improve spatial prediction accuracy of soil potentially toxic elements

Effective environmental management and contamination remediation require accurate spatial distributions and predictions of potentially toxic elements (PTEs) in the soil. However, no single method has been developed to predict soil PTE accurately. This study evaluated the advanced geostatistical method of empirical Bayesian kriging regression prediction (EBKRP), machine learning algorithm of random forest (RF), and hybridized model (RF-EBKRP) to predict and map PTEs content in greenspace soils. As identified by RF, soil organic carbon, organic matter, total (nitrogen, phosphorus, and potassium), topographic features, and urban functional types were used as significant covariates to improve the prediction accuracy of PTEs in soil. The model prediction performance was evaluated using the root mean square error (RMSE), mean absolute percentage error (MAPE), and coefficient of determination (R2). Results showed that RF performed much better than EBKRP in predicting soil PTEs, with lower prediction errors and a higher R2. The RMSE, MAPE, and R2 values for the RF model were 0.25–85.32 mg/kg, 3.86–25.40%, and 0.77–0.90, respectively, while the values for the EBKRP method were 0.51–99.03 mg/kg, 5.42–32.13%, and 0.40–0.66. Moreover, the RF-EBKRP method produced more accurate spatial predictions and distributions of PTEs than the individual models, with R2 improvements of 122.5% for EBKRP and 15.58% for RF. The better performances of RF-EBKRP are due to its incorporation of various covariates and its ability to handle complex nonlinear relationships between soil PTEs and covariates. In the end, a hybrid RF-EBKRP method is a promising approach to improving the spatial distribution map of soil PTEs.


Introduction
Urban soil contamination with potentially toxic elements (PTEs) such as Cd, As, Cr, Ni, Hg, Pb, and Zn has elicited significant global concern because of its potential environmental and human health threats (Praveena et al. 2015;Amari et al. 2017). The threat of PTEs to the environment and human health is due to their toxicity, persistence, bioaccumulation (Amari et al. 2017), easy enrichment, and poor mobility in the soil (Sundaramanickam et al. 2016). In addition, urban soils are enriched with a substantial level of PTEs relative to the natural reference values (Liu et al. 2006;Miao et al. 2008;Adedeji et al. 2019). The high levels and variations in urban soil pollution are primarily due to various PTE sources in the cities. It can be categorized as point source or non-point source pollution. Point source pollution, such as from vehicles and ignition of fuels (Minguillón et al. 2014), emissions from factories, industries, power plants, and transportation (Dao et al. 2014;González-Guzmán et al. 2022), and the non-point sources related to rainfall or stormwater runoff, urban sprawl (Martínez and Poleto 2014;Walaszek et al. 2018;Zhi et al. 2018) contribute to the enrichment of PTEs in urban soils. City expansion or urbanization also enhances the levels of PTEs in the soil due to the increased area of urban impervious surfaces, such as paved roads, urban sprawl, roofs, and exterior walls of residential and industrial facilities, on which various contaminants accumulate (Zhi et al. 2018;Jeong et al. 2020). Moreover, rapid industrialization, urbanization, population growth, and resource exploitation in many low and middleincome countries can increase point and non-point PTEs sources, resulting in significant environmental degradation (Zhi et al. 2018). Therefore, understanding the efficient and accurate spatial distributions of PTEs in the soil is crucial for effective management and contamination remediation (Song et al. 2019).
The spatial prediction and distribution of PTEs in the soil using conventional techniques have been laborious and time-consuming. However, Digital soil mapping (DSM) in contemporary times is the utmost effective technique to map the spatial distribution of PTEs, types of soil, and soil properties by using field and laboratory observational methods coupled with spatial and non-spatial soil inference systems (McBratney et al. 2003;Minasny and McBratney 2016). Numerous geostatistical methods and machine learning algorithms have been developed to predict the spatial distribution of soil properties (including PTEs) within the DSM framework (McBratney et al. 2003). For example, ordinary kriging (OK), simple kriging, universal kriging, co-kriging (CK), regression kriging (RK), and empirical Bayesian kriging (EBK) are widely used geostatistical methods in DSM and environmental sciences (McBratney et al. 2003;Webster and Oliver 2007;Jiang et al. 2017;Gribov and Krivoruchko 2020). The OK is one of the most robust and popular geostatistical interpolation techniques used in DSM and the spatial distributions of PTEs in soil by considering spatial dependency and stationarity (Webster and Oliver 2007). The CK and RK are common geostatistical tools considering the relationship between primary and secondary variables (Mallik et al. 2020).
However, each geostatistical method makes many assumptions and limits or does not include various explanatory variables influencing the spatial prediction and distribution of PTEs in soils. For example, OK considers spatial autocorrelation and data stationarity (Webster and Oliver 2007) and excludes explanatory variables in the spatial predictions. The RK and CK try to overcome some OK limitations by detrending the data with linear regression and incorporating some covariates. Universal kriging (UK), one of the best linear unbiased predictors, handles the random variation by incorporating the trend model into the kriging method (Lark et al. 2006;Keskin and Grunwald 2018). However, the relationship between the response variables and the covariates in the RK (Shi et al. 2009;Sun et al. 2012), the three maximum number of covariates in CK (Wackernagel 1994;Giraldo and Herrera 2020), and the failure to find a spatial variance model in the UK (Lark et al. 2006), are further limitations to obtaining the acceptable predictions of soil PTEs. On the other hand, if any of the above premises or assumptions are not met, the outcome values by OK, RK, UK, and CK might be suboptimal, and predicting values at non-sampled locations lacks reasonable accuracy (Pilz and Spöck 2008).
EBK regression prediction (EBKRP) is an advanced geostatistical technique that integrates kriging with regression methods to make predictions more precise than either regression or kriging can achieve independently (Krivoruchko and Gribov 2019;Gribov and Krivoruchko 2020). The EBKRP method is also estimated regionally and accounts for regional effects (Gribov and Krivoruchko 2020). For example, the associations between the explanatory and dependent variables may change in different areas, and the EBKRP can precisely model these local changes. Nevertheless, despite the many advantages of EBKRP, it has several shortcomings. For example, the EBKRP does not identify the independent variables strongly associated with the response variables and the significance of the explanatory variables that influence the predicting variables. Therefore, random forest (RF) is a machine learning technique for the spatial prediction that solves the limitations associated with the EBKRP method.
The RF model is more promising than other machine learning approaches for accurately estimating and predicting soil PTEs because of its insensitivity to noise features, resistance to overfitting, and unbiased error rate measurement (Breiman 2001). RF can also be used to identify the influences and associations of explanatory variables. On the other hand, the method of RF algorithms did not clarify how the predictions were made, and some studies have labeled RF as a "black box" algorithm (Taghizadeh-mehrjardi et al. 2016;Hengl et al. 2018;Minasny et al. 2018). Furthermore, the RF algorithm only considers the relationships between the response and auxiliary variables, ignoring the impact of nearby observed data (spatial autocorrelation) (Guo et al. 2015). Hence, no single method has shown better performance for soil PTE prediction and distribution studies (Xiang et al. 2020); thus, combining different techniques can neutralize their weaknesses.
A combination of different techniques, hereafter called hybrid methods, refers to combining two conceptually different approaches to model the spatial variation of soil properties (Mirzaee et al. 2016). The hybrid methods were initially introduced by combining geostatistics with additional information, such as kriging with external drift and CK with auxiliary variables (Tziachris et al. 2019). Regression kriging was another hybrid method that gained popularity over time (Hengl et al. 2004). Recently, hybrid methods have been increasingly used to model geostatistical trends with machine learning models to enhance the prediction accuracy of interpolation methods (Matinfar et al. 2021). For example, it includes regression-kriging (RK) and artificial neural networks (RK-ANN) (Mirzaee et al. 2016), and RF with OK (RF-OK) (Tziachris et al. 2019).
Several studies have reported that hybrid models improve prediction accuracy (Dai et al. 2014;Mirzaee et al. 2016;Tziachris et al. 2019;Matinfar et al. 2021). However, few studies have been conducted on the performance of EBKRP over machine-learning approaches (Requia et al. 2019;Mallik et al. 2020). Specifically, no well-documented studies have reported the prediction performance of EBKRP, RF, and the hybrid methods of RF-EBKRP for the spatial prediction of PTEs in the soil. Furthermore, studies on the relationship with covariates using geostatistical techniques and machine learning models have not been conducted entirely (Hengl et al. 2018). It is also critical to identify the key auxiliary variables that directly correlate with their target PTEs to explain the spatial variability of soil PTEs in the prediction process. For example, soil fertility indicators for managing urban green space areas, such as applying fertilizer, irrigation, herbicides, and pesticides, and the functional types (industrial, commercial, medical, residential areas, etc.) influence the spatial predictions of PTEs in urban green space soils. Moreover, topographic factors for designing and improving greenspace landscapes, such as slopes, aspects, and elevations, play a role in predicting PTEs in greenspace soil. Therefore, the objective of this study was to identify the roles of auxiliary variables for predictive model performance, compare the prediction capabilities of the EBKRP and RF, and develop hybrid RF-EBKRP models for estimating the five soil PTEs (Pb, Cu, Zn, Cr, and Cd).

Study area descriptions
Shanghai is a densely populated metropolitan city in eastern China, located at 31.14° N and 121.29° E (Fig. 1). It is a coastal city that covers 6340.5 km 2 , of which 6218.65 km 2 are lands (Shi et al. 2008); the municipality covers 0.06% of China's total territory. The city has a subtropical monsoon climate with annual average rainfall and temperatures of 1122 mm and 15.8 °C. The main soil types were paddy and coastal saline.
Shanghai plays several vital roles in the country's main economic, financial, trade, and shipping centers. It is also one of China's largest industrial centers, with over 10,000 factories and industrial enterprises (Wang et al. 2009). The urban area is divided into three zones: the city center and inner and outer suburbs. The present study focused on the city center (Fig. 1). Many green space types and densities, including parkland, protective green space, square green space, subsidiary green space, etc., are mainly concentrated in the city center rather than in other parts of the city area (Shanghai Municipal Government (SMG) 2018). In addition, this area has experienced rapid industrialization and urbanization, resulting in the accumulation of soil PTEs, which is a significant concern.

Soil sampling and analysis
Two hundred surface soil samples (0-20 cm) were collected by random sampling strategies from the green space areas in 2018. The geographical locations of soil sampling points are shown as black dots in Fig. 1. Five sub-samples were collected around each sampling location using a soil corer (2.5 cm in diameter) and mixed thoroughly to obtain a representative composite sample at each location. Next, the composite samples were air-dried, and visible plant roots and residues were removed. Once the air-dried soils were screened, they were ground to pass through a 0.15 mm nylon mesh sieve to ensure the complete digestion of soil samples.
For each composite soil sample, 0.5 g of soil was digested with HNO 3 , HF, and HClO 4 as stated in the EPA 3052 method (EPA 1996). Then, the PTEs, such as Cu, Zn, Cd, Cr, and Pb, were measured using standard inductively coupled plasma mass spectrometry (ICP-MS, NexION 300X, USA). A sequence of soil quality assurance, quality control, and geochemical reference materials, supplied by the National Research Center for Certified Reference Materials of China, was also checked. The limit of detection (LOD) and limit of quantification (LOQ) for each PTEs and verification are described in a previously published work . The LODs of Cr, Cu, Zn, Cd, and Pb were 0.47 mg/kg, 0.25 mg/kg, 0.70 mg/kg, 0.01 mg/kg, and 0.30 mg/kg, respectively. The LOQ of the five PTE mentioned above was four times their respective LOD.

Covariates selections and preparations
Explanatory variables play an essential role in predicting PTEs in the soil by indirectly affecting the distribution of elements (Maas et al. 2010;Kheir et al. 2014). The covariates used in this study could be categorized into soil fertility indicators, topographic features, and urban functional types. activities, such as fertilizers, irrigation, and herbicides, for managing greening areas ). Soil fertility indicators such as electrical conductivity (EC), pH, total nitrogen (TN), total phosphorus (TP), total potassium (TK), and soil organic carbon (SOC) were analyzed using standard measurement methods from 200 collected soil samples. For example, pH was determined using a pH meter The detailed descriptions and preparations of each explanatory or auxiliary variable are explained below:

Soil fertility indicators
The soil fertility indicators were selected as auxiliary variables owing to the influence of PTEs by anthropogenic The DEM data was used to generate slope and aspect data ( Fig. 2B) by the ArcGIS Pro spatial analyst tools.

Urban functional types
Studies demonstrated that different urban functional areas, such as industrial, commercial, medical, and residential, were used as covariates due to their better reflection of the spatial variations of PTEs in soil (Shi et al. 2021;Wang et al. 2021). This study used Sentinel-2 A multispectral images from Copernicus to combine and classify urban functional types. Sentinel-2 A 10 m resolution images were downloaded from https://scihub.copernicus.eu/. Then, the supervised machine learning classification method was used by ArcGIS Pro and followed the basic classification procedures such as multiple image preprocessing steps, training sample selection, training, classifying, and assessing accuracy (Jim et al. 2003). Finally, five urban functional types (residential areas, industrial areas, green space areas, public service areas, and water bodies), with an overall accuracy (Kappa index value) of 0.875, were used as covariates for the spatial predictions of soil PTEs (Fig. 2 C).
(soil: water = 1:5), and EC was measured using a conductivity meter standardized with a salt solution (Smith and Doran 1996). TN and TP were calculated by the Kjeldahl (Bremner 1960) and Olsen methods (Olsen et al. 1954). A flame atomic absorption spectrophotometer and loss on ignition at 550 °C were used to determine the TK and SOC contents (Bremner and Jenkinson 1960). Finally, due to the input covariates in EBKRP tools in the raster format, the soil fertility indicator data was transformed into a raster by the EBK interpolation method (Fig. 2 A). We used the EBK interpolation techniques because EBK considers the error of the semivariogram model estimation, unlike the other kriging techniques (Krivoruchko and Gribov 2019; Gribov and Krivoruchko 2020).

Topographic features
Among numerous terrain or topography parameters, digital elevation models (DEMs), slopes, and aspects are the most influential factors for PTEs in the soil (Qiao et al. 2017;Ballabio et al. 2018). A 30 m resolution DEM was downloaded from the Resources and Environmental Scientific Data Center, Chinese Academy of Sciences (http://www.resdc.cn).

Random Forest (RF) for spatial predictions and variable importance measures
RF is a classification and regression method based on the aggregation of many decision trees, first invented by (Breiman 2001) and recently available in the literature (Prasad et al. 2006;Biau and Scornet 2016). In addition, several studies have proven that it is one of the best machine learning techniques currently available, and its detailed mathematical formulation has been reported elsewhere (Breiman 2001;Prasad et al. 2006;Boulesteix et al. 2012;Vaysse and Lagacherie 2015;Olson et al. 2017;Nussbaum et al. 2018).
This study used the RF method with the "randomForest" package in the R software (R CoreTeam 2021) to address two main tasks: (1) to create a prediction and distribution map of PTEs using supervised machine learning. (2) Assess and rank explanatory variables according to their ability to predict response variables (PTEs). The latter are called variable importance measures, which are directly computed for each predictor within the RF algorithm. As a result, only one or several other predictors can accurately predict the response variables (PTEs) (Fig. 2), and their rank levels were expressed as the relative importance in percentage (Fig. 4).
RF is used for prediction by aggregating many decision trees. First, each decision tree is built using randomly

Spatial modeling and prediction
Three approaches were used to predict the spatial distribution of soil PTE content: the EBKRP, RF, and hybrid method (RF-EBKRP). Below are detailed descriptions of each technique.

EBKRP interpolation method
The EBKRP is an advanced geostatistical interpolation method that combines EBK with regression analysis (Krivoruchko and Gribov 2019; Gribov and Krivoruchko 2020) provide a detailed mathematical description of the EBKRP model. The EBKRP models automatically solve the most challenging features of the kriging model. Other kriging methods require manual adjustment of parameters to obtain accurate results; however, EBKRP automatically adjusts the parameters needed through the subsetting and simulation processes (Krivoruchko and Gribov 2019).
The EBKRP in this study followed the following procedures: (1) The raster's form of covariates was changed into their principal components to reduce multicollinearity problems (explanatory variables associated with each other) (Gribov and Krivoruchko 2020).
(2) The converted covariates in the form of principal components were used in the regression model (Fig. 2). Once the input covariates were organized and prepared based on the model requirements, the final step was to select an appropriate semivariogram model and perform a log-data transformation (if required) of the predictive variables (soil PTEs) (Table S2). A good testing data. Before judging how well the model worked, the Moran's I Index and Semivariogram (Sill to nugget ratio) were used to test the assumption that the training and validation data sets were independent of each other. As a result, the features of the validation data were randomly distributed and had moderate to weak spatial dependency (Table S1). The model performance and good estimators were determined using different performance metrics (PMs) or accuracy statistics using the testing data sets. The model PMs are frequently measured using the root mean squared error (RMSE), mean absolute percentage error (MAPE), and R 2 . However, each PM has a shortcoming in evaluating the model efficiency. For example, the MAPE and R 2 are more sensitive to extreme values (Chai et al. 2014;Bhagat et al. 2019), and the RMSE is more prone to outliers (Chai et al. 2014). Consequently, multiple PMs are used to bridge the gaps in each metric (Tofallis 2015;Morley et al. 2016). Table 1 presents the mathematical definitions and descriptions of each metric.

Statistical software used
The spatial distribution map and prediction of soil PTEs by EBKRP were created using ArcGIS Pro (Ver. 2.7, Esri Inc., https://www.esri.com/en-us/arcgis/products/arcgis-pro/ overview). Descriptive statistics, a Spearman's rank correlation, and prediction of soil PTEs by RF were performed using open-access R Statistical Software (Ver. 4.1.1), (R CoreTeam 2021). The significance P-values for the Spearman correlation rank analysis were 0.05 and 0.01. Table 2 shows the summary statistics and mean values for the five PTEs and other covariates in the urban greenspace soils. Pb, Cu, Zn, Cr, and Cd concentration ranges were 17-175.75, 15-225.44, 59-798.87, 50.50-104.02, and 0.06-3.68 mg/kg, respectively. Table 2 also shows the PTE values compared to the soil background values in Shanghai generated portions of the original (training) data. For validation purposes, 25% of the data were excluded from training data sets. Once the data is categorized as testing and training data, the second step is to train the model by adjusting the forest parameters, such as the number of trees (Ntree), the number of variables tried (Mtry), and the minimum leaf size (Nodesize). The default value for Ntree is 100. Increasing Ntree will result in a more accurate prediction, but it will take a long time to calculate. Mtry refers to the number of explanatory variables selected randomly at each split tree and is usually determined by the square root of the total variables. Finally, Nodesize is the minimum amount of training data required to continue the tree growth process. After many trials and errors, Ntree, Mtry, and Nodesize are defined as 500, 5, and 5, respectively. Once the best-trained models were selected based on out-of-bag (OOB) errors, variable importance, root mean square error and R 2 , predictions of unknown locations of soil PTEs were predicted by RF in the form of raster features.

Hybrid method
Hybrid methods have been developed to examine the spatial variations in residuals to improve the predictions of un-sampled locations (Matinfar et al. 2021). The hybrid RF-EBKRP method used in this study combines a non-spatial approach (RF) and spatial interpolation techniques (EBKRP) to improve the accuracy and reliability of soil PTEs prediction. The implementation procedure followed three main stages, and the overall methodological approaches are shown in Fig. 3. First, the trained RF model was obtained, and the prediction was performed as explained in Sect. 2.4.2. Subsequently, a residual by RF was generated as follows (Eq. 1): Where r RF (x i ) is the residual generated by the RF model at location x i , Z (x i ) is the measured PTE value, and y RF (x i ) is the value predicted by RF. In the second stage, the residuals or error terms generated by RF (r RF (x i )) were imported to ArcGIS Pro 2.7 software, and the prediction was carried out using the EBKRP method. Finally, the final estimated soil PTEs contents ŷ(x i ) by hybrid RF-EBKRP methods was obtained in Eq. 2.
ŷ (x i ) = y RF (x i ) + r EBKRP (xi) (2) Where y RF (x i ) is the PTE value estimated by RF at location x i, and r EBKRP (xi) is the RF residual value calculated by the EBKRP method.

Model validation and performance
For validation purposes, the data for each prediction model was randomly divided into 75% training and 25% x i and y i denote the measured and predicted value at location i, respectively, and n stands for the number of observations MAPE µ a is the mean value of the x values elevation features were significantly correlated with different soil PTEs (Table 3).

Key factor for PTEs predictions improvement
RF is used to identify the key factors that influence the spatial distribution of PTEs in soil and can help improve predictions. The relative importance (in percentage) for each covariate of soil PTEs is shown (Fig. 4). Except for Cr, soil fertility indicators, such as OM and SOC concentrations, were the most effective predictors for explaining the spatial variation of all PTEs. Total N and K were ranked first and second in their influence on soil Cr. Other soil fertility indicators, such as EC, total P, C/N, and pH, were classified as having moderate to weak influences on the spatial distribution and prediction of PTEs in green space soil (Fig. 4).
The topographical covariates (elevation, slope, and aspect) also contributed to improving the PTEs predictions in urban green space soils. Slopes had a minor impact on Cr distributions and dynamics in greenspace soils. The urban functional types, including residential, industrial, public services, etc., better explained the spatial variations and distributions of soil Pb, Zn, and Cd, and lower influence soil Cu and Cr.

Spatial prediction and distribution of soil PTEs
The spatial distribution maps of the five soil PTEs determined by the three methods are shown in Fig. 5. The maps generated by EBKRP showed a relatively smooth and noisy surface, whereas the maps produced using RF and RF-EBKRP showed distinct geographical distributions. The discrete geographical patterns of PTEs in RF-EBKRP and (Wang and Luo 1992). As a result, soil samples' Pb, Cu, Zn, and Cd concentrations exceeded background values by 41.34%, 27.67%, 43.61%, and 48%, respectively. However, the sampled Cr value was 5.75% lower than the reference value (Table 2). Similarly, the coefficients of variation (CV) showed that Cr (13.61%) had lower values and Cd had higher values (116%). The soil fertility status of urban greenspaces was assessed using TN, TP, TK, OM, C/N, and SOC, which averaged 1. 39, 0.84, 18.31, 43, 26.44, and 24.92 g/kg, respectively (Table 2). Soil alkalinity and salinity levels were comparatively higher, with EC and pH values averaging 0.26 mS/cm and 8.11, respectively. Except for Cr, all PTEs had higher values for kurtosis and skewness, which showed that the data were not normally distributed. PTEs with high SD values also showed non-normal distributions. So, the PTEs that didn't follow a normal pattern were turned into log transformations before the modeling and prediction analysis (Table  S2).

Spearman's rank correlation analysis
Due to the non-normal distribution of the original data, Spearman's correlation rank analysis was used to determine the relationship between PTEs and other explanatory variables ( Table 3). Statistically significant correlations (p ≤ 0.01) were observed between the PTEs. Furthermore, a statistically significant correlation was found between PTEs and most soil fertility covariates. For example, SOC, TN, TP, and OM were positively correlated with all the PTEs. In contrast, soil pH was negatively associated with other PTEs except for Pb and Cr. Except for Cr and Cd, the topographic  Fig. 6). Overall, the three interpolation methods predicted lower Pb, Cu, Zn, and Cd concentrations in the northwestern parts of the area and higher Pb, Cu, and Zn concentrations on the eastern sides of the study area. Similarly, all models generated high and low Cr concentrations in the southern and central regions, respectively.

Model performance in predicting soil PTEs
The difference between observed soil PTEs and prediction data in testing data sets (25% of sampled data) was used to evaluate the prediction accuracy and validation of the various models (Table 4; Fig. 6). Three metrics were independently calculated to assess the accuracy of each model: RMSE, MAPE, and R 2 . Generally, the hybrid RF-EBKRP model was the best prediction method, followed by the RF model. The RF method decreased RMSE by 0.26-13.1 mg/kg and MAPE by 0.94-8.02% on average compared to the EBKRP interpolation method ( Table 4). The RF-EBKRP model reduced RMSE by 0.47-77.63 mg/kg and 0.21-63.92 mg/kg to predict soil PTEs by EBKRP and RF. Similarly, the RF-EBKRP method decreased MAPE by 2.77-17.06% and 1.21-12.43% on average compared to the EBKRP and RF methods.

Important factors improving soil PTEs predictions
The accumulation and distribution of PTEs in urban green space soil are influenced by various factors, including agrochemicals, wastewater irrigation, traffic and industrial sources, and complex adsorption mechanisms by soils (Weng et al. 2002;Luo et al. 2012;Zhang et al. 2021b). However, it is not straightforward to determine the exact relationship between soil PTEs and their potential influencing factors. RF is a machine learning algorithm used to discover the relative importance of possible factors affecting the dynamics of PTEs in the soil (Breiman 2001). Urban functional types, including residential, industrial, public services, etc., are identified as the main contributors to Pb, Zn, and Cd in the Table 3 Spearman's correlation coefficient results within and between PTEs and covariates Similarly, other studies on Cd indicated that soil type significantly affected the distribution and content of Cd in the soil (Cao et al. 2017). Other agrochemicals, such as pesticides and herbicides, to manage trees and grasses in the gardens also contributed to the availability of PTEs in urban green space soils, such as Cd, Zn, and Cu (Mico et al. 2006;Zhang et al. 2021b). SOC can also influence the distribution of PTEs by retaining the formation of the organic-metal complex (Shi et al. 2013;Hong et al. 2019). Topographic covariates also play an essential role in predicting soil PTEs by affecting the distribution of PTEs in the soil (Kheir et al. 2014). For example, elevation and slope are the most influential topographic factors for PTEs distribution and content in the soil (Qiao et al. 2017;Ballabio et al. 2018). The results derived by RF indicated that elevation, slope, and aspect were the essential factors contributing to the prediction accuracy improvement in Pb, Cu, and Zn (Fig. 4), and lower prediction errors in the models. They largely determine the concentration, mobility, and availability of PTEs in the soil by affecting runoff, drainage, and soil erosion (Dubovik and Dubovik 2016;Liu et al. 2020). In addition, it may influence the migration of PTEs from atmospheric deposition into the soil by affecting their distribution characteristics (Qiao et al. 2017;Dubovik and soils (Fig. 4). Some studies found that using urban function type as a covariate reflects the diversity of human activities in urban environments and improves digital soil mapping by prediction models. For example, studies on digital mapping of Zn in urban topsoil using urban functional types as predictor variables better explained the spatial variations of Zn in the soil by RF models (Shi et al. 2021). Other studies on heavy metal pollution in urban river sediment also found that heavy metal occurrence strongly correlated among the urban functional areas (Wang et al. 2021).
The soil fertility covariates identified by the RF and Spearman correlation coefficients confirmed that soil fertility covariates were strongly associated with soil PTEs. The higher variations in CV for the studied PTEs also imply that the elements were mainly influenced by anthropogenic activities, such as fertilizers, irrigation, and herbicides, for managing greening areas ). Moreover, the considerably higher Pb, Cu, Zn, and Cd values in greenspace soils than those of the reference values suggested that a considerable number of PTEs were sourced from practical practices to improve the poor conditions of greenspace soils. For example, studies (Huang et al. 2007) showed that Cd accumulation in agricultural soils is associated with agrochemicals and organic manure for soil fertility Compared with EBKRP, the RF and RF-EBKRP hybrid methods achieve good predictions and interpretability of soil PTEs content mapping.
The better prediction and mapping of soil PTEs by RF could be because the RF model is less dependent on sample size and can model nonlinear relationships in the dataset (Khaledian and Miller 2020). Furthermore, the RF method overlooks geostatistical methods' stationarity and variogram assumptions (Hengl et al. 2018). RF is also used to provide variable importance measures for each predictor and simplifies the model interpretation (Behrens et al. 2010;Ließ et al. 2012;Khaledian and Miller 2020). Numerous comparisons of RF with other geostatistical methods have also demonstrated the ability of the RF model to predict the dynamics of soil properties (Vaysse and Lagacherie 2015;Hengl et al. 2018;Huang et al. 2019;Zhang et al. 2021b). For instance, Hengl et al. (2018) found that RF can result in accurate and unbiased predictions in different versions of Dubovik 2016) discovered that aspect is a significant factor controlling the distribution of the majority of heavy metal bulk, mobile, and acid-soluble compounds.

Evaluating models' predictive capabilities
PTEs contamination in green space areas' soil is a major environmental and health issue because people can quickly contact the soil during socializing and gathering after work in green space areas (Manta et al. 2002). Therefore, understanding the efficient and accurate spatial distribution information about PTEs in the green space soil is crucial to effective management and contamination remediation. However, because of multiple sources and spatial variability of PTEs, getting an accurate prediction and mapping of soil PTEs is difficult by interpolation techniques (Ha et al. 2014). Thus, a better methodology should be developed to improve the accuracy of PTEs content mapping in the soils.

Conclusions
Quantifying the effect of auxiliary information and evaluating the efficiency of interpolation methods are crucial for increasing the accuracy of soil PTE estimation. Soil fertility characteristics (SOC, OM, TN, TP, and TK) and urban functional types (residential, industrial, public services, etc.) are important covariates that enhance soil PTE distribution and prediction maps by RF and EBKRP methods. Similarly, topographic indicators, including elevation, slope, and aspect, were ranked as crucial factors contributing to improving Pb, Cu, and Zn prediction accuracy in green space soil. Furthermore, machine learning methods (RF) enhance the precision of soil PTE prediction more than advanced geostatistical techniques (EBKRP). The improvement was 21-50% for RMSE, 36-92% for R 2 , and 8-32% for MAPE. However, hybrid methods (RF-EBKRP) increased accuracy by comparing the individually predicted models with 122% and 15% in EBKRP and RF, based on R 2 , respectively. In conclusion, in addition to incorporating covariates into the models, combining kriging residuals with the machine learning method (RF-EBKRP) resulted in a promising approach for improving the distribution accuracy and mapping of PTEs in the soil.

Acknowledgements
The authors express their gratitude to the organizations supporting this work. We also thank the CAS_TWAS presidential fellowship international doctoral program. kriging. Moreover, studies on predicting soil organic matter also reported that the machine learning method (RF) improved the prediction accuracy by (250% increase in R 2 ) over the OK geostatistical interpolation techniques (Tziachris et al. 2019).
However, the RF algorithm did not show functional relationships between the target and predictor variables nor explain how the predictions were made. Because of these drawbacks, some studies referred to RF as a "black-box" algorithm (Taghizadeh-mehrjardi et al. 2016;Hengl et al. 2018;Minasny et al. 2018). However, by combining the residual results of EBKRP with the predictive values by the RF method, the limitations associated with individual models can be overcome and result in a better PTEs spatial distribution map. Therefore, this study demonstrated the RF-EBKRP hybrid method outperformed individual models for predicting soil PTEs. Several studies in digital soil mapping have confirmed that combining geostatistical methods with machine learning methods results in better prediction performance than individual models (Dai et al. 2014;Guo et al. 2015;Mirzaee et al. 2016;Song et al. 2017;Tziachris et al. 2019;Matinfar et al. 2021). For example, Matinfar et al. (2021) reported that combining RF and OK provided a more accurate prediction of SOC than the separated methods. Similarly, other studies found that a hybrid artificial neural network (ANN-OK) predicted SOC better than the individual ANN and OK prediction methods (Mallik et al. 2020). The same SOC studies, however, revealed that the advanced geostatistical-based EBKRP method outperformed the hybrid ANN-OK methods (Mallik et al. 2020).
Generally, unlike the individual models, the proposed hybrid RF-EBKRP methods have the benefits of linking the predicting variables to various covariates and complex nonlinear relationships and producing better accurate soil PTEs distribution maps. Another advantage of these methods is that they consider more significant spatial and non-spatial relationships between multiple variables than individual machine learning and geostatistical models (Matinfar et al. 2021). Moreover, the implementation procedures of the combined RF-EBKRP methods were computationally straightforward, extending only one more step of the individual models. The preprocessing stages, such as training predictive models, descriptive analysis of original data, and