Digital mapping and spatial characteristics analyses of heavy metal content in reclaimed soil of industrial and mining abandoned land

The reclaimed soil properties of industrial and mining wasteland have strong spatial specificity. The paper aimed to screen out a hybrid multifractal and kriging (Named as Mkriging) method for digital mapping and scientifically reveal the spatial distribution characteristics in view of heavy metal in reclaimed soil of industrial and mining abandoned land. The results of the study showed that for reasons of history and reclamation, the original samples of heavy metals in reclaimed soil of industrial and mining abandoned land showed a very large range and variation degree, the C0/(C0 + C1) values of different heavy metals basically were all greater than 50%, random factors played a dominant role. The five kinds of heavy metals in reclaimed soil were in the following descending order in terms of homogeneity: Cd, As, Hg, Ni and Cr. Compared with ordinary Kriging method, the relative improvement of root mean squared errors of elements Cd, Hg, As, Cr and Ni based on Mkriging were 95.28%, 61.74%, 78.54%, 82.51% and 83.58% respectively. The higher the fractal degree of heavy metals in reclaimed soil was, the higher the prediction accuracy will be. Mkriging method is more suitable for spatial prediction of heavy metals in reclaimed soil of industrial and mining abandoned land.

Industrial and mining abandoned land is a special type of land space damaged due to industrial and mining activities. Its reclamation and utilization play an important role in ameliorating ecological environment, optimizing the development and layout of land space and promoting resource conservation and ecological civilization. The digital mapping and spatial characteristics analyses of soil properties in the reclaimed area are an important link of reclamation monitoring.
Due to the long indigenous smelting sulfur, the industrial and mining abandoned land in the southwestern region of China were seriously polluted by heavy metals, which is one of the obstacles holding back the progress of ecological civilization, leads to the reduction of yield and quality of the crops on the reclaimed land and threatens the safety of the ecosystem and humans in the reclaimed area 1 . A lot of research has been done on heavy metals in soil at home and abroad  . Currently, the research focuses on the non-restructured soil in cities, vegetable farms, orchards, and farmland near a mining area 6,7,11,20 , and analyzes and evaluates the pollution condition in soil and their impact on safe land use based on soil environment standard or regional geological background from the perspective of sampling points  . There is a complicated historical background in the wasteland of industrial and mining. There are significant differences in the reclamation project. The heavy metal content in the rebuilt soil is affected by the terrain and human factors. Even in a very small distance, there may be a great change. The variability and variability are strong. A conventional and regular Kriging method requires minimum error variance and has a very strong smoothing effect, which will reduce or eliminate the areas with strong spatial variability and makes the whole study area tend to a median (or interval). It is impossible to realize the conversion of heavy metals from the point to surface in the reclamation soil, to influence the follow-up monitoring, evaluation and selection of management and protection measures. Therefore, this article aims to find a digital mapping method that can not only highlight specific values and effectively maintain and reproduce the spatial structure of original data-based variables but also guarantee the overall accuracy and model fitting effect of the prediction result. Mkriging method is in fact a kind of extended interpolation of moving weight average, which is multiplied by a singularity index as correction factor since interpolation of original moving weight average. It overcomes the impact of smoothing effect, and properly maintains the spatial variation degree. At present, the multifractal theory has been applied in mineral research, agricultural research and other nonlinear spatial analysis [25][26][27][28][29][30][31] , but its application in digital mapping and spatial distribution characteristics analyses of heavy metals in soil, particularly in reclaimed soil of industrial and mining abandoned land, is seldom reported. This article chose heavy metals in reclaimed soil of industrial and mining abandoned land in the southwestern region of China as a study object and carried out the following research: (1) Using Mkriging method, the spatial mapping of heavy metals in reclaimed soils from industrial and mining wastelands was realized on the basis of variogram analysis and multifractal analysis; (2) Reveal the spatial distribution characteristics of heavy metal content in reclaimed soil of industrial and mining abandoned land. The above research will provide methodological guidance for spatial prediction of heavy metal content in reclaimed soil from industrial wasteland.

Overview of the Study Area and Data Processing
Overview of the study area. The study area lies in the southwestern region of China, and began mining and beneficiation of sulfur ore in the late 1950s. In 2004, the mine was bankrupt, restructured and handed over to local township government for management. After 40 years' mining and smelting, solid waste (such as sulfur slag) accumulation is very large. The digging, occupation and pollution of the mining area are serious, and the ecological environment of the whole mining area is harsh. The total reclaimed area is 266.49 hm 2 . The work of reclamation was initiated in 2013 and completed with acceptance in 2014. The study area, with the altitude being 500~1100 m, is in the transition belt from the southern margin of Sichuan Basin to Guizhou Plateau, and within Chishui River basin. The study area bears the climatic features of Sichuan Basin and Guizhou Plateau. The type of soil is yellow soil.
Data acquisition and processing. Soil samples were taken in July 2016. A total of 58 samples were selected to represent a range of reclamation measures and other factors (Fig. 1). The interval between all the samples is less than 500 m, which is within the range of the semi-variogram of each heavy metal (Fig. 2). The sampling depth was 0~20 cm. The locations of the sampling sites were recorded using a global positioning system (GPS) receiver (Garmin3790T). Three to five soil samples were collected at depth 0-25 cm over a circle of radius 10 m surrounding the specified sampling location and mixed thoroughly. About 1.5 kg soil was collected from each mixed sample. According to the relevant standards such as Soil Quality-Determination of Lead, Cadmium-Graphite Furnace Atomic Absorption Spectrophotometry (GB/T17141-1997), Soil and Sediment-Determination of Mercury, Arsenic, Selenium, Bismuth, Antimony-Microwave Dissolution/Atomic Fluorescence Spectrometry (HJ 680-2013), and Soil and Sediment-Determination of Mercury, Arsenic, Selenium, Bismuth, Antimony-Microwave Dissolution/Atomic Fluorescence Spectrometry (HJ 680-2013) and so on, five types of heavy metals in the reclaimed soil, including cadmium (Cd, mg/kg), arsenic (As, mg/kg), chromium (Cr, mg/kg), mercury (Hg, mg/kg) and nickel (Ni, mg/kg) were tested.

Methods of digital mapping and spatial distribution characteristics analyses. Geostatistics.
Geostatistics is a science based on the theory of regionalized variable and using variogram function as a main tool to study natural phenomena. The ordinary Kriging method uses spatial autocorrelation to estimate the value of the estimated point in the minimum variance. where, z(x 0 , y 0 ) is the value of the estimated point, n is the number of adjacent related points, and λ i is the kriging weights. γ(x i − x j ) is the value of variation function between x i and x j and μ is the Lagrange multiplier. By solving the above linear equations, all the weights λ i and Lagrange multipliers μ can be obtained. This article chose ordinary Kriging as a control method. For more information, please refer to related literature [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47] .
Multifractal analysis. Multifractal analysis method firstly obtains the probability distribution of unevenly distributed spatial variables (Heavy metal of reclaimed soil in the article) by box-counting method through a computer.
As for any given real number q, a partition function χ q (l) may be defined, and probability μ i (l) is weighted and summated through the q th power.
Its mathematical expression is:  Fig.  2 were model results that provide the best fit for the data points. Binned values are shown as circles, triangles, meters, and crosses, and were generated by grouping (binning) empirical semivariogram points together using square cells that are one lag wide. Nugget(C 0 ) represents spatial variability due to random errors such as measurement errors and microscale processes; Partial sill(C 1 ) is the structural variance, which indicates the spatial structure caused by spatial correlation. where l is measurement scale (side length of the box); N(l) is the number of boxes; μ i (l) is probability distribution of spatial variable. If spatial variable bears a fractal feature, then χ q (l) ≈ l τ(q) , i.e., partition function χ q (l) and l have a relation of power function. τ(q) generally refers to a quality index. The singular exponent α and multifractal spectrum f(α) can be further calculated.
The τ(q) and q are transformed into f (α) and α through the Eqs 4 and 5. Multifractal spectrum width Δα and Rayleigh dimension D q can be further calculated through the following formula [48][49][50][51] .
where α max and α min is the maximum and minimum value of the singular exponent, respectively.
Mkriging. Considering the spatial distribution characteristics of the properties of reclaimed soil, we combined fractal theory with ordinary Kriging (OK) and conducted spatial prediction of quality elements of reclaimed land by Multifractal Kriging (Named as Mkriging) method to eliminate the impact of the smoothing effect on prediction result. Mkriging is an extended moving weighted average interpolation method. It uses the result of moving weighted average times a factor relevant with measurement scale and singularity exponent as an estimated value of the regionalized variable. Usually, the mean aggregation of spatial variable varies with measurement scale 28,[52][53][54][55][56] .
The relation between the measure and the dimension of a small square with a length of l is like a formula (8).
Taking the square with a length of l as the center, the measure in the square of Nl is: N Simultaneous Equations (8) and (9) can be obtained: If the mean value in the square of Nl is estimated by formula (1), the multidimensional fractal Kriging estimation formula in two-dimensional case can be obtained. Analysis on the effect of spatial digital mapping. We evaluated the mapping effect from three aspects: reproducibility, prediction accuracy and model simulation effect. We used the coverage diagram and coverage ratio of specific value (CRSV) to reflect the reproducibility. We formed a maximum (minimum) specific value group consisting of 20% of the maximum (minimum) values of measured values and predicted values. The CRSV value is the coincidence percentage of coverage of minimum (maximum) specific value of measured values and predicted values. The higher the percentage is, the better the effect of spatial prediction will be 54 . A cross validation method was adopted to validate the simulation accuracy and model fitting effect. Root mean squared errors (RMSE) and mean squared deviation ratio (MSDR) were adopted to validate the spatial prediction accuracy and model simulation effect. They are defined as: value is, the more accurate the prediction result will be. MSDR was used to evaluate the degree of fitting of theoretical variation function. The closer the MSDR value is to 1, the more accurate the fit variation function will be 45,47 .

Results and Analyses
Variogram analysis. The typical statistical characteristic values, optimum variogram functions and their relevant parameters of five heavy metals Cd, Hg, As, Cr and Ni in real samples were obtained based on the SPSS20.0 and ArcGIS10.2.
The full-sample amplitudes of heavy metal elements Cd, Hg, As, Cr and Ni in reclaimed soil were 6.01 mg/kg, 0.56 mg/kg, 28.74 mg/kg, 461.27 mg/kg and 145.70 mg/kg. The coefficients of variation (CV) of the above five kinds of heavy metals were 97.34%, 51.92%, 46.66%, 50.54% and 37.91% in turn, and the CV value of Cd was above 90%. There is mutability and strong spatial variability of heavy metal content of reclaimed soil of industrial and mining abandoned land from the variation and CV. This was mainly due to the difference of types and degrees of damage caused by the mining activities, and the difference of the reclaim measures and standards adopted in the process of reclamation.
Referring to the Technical Specification for Soil Environmental Monitoring, the variability coefficient and relative deviation are used to determine whether the sample points meet the research requirements. The formula is as follows.
where, n is the number of samples; t is the value of t at a certain degree of freedom at a given confidence level (generally 95% for soil monitoring); CV is the coefficient of variation (%) obtained from the collected data; m is the acceptable relative deviation (%); and soil monitoring is generally limited to 10-20%. According to the analysis results, the average coefficient of variation of all monitoring indicators was 56.87%, and the acceptable relative deviation m was taken as the minimum value of 10%. Through calculation, we need to set up sampling points 49. Based on the above results, combined with the variogram range and the area of the study area, the 58 sampling points meet the requirements. Table 1 indicated that different types of heavy metals have different statistical distributions. The distribution characteristics of heavy metals in reclaimed soil decided that it was not advisable to adopt ordinary Kriging with strong smoothing and central effects (refers the situation that the smoothing action of the model makes values of simulation converge towards a narrow range), which is unable to meticulously depict the mutation law of some zones.
Nugget/sill ratio (C 0 /(C 0 + C 1 )) refers to the ratio of spatial heterogeneity caused by random parts to total variation of the system and can usually be used to measure spatial correlation of variables. The smaller the ratio is, the stronger the spatial correlation will be. If the ratio is smaller than 25%, it means the variable has strong spatial correlation; if the ratio is 25~75%, it means moderate spatial correlation; if the ratio is greater than 75%, it means weak spatial correlation 43,44 . The nugget/sill ratios of heavy metals Cd, Hg, As, Cr and Ni in reclaimed soil were 50.32%, 47.15%, 63.88%, 66.29% and 65.06%, respectively, close to 75%, showing moderate spatial autocorrelation (Fig. 2). C 0 , C 1 and C 0 /(C 0 + C 1 ) indicate that except Hg, spatial variability (C 0 ) values brought by measuring error, microscale process and other random parts were all greater than structure variance, i.e., C 0 /(C 0 + C 1 ) values were all greater than 50%, and random factors played a dominant role. Spatial variability of heavy metals mainly stems from soil covering, improvement of soil fertility, soil pH conditioning measures and random factors.
Fractal Analysis. The multifractal analysis was used to calculate fractal parameters (Table 2), and to draw corresponding Rayleigh spectrum and multifractal spectrum of different heavy metals in reclaimed soil (Fig. 3). The above figures and related parameters may quantitatively and intuitively reveal the local characteristics of multifractal such as spatial heterogeneity, heterogeneity and shape characteristics of heavy metals in reclaimed soil, providing a basis for determining whether spatial prediction may be conducted by Mkriging or not.
The process dividing fractal into hierarchical zones by setting q at different values is called multifractal analysis. In theory, the larger the value range of q was, the better the result of multifractal analysis will be. In actual computation, however, with the increase of value q, workload will be multiplied and when q is raised to a specific degree, D q , α and f(α) will not change with the increase of q. Through multiple experimental simulation and in view of related literature 28  The larger the variation of Rayleigh dimension (ΔD q ) is, the larger the multifractal spectrum width (Δα) will be. If Rayleigh spectrum is a straight line (ΔD q → 0), multifractal spectrum will be concentrated to one point and the study object will show single fractal 52 . Figure 3(a) indicated when q → +∞, the maximum probability plays a decisive role; when q → −∞, the minimum probability plays a decisive role. None of the five kinds of heavy metals in reclaimed soil showed a linear distribution, suggesting these elements have a multifractal feature and are suitable to adopt Mkriging method for spatial digital mapping 28,52 . Figure 3(b) and Table 2 indicated ΔD q of the five kinds of soil heavy metals was consistent with Δα in the multifractal spectrum. Cd has the largest ΔD q (ΔD q = 1.4461), followed by As (ΔD q = 1.0016), Hg (ΔD q = 0.8309), Ni (ΔD q = 0.6526) and Cr (ΔD q = 0.4687); Δα width also shows the same behavior, Cd (Δα=1.571) > As (Δα = 1.1087) > Hg (Δα = 0.928) > Ni (Δα = 0.7857) > Cr (Δα = 0.5485). The wider multifractal spectrum Δα means more irregular spatial distribution of heavy metals in this soil and more fractal. The five heavy metals in reclaimed soil of industrial and mining abandoned land are in the following descending order in terms of homogeneity: Cd, As, Hg, Ni and Cr. Digital mapping. Spatial mapping of heavy metal content in reclaimed soil was conducted based on ordinary Kriging and Mkriging Considering the characteristics of variation function, the setting of the search radius of Mkriging method varied accordingly. For different heavy metals in reclaimed soil, five boundary points are selected in 500~3000 m, the side length is set as 4000 m and the observation scale is 0.3~0.6 (empirical value). The output grid is 10 m × 10 m (Fig. 4). Figure 4 revealed that regardless of spatial mapping methods, the same heavy metal in reclaimed soil basically had a consistent predicted spatial distribution pattern. The heavy metal content in reclaimed soil of the study area was high in the middle part and low in the periphery. This pattern was caused by mining, reclamation and other human activities, as well as topography, landforms and other natural factors. Damage by mining and reclamation activities were the major influence factors. Higher concentrations of heavy metals are found in areas with lower topography and more sulfur slag heap. Historically, as a subsidiary area of mining, heavy metal content was relatively low due to its high topography and far away from the sulphur slag heap. Figure 4 indicated that based on Mkriging method, the fluctuation was stronger, and the depiction of local specific values was more meticulous. The ranges of predicted values with Mkriging of elements Cd, Hg, As, Cr and Ni in reclaimed soil were 0.08-6.25 mg/kg, 0.04-0. 57 Table 2. Multifractal characteristics of different kinds of heavy metal. Note: α max and α min is the maximum and minimum value of the singular exponent, respectively. D qmin and D qmax is the maximum and minimum value of the Rayleigh dimension, respectively, ΔD q is their difference, representing the degree of variation of Rayleigh dimension. Δα expresses multifractal spectrum width.  Effect evaluation. The effects of different methods on spatial prediction of heavy metal content in reclaimed soil of industrial and mining abandoned land were comprehensively evaluated form three aspects including specific value reproduction effect, prediction accuracy and model simulation effect evaluation.
Comparison of reproduction effect of specific values. The difference in damage background and reclamation measures resulted in mutability and variability of spatial distribution of heavy metals. At present, the mostly frequently used spatial mapping method was ordinary Kriging method, but ordinary Kriging method alone had a smoothing effect and was unable to highlight the mutability of quality elements of reclaimed land resulting from the difference of reclamation measures. Therefore, it was crucial for a spatial digital mapping method to accurately depict the mutability of heavy metals. We used the coverage diagram and coverage ratio of specific value (CRSV) to describe the coincidence of different methods between predicted and measured specific values. We chose 20% of the minimum or maximum measured and predicted values as statistical magnitudes. CRSV value showed the proportion of their overlap. The minimum CRSV value of the same heavy metal element in reclaimed soil showed obvious differences among different methods. Regardless of the types of heavy metals in reclaimed soil, the minimum CRSV values with Mkriging method were 83.31%, 91.67%, 91.67%, 74.7% and 91.67% for five heavy metals Cd, Hg, As, Cr and Ni, respectively (Fig. 5). The measured and predicted values of Hg, As and Ni had the highest degree of coincidence. The predicted values based on Mkriging could more effectively reproduce the minimum specific values of the measured values. In comparison, kriging method showed a lower degree of coincidence in the minimum specific values of different heavy metals, and the CRSV value of element Cd did not exceed 60% though highest, while the lowest CRSV of element Cr was only 25%.
Heavy metal elements Cd, Hg and As in reclaimed soil have similar spatial distribution patterns, and showed a feature of being concentrated and mainly distributed in the western part and northern part of the area, i.e., in the perimeter zone of the study area, consistent with the pattern of the spatial distribution diagram as shown in Fig. 4. Before reclamation, these regions were mostly the land for auxiliary facilities such as factories and dormitories of mining and the pollution was less serious. During reclamation, soil covering and other protective measures were adopted.
There are greater differences between the two spatial prediction methods on the reproduction of the maximum specific value (Fig. 6). The maximum predicted values of elements Cd and Ni fully coincide with measured values, and the CRSV value was as high as 100%. The CRSV values were 100%, 75%, 83.33%, 91.44% and 100% respectively, while those of heavy metal elements Cd, Hg, As, Cr and Ni in soil based on ordinary Kriging were 50%, 41.67%, 25%, 50% and 41.67% respectively, mostly only one half of the values based on Mkriging method. The CRSV value of element As was less than one third of the value based on Mkriging method. As far as the reproduction of prediction results to specific values (mutability) was concerned, the spatial digital mapping of heavy metals in reclaimed soil of industrial and mining abandoned land based on Mkriging was more scientific and comprehensive.   Scatter diagrams of the measured and predicted values of different heavy metals in reclaimed soil under different methods were drawn (Fig. 7) and the linear relations of predicted values under different methods were fit to compare the degree of deviation from 1:1 bisector 43,45,46 . Figure 7 indicated that the same method showed different prediction accuracy and model fitting effects among different heavy metals in reclaimed soil, and the two methods showed a more significant difference in spatial prediction accuracy and model fitting effect of a same heavy metal in reclaimed soil. Regardless of the types of heavy metals, the predicted values based on kriging method all tended to be the same (close to the average of original data) and the fitted lines tended to be parallel, i.e., showed central and smoothing effects. In addition, the fitted lines of the predicted values of heavy metals in reclaimed soil based on kriging were blow 1:1 bisector. A larger included angle (acute angle) between the fitted line and bisector meant a greater downward deviation. When it was close to horizon, the prediction effect was the poorest. The fitted line of element Ni was basically parallel with the x-axis, and the predicted values tend to the same. Mkriging method reproduces specific values more satisfactorily 28,30 , the fitted line is closer to 1:1 bisector, and the fitted line of Cd basically coincides with the bisector. The fitted lines of predicted values of all heavy metals in reclaimed soil with Mkriging method were all above 1:1 bisector, and predicted values were generally larger, but the degree of deviation was much smaller than OK method. The relation between the fitted curve and bisector also revealed Mkriging was more capable to reflect the spatial structure of original regionalized variables and the mutual relations of variables.

Comparison of prediction accuracy and model simulation
The Kriging were larger, which indicated that the prediction accuracy is low and the fitting effect is poor. Mkriging was just the opposite, The RMSE and MSDR values were generally smaller, showing higher prediction accuracy and a better model fitting effect (Fig. 7). Compared with Kriging method, the RIs of RMSE of elements Cd, Hg, As, Cr and Ni with Mkriging were 95.28%, 61.74%, 78.54%, 82.51% and 83.58% respectively, and the spatial prediction accuracy of heavy metal element Cd in reclaimed soil based on Mkriging method is almost doubled.
Combining CRSV, RMSE and MSDR can not only reflect the characteristics of original samples and maintain the spatial variability and mutability of original samples but also guarantee prediction accuracy and model fitting effect of different methods. After comprehensive consideration of reproduction of specific values, RMSE and MSDR as well as measured and predicted scatter diagrams, the effect of spatial prediction based on Mkriging method was obviously higher than OK method for soil heavy metals in the study area. OK method was a low-pass filter and unable to rebuild the high frequency, local and weak signals of original signals 28,50 . The deletion of local variation information may result in the loss of some useful data 28,48,50 .
Mkriging method had a better overall effect, which may highlight specific values, and more satisfactorily maintain and reproduce the spatial structure of original data-based variables 28,48 . The effect of spatial prediction with Mkriging also showed certain difference among different kinds of heavy metals, and element Cd showed the best effect, followed by Hg, As, Ni and Cr, which was consistent with the fractal degrees, that was to say, the larger the fractal degree is, the more accurate the prediction will be. Element Cr showed a relatively poor prediction effect because it probably had a feature of tending to single fractal. Mkriging method may be used as one of the scientific digital mapping methods for heavy metals in reclaimed soil of industrial and mining abandoned land.

Conclusions and Discussions
We revealed the advantages of Mkriging method, provides methodological guidance for spatial prediction of heavy metal content in reclaimed soil of industrial and mining abandoned land.
(1) For reasons of history and reclamation, the original samples of heavy metals in reclaimed soil of industrial and mining abandoned land showed a very large range and variation degree, the C 0 /(C 0 + C 1 ) values of different heavy metals basically were all greater than 50%, random factors played a dominant role, and spatial variability mainly steamed from soil covering, conditioning measures of soil pH and other random factors. Such characteristics of original data decided the Kriging method cannot be used for spatial digital mapping. (2) Based on Mkriging method, the fluctuation was stronger, the depiction of local specific values was more meticulous and closer to the ranges of real samples of heavy metals in reclaimed soil. The Mkriging method was more capable to reflect the spatial structure of original regionalized variables and the mutual relations of variables. In comparison, the predicted values with Kriging are closer to the average of all real samples of heavy metals in reclaimed soil and had obvious smoothing and central effects. (3) After comprehensive consideration of CRSV, RMSE and MSDR as well as measured and predicted scatter diagrams, regardless of the types of heavy metals in soil, the effect of spatial digital mapping with Mkriging was obviously higher than Kriging method. Mkriging method may be used as one of the scientific methods predicting elements Cd, Hg, As, Cr and Ni and other heavy metals in reclaimed soil of industrial and mining abandoned land. Mkriging method may highlight specific values, and more satisfactorily maintain and reproduce the spatial structure of original data-based variables. The larger the fractal degree of a heavy metal in reclaimed soil is, the more accurate the prediction will be. The Mkriging method may interpolate irregularly distributed time-space signals into regularly distributed signals, which may extract the high frequency, local and weak signals of spatial signals. It was estimated that the process parameters may be used as characteristic parameters for mode recognition. The convolution filtering theory was applied to quantitatively derive the low-pass filtering feature of geostatistics. During interpolation, high frequency, local and weak signals were lost. After defining the measurement scale and measure of spatial signals, the article realized multidifractal interpolation, which retained more high frequency information in the system 56 (4) The combination of CRSV with RMSE and MSDR is an effective way to evaluate the effect of spatial mapping, which can not only reflect the spatial variability and mutability of original samples but also clearly state the prediction accuracy and model fitting effect of different methods. (5) The analysis of this article also indicated that Mkriging showed a good overall prediction effect of heavy metals in reclaimed soil of industrial and mining abandoned land, but its prediction effect on element Cr was poor probably due to the small sample size. A more suitable method needed to be selected in combination with other relevant methods, including Empirical Bayes Kriging method.