## 2 Methods  
### Methods for Map Construction Based on Spatial Databases  
Pollutant concentrations can vary greatly across the surface of the earth, so it is very difficult to obtain an accurate spatial distribution of heavy metals. The presence of a certain proportion of samples that exceed the established regulatory threshold has been a classic method of characterizing the degree of soil contamination [47; 48]. However, the classical estimation method has many limitations whichprobably cause errors or uncertainties in pollution assessments. The classical statistical method usually requires that data be subject to a number of assumptions: independence of observations, exact or approximate normality of observations, and large and repeated sampling. Mapping the spatial distribution of pollutants requires the use of spatial interpolation methods. The most commonly used interpolation methods are inverse distance weighting (IDW), kriging, and spline. Interpolation accuracy is related to the precise definition of the polluted area and its boundaries. Consequently, this directly affects the accuracy of pollution assessment. There have been many studies of the performance of the spatial interpolation methods mentioned above, but the results are not clear-cut [49]. Some of them found that the kriging method performed better than IDW [50], while others showed that kriging is no better than alternative methods [51]. The first considered interpolation method of kriging is based on the assumption that the interpolated parameter can be considered as a regionalized variable. As with IDW, the kriging estimation is given by a linear combination of the observed values with weights. Depending on the stochastic properties of random fields, different types of kriging are applied. The type of kriging defines the linear constraint on weights implied by the unbiased condition. There are several types of kriging, including simple kriging, ordinary kriging, universal kriging, etc., and ordinary kriging is the most commonly used method. The ordinary kriging weights are derived from the kriging equations using the semivariance function. The parameters of the semivariance function and the nugget effect can be estimated using the empirical semivariance function [52]. The unbiased assessment of the semivariance function is half the average squared difference between paired data values:  
\[\gamma(p)=\frac{1}{2N(h)}\sum\limits_{i=1}^{N(h)}\left[z(x_{i})-z(x_{i}+h) \right]^{2} \tag{1}\]  
where \(\gamma(h)\) is the semivariance value at distance interval h,N(h) is the number of sample pairs within the distance interval h, and \(z(x_{i}+h)\) and \(z(x_{i})\) are sample values at two points separated by the distance interval h.  
The next kind of interpolatorsareradial basis functions (RBF). This is the name given to a large family of exact interpolators that use a basic equation that depends on the distance between the interpolated point and the sampling points. The RBF prediction value can be expressed as the sum of two components [53]:  
\[Z(x)=\sum\limits_{i=1}^{m}a_{i}f_{i}(x)+\sum\limits_{j=1}^{n}b_{j}\psi(d_{j}) \tag{2}\]  
where (d\(j\)) shows the radial basis functions and d\(j\) is the distance from sample site to prediction point x. f\({}_{j}\)(x) is a trend function, a member of a basis for the space of polynomials of degree \(<m\). The coefficients a\({}_{i}\) and b\({}_{j}\) are calculated by means of the resolution of the following system of n + m linear equations; n is the total number of known points used in the interpolation. Polynomial interpolation is a process of finding a formula (often a polynomial) whose graph will pass through a given set of points. Global polynomial interpolation corresponds to the polynomial of the entire surface, while local polynomial interpolation can be seen as a combination of global polynomial methods and a moving average procedure. Instead of fitting the polynomial to the entire dataset, it is fitted to a local subset defined by the window, as in the moving average model.