Next Article in Journal
Hypertension Knowledge-Level Scale (HK-LS): A Study on Development, Validity and Reliability
Previous Article in Journal
Social Aspects of Suicidal Behavior and Prevention in Early Life: A Review
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Combining Geostatistics with Moran’s I Analysis for Mapping Soil Heavy Metals in Beijing, China

1
Beijing Academy of Agriculture and Forestry, Beijing 100089, China
2
Department of Environmental Engineering, Taiyuan University, Taiyuan 030009, China
3
Department of Land Resources and Management, College of Natural Resources and Environment Science, China Agricultural University, Beijing 100193, China
*
Author to whom correspondence should be addressed.
Int. J. Environ. Res. Public Health 2012, 9(3), 995-1017; https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph9030995
Submission received: 4 January 2012 / Revised: 4 March 2012 / Accepted: 6 March 2012 / Published: 19 March 2012

Abstract

:
Production of high quality interpolation maps of heavy metals is important for risk assessment of environmental pollution. In this paper, the spatial correlation characteristics information obtained from Moran’s I analysis was used to supplement the traditional geostatistics. According to Moran’s I analysis, four characteristics distances were obtained and used as the active lag distance to calculate the semivariance. Validation of the optimality of semivariance demonstrated that using the two distances where the Moran’s I and the standardized Moran’s I, Z(I) reached a maximum as the active lag distance can improve the fitting accuracy of semivariance. Then, spatial interpolation was produced based on the two distances and their nested model. The comparative analysis of estimation accuracy and the measured and predicted pollution status showed that the method combining geostatistics with Moran’s I analysis was better than traditional geostatistics. Thus, Moran’s I analysis is a useful complement for geostatistics to improve the spatial interpolation accuracy of heavy metals.

1. Introduction

Heavy metal pollution in agricultural soils is becoming an urgent problem worldwide, due to increasing intensive anthropogenic activities, such as the discharge of wastes from metal processing plants, burning of fossil fuels and pesticide use. Excessive accumulation of heavy metals in agricultural soils can also be a source of pollution of surface and ground waters, living organisms, sediments, and oceans. Thus, mapping the spatial distribution of heavy metals in soils is critical for risk assessment of potential environmental pollution and for establishing protocols for pollution remediation, in particular, for China, with the recent three decades of intense economic development, the soil heavy metals’ pollution is now in a high risk state.
Geostatistics, providing a technique of semivariance to quantify the spatial patterns of soil parameters, is being increasingly adopted for spatial pattern analysis of heavy metals [1,2]. One criterion for evaluating the spatial dependence of spatial variables is imparted by the nugget/sill ratio [3,4], without the significance test of spatial dependence. In geostatistics, the active lag distance specifies the range over which semivariance can be calculated, which is usually about half of the maximum separation distance. However, this is only an empirical method [5]. In practice, due to the complexity of spatial data, it is difficult to perform iterative trial runs to determine a suitable active lag distance in order to generate a relatively stable and better-fitting theoretical semivariance without a priori knowledge [6]. In addition, spatial outliers in a dataset often make the semivariogram exhibit erratic behavior [7]. As a result, outliers are often deleted from spatial predictions [6,7]. However, in soil heavy metal evaluations, it is dangerous to ignore outliers, as these may actually represent potentially severely pollution areas [8].
Spatial autocorrelation analysis is another alternative method that has been widely used to explore the spatial pattern of variables in many fields [9,10,11]. Moran’s I is used to estimate the strength of the spatial correlation, and the significance of the spatial correlation can be tested [7,11]. Furthermore, spatial autocorrelation analysis can identify spatial clusters (positive autocorrelation) and spatial outliers (negative autocorrelation) of a regionalized variable [12]. Consequently, Huo et al. adopted Moran’s I to describe the spatial distribution pattern of heavy metals in Beijing cultivated soils [8].
Although spatial autocorrelation analysis cannot be used for estimation of unsampled areas, we recognized it could provide useful information for spatial variable mapping, if combined with a kriging method, for the production of high quality distribution maps. Therefore, taking heavy metals in Beijing agricultural soils as a case study, the primary objectives of this research were: (1) to compare spatial autocorrelation analysis and geostatistics for identifying the spatial pattern of heavy metals; and (2) to use the Moran’s I analysis results as the a priori knowledge for geostatistical interpolation, and to evaluate their effects on the estimation accuracy of heavy metal kriging mapping, focusing particularly on their influence on pollution status estimations.

2. Materials and Methods

2.1. Soil Sampling and Analysis

To investigate the pollution status of heavy metals in Beijing agricultural areas, a large-scale soil sampling project was conducted after the crop harvest in the autumn of 2006. According to the agricultural land distribution and land use type maps of Beijing, a non-uniform distribution of the stratified sampling technique was adopted to collect samples and ensure the representativeness of samples. The sampling strategy was divided into three steps to collect a total of 1,018 samples. First, 231 soil samples were collected from the entire study area, with uniform sampling being the low sampling density (C). Secondly, another 360 soil samples were added from areas with more agricultural soils to create the medium sampling density (M). Third, 427 soil samples were further collected on the basis of the two previous samplings and the agricultural soils to make a high sampling density (F). Figure 1 shows the distribution of soil samples at the three sampling densities.
Figure 1. The distribution of soil samples at three levels.
Figure 1. The distribution of soil samples at three levels.
Ijerph 09 00995 g001
For each sample, five surface soil (0~20 cm) sites were sampled within 10 × 10 m square areas and then mixed. A Global Positioning System was used to precisely locate each sampling position (latitude and longitude); and a total of 1 kg of mixed soil per sample was collected. All soil samples were collected using a stainless steel spade and a scoop made from bamboo and then stored in polyethylene bags. The soil samples were air-dried, crushed in an agate mortar, and then passed through a 100-mesh nylon sieve. The concentrations of eight heavy metals, including Cr, Ni, Cu, Zn, As, Cd, Pb, and Hg, were analyzed in the soil samples following the Chinese Environmental Quality Standard for Soils (GB15618-1995). After digesting the samples with a mixture of HCl, HNO3 and HClO4, the Cr, Ni, Cu, and Zn concentrations were analyzed by flame atomic absorption spectrophotometry, Pb and Cd were analyzed by graphite furnace atomic absorption spectrophotometry, and the As concentration was determined by potassium borohydride-silver nitrate spectrophotometry. In addition, the Hg concentration was analyzed by cold atomic absorption spectrophotometry after the samples were digested with a mixture of H2SO4, HNO3 and KMnO4. During processing, all samples were handled carefully to avoid input or loss of trace elements during preparation and analysis.

2.2. Spatial Autocorrelation

Spatial autocorrelation is an assessment of the correlation of a variable in reference to spatial location of the variable, which is a match between location similarity and attribute similarity [13]. Moran’s I is the more popular test statistic for spatial autocorrelation. Global Moran’s I examines whether spatial correlation exists or not over an entire region, and it is calculated as follow as [14]:
Ijerph 09 00995 e001
where Ijerph 09 00995 i003 is the number of observations of the whole region, Ijerph 09 00995 i004 and Ijerph 09 00995 i005 are the observations at locations of Ijerph 09 00995 i006 and Ijerph 09 00995 i007, Ijerph 09 00995 i008 is the mean of Ijerph 09 00995 i009, and Ijerph 09 00995 i010, an element of spatial weights matrix Ijerph 09 00995 i011, is the spatial weight between locations of Ijerph 09 00995 i006 and Ijerph 09 00995 i007. The value of Moran’s I generally varies between 1 and −1. Positive autocorrelation in the data translates into positive values of Moran’s I; negative autocorrelation produces negative values. No autocorrelation results in a value close to 0 [11].
The selection of neighbors is formally specified in the weights matrix Ijerph 09 00995 i011, which depicts the relationship between an element and its surrounding elements. A distance-based weight matrix was adopted in this study, and each distance class is specified as a threshold distance, such that all locations within the given distance are considered to be “neighbors” (the value not equal to zero in the matrix) in the distance-based weight matrix. Usually, normal approximation as a precondition, global Moran’s I can be standardized to Ijerph 09 00995 i012, and Ijerph 09 00995 i012 is calculated as [14,15]:
Ijerph 09 00995 e002
where Ijerph 09 00995 i019 is the sum of all weights located in the row Ijerph 09 00995 i006, Ijerph 09 00995 i020 is the sum of all weights in the column Ijerph 09 00995 i006. The threshold of 1.96 was applied to test the significance level of Ijerph 09 00995 i012. If Ijerph 09 00995 i012 was greater than 1.96 or smaller than −1.96, it implied that the spatial autocorrelation was significant [7].
The spatial correlogram is a graph where the Moran’s I is plotted in ordinate, against distances among localities (in abscissa). According to Legendre and Fortin, the spatial correlogram can be standardized into a standardized correlogram, in which the ordinate is the standardized Moran’s I, Z(I) [16]. The shape of the standardized correlogram provides indications about the spatial pattern (spatial clusters and spatial outliers) and spatial correlation distance of a variable [9,16]. However, the standardized correlogram often has one or more positive correlation ranges. Zhang et al. explained that the closer positive correlation range represents the average size of the zone of spatial clusters, that is, the spatial correlation distance [17].
Local Moran’s I is a local test statistic for spatial autocorrelation, which is used to identify the locations of spatial clusters and spatial outliers. It is computed as follows:
Ijerph 09 00995 e003
The notations in Equation (3) are as described for Equation (1), but the corresponding values are from the local neighboring region. For more details of Moran’s I principles and methods, see the references [12,14]. With the local Moran’s I statistic analysis, five categories of local spatial autocorrelation can be distinguished. Two of these are spatial clusters, including high values surrounded by high values (High-high), and low values surrounded by low values (Low-low) types. Two are spatial outliers, including high values surrounded by low values (High-low) and low values surrounded by high values (Low-high). The last type is the spatial randomness that is without significant spatial patterns at corresponding weight matrix.

2.3. Geostatistics Method

Geostatistics uses the technique of semivariance to measure the spatial variability of a regionalized variable, and provides input parameters for the spatial interpolation of kriging [18]. It relates the semivariance, γ(h), which is computed as half the average squared difference between the components of data pairs [19]:
Ijerph 09 00995 e004
where Ijerph 09 00995 i024 is the number of data pairs separated by a distance Ijerph 09 00995 i025, and Ijerph 09 00995 i026 is the measured value of the variable Ijerph 09 00995 i027 at location of Ijerph 09 00995 i028, and Ijerph 09 00995 i029 is the measured value of the variable Z at location Ijerph 09 00995 i030. For irregular sampling, it is rare for the distance between the sample pairs to be exactly equal to h, therefore, the lag distance h is often represented by a distance band.
A variogram plot can be acquired by calculating variogram at different lags. Data pairs were grouped into lag “bins” and Equation (4) was used to calculate the variogram for that bin. The mean lag of all the pairs in a particular bin was used as the representative lag for that bin.
The variogram plot is fitted with a theoretical model, such as spherical, exponential, Gaussian, linear and power models. In this study, the exponential model was selected.
The exponential function is:
Ijerph 09 00995 e005
where Ijerph 09 00995 i033 is the nugget variance ( Ijerph 09 00995 i034), represents the experimental error and field variation within the minimum sampling spacing. Typical, the semivariance increase with increasing lag distances to approach or attain a maximum value or still ( Ijerph 09 00995 i035) equivalent to the population variance. C is the structural variance and a is the spatial range across which the data exhibit spatial correlation.
Due to the complexity of spatial data, its spatial variability usually needs be described using two or more theoretical semivariances. This is the so-called nested model, which is described by the following equation [20]:
Ijerph 09 00995 e006
where Ijerph 09 00995 i037 is the nugget value of the nested model, which is usually considered to be the spatial variability that cannot be described at the smallest sampling scale, and Ijerph 09 00995 i038 is the semivariance at different scales.
The fitted model provides information about the spatial structure as well as input parameters for kriging interpolation. Kriging is a linear interpolation technique that provides a linear unbiased estimate for spatial variables, which can be depicted as follows:
Ijerph 09 00995 e007
where Ijerph 09 00995 i040 is the predicted value at location of Ijerph 09 00995 i041, Ijerph 09 00995 i042 is the measured value at location of Ijerph 09 00995 i028, Ijerph 09 00995 i003 is the number of sites within the search neighborhood used for the estimation. Contrary to other methods (such as inverse distance weighted), the weighting function Ijerph 09 00995 i043 is no longer arbitrary; rather, it is calculated based on the parameters of the semivariogram model. To ensure that the estimate is unbiased, the weights need to sum to one:
Ijerph 09 00995 e008
and the estimation errors (or kriging variances) need to be minimized.
With wide and increasing applications of spatial interpolation methods, there is a growing concern about their accuracy and precision. Accuracy of spatial interpolation was evaluated through cross-validation approach. Commonly used error measures include: mean error (ME), mean absolute error (MAE), mean squared error (MSE) and root mean squared error (RMSE). Willmott suggests that MAE and RMSE are among the “best” overall measures of model performance [21]. RMSE provides a measure of error size, but it is sensitive to outliers as it places a lot of weight on large errors [22]. MAE provides an absolute measure of the size of the error, and it is less sensitive to extreme values [23].
MAE is calculated as:
Ijerph 09 00995 e009
RMSE can be calculated as:
Ijerph 09 00995 e010
Because MAE does not reveal the magnitude of error that might occur at any point, MSE will be calculated [24] as:
Ijerph 09 00995 e011
where Ijerph 09 00995 i042, Ijerph 09 00995 i048 is the observed and predicted value at location of Ijerph 09 00995 i006, respectively. Smaller MAE values indicate few errors. The RMSE can be used to compare different methods by seeing how closely predicted values match the observed values, the smaller the RMSE the better. Squaring the difference at any point gives an indication of the magnitude. Smaller MSE values indicate more accurate point-by-point estimation.

2.4. Data Treatment with Computer Software

Soil samples were stored using the ArcView 3.2 software to create a spatial database. The spatial autocorrelation analysis was conducted using Geoda095i software. The experimental semivariance models were constructed using GS+5.3, while kriging was performed using the geostatistical analyst extension of ArcGIS 8.3.
As in conventional statistics, a normal distribution for the variable under study is desirable in linear geostatistics. Even though normality may not be strictly required, serious violation of normality, such as too high skewness and outliers, can impair the variorum structure and the kriging results. It is often observed that environmental variables are lognormal [1], and data transformation is necessary to normalize such data sets.
The normality tests of the eight heavy metals for the 1,018 samples were performed as described by Huo et al. [25]. Compared with the raw data sets, the log-transformation significantly reduced the skewness values of data sets towards “0”. However, because of data sets with many duplicate values or multi-modals, Cu, As, Cd, and Pb still cannot pass the Kolmogorov-Smirnov test for normality (K-S p) after log-transformation. As a result, further analysis focused on Cr, Ni, Zn, and Hg. In order to compare with geostatistics, the data after log-transformation were used in the Moran’s I analysis.

3. Results and Discussion

3.1. Spatial Pattern Analysis of Heavy Metals Using the Spatial Autocorrelation Analysis

In general, the higher the absolute value of Moran’s I is, the stronger a spatial autocorrelation exists, and the larger the absolute value of standardized Moran’s I is, the more significant a spatial structure exists. Figure 2 and Figure 3 represent the raw and standardized spatial correlograms of Cr, Ni, Zn, and Hg, respectively. Table 1 and Table 2 show the spatial autocorrelation characteristics of the four heavy metals, and the critical distance of weight matrix is the distance where the Moran’s I and the standardized Moran’s I, Z(I) reached a maximum, respectively.
Figure 2. Raw spatial correlograms of heavy metals (a) Cr, Zn, (b) Ni, Hg.
Figure 2. Raw spatial correlograms of heavy metals (a) Cr, Zn, (b) Ni, Hg.
Ijerph 09 00995 g002
Figure 3. Standardized spatial correlograms of heavy metals (a) Cr, Zn, (b) Ni, Hg.
Figure 3. Standardized spatial correlograms of heavy metals (a) Cr, Zn, (b) Ni, Hg.
Ijerph 09 00995 g003
The advantage of the standardized Moran’s I is that it can compare the significant spatial patterns of different variables or of the same variable with different calculating parameters. At the global level, Table 1 indicates Hg did not pass the significance test for the global standardized Moran’s I (1.96), and Zn was just at the significance level. Table 2 shows the four metals all pass the significance test, and Cr had the largest spatial dependence, followed by Ni, Hg. Zn had the weakest spatial structure among the four metals.
Compared with the local Moran’s I statistical analyses in Table 2, there were three major variations in Table 1. The first involved spatial clusters and High-low spatial outlier. Their percent considerably decreased for the four heavy metals, particularly Cr and Hg, whereas the absolute value of both Moran’s I and its standardized value for the three types increased distinctly for the four heavy metals. The second related to the spatial randomness (No significant spatial pattern), which increased considerably in percentage for the four heavy metals. These changes implied that spatial clusters and High-low spatial outlier became stronger and more significant; thus, the cores of these local spatial types were highlighted by reducing their structure percent. The third difference was that the Low-high spatial outlier became less significant.
Table 1. Spatial autocorrelation characteristics of the four heavy metals at global and local levels based on the distance where the Moran’s I reached maximum.
Table 1. Spatial autocorrelation characteristics of the four heavy metals at global and local levels based on the distance where the Moran’s I reached maximum.
Heavy metals Local spatial correlation typeGlobal
No significanceHigh-highLow-lowLow-highHigh-low
CrMoran’s I0.08720.84980.7035−0.1428−0.27550.4801
Standardized Moran’s I0.70406.78985.6228−1.1319−2.19063.8396
Percent of spatial types56.0914.3422.23.054.32-
NiMoran’s I0.16610.99940.7527−0.1724−0.43860.3173
Standardized Moran’s I1.13186.77565.1044−1.1608−2.96352.1558
Percent of spatial types69.947.0712.487.962.55-
ZnMoran’s I0.11650.73150.8123−0.1152−0.7820.2924
Standardized Moran’s I0.78714.90505.4464−0.7650−5.23001.9648
Percent of spatial types66.78.7413.467.563.54-
HgMoran’s I0.07420.99580.7823−0.228−0.39710.2725
Standardized Moran’s I0.52236.92795.4441−1.5775−2.75291.9009
Percent of spatial types67.789.6311.38.352.95-
Table 2. Spatial autocorrelation characteristics of the four heavy metals at global and local levels based on the distance where the standardized Moran’s I reached maximum.
Table 2. Spatial autocorrelation characteristics of the four heavy metals at global and local levels based on the distance where the standardized Moran’s I reached maximum.
Heavy metals Local spatial correlation typeGlobal
No significanceHigh-highLow-lowLow-highHigh-low
CrMoran’s I0.02960.39060.3731−0.2065−0.16560.3333
Standardized Moran’s I0.39515.05614.8293−2.6525−2.12454.3158
Percent of spatial types28.0919.4532.816.1913.46-
NiMoran’s I0.03510.46430.4287−0.1674−0.21480.2444
Standardized Moran’s I0.36854.74864.3857−1.6979−2.18222.5046
Percent of spatial types55.3014.8320.735.014.13-
ZnMoran’s I0.08590.46290.5584−0.2667−0.43080.2367
Standardized Moran’s I0.74223.96344.7794−2.2698−3.67222.0308
Percent of spatial types56.5813.1618.275.806.19-
HgMoran’s I0.00640.58120.3923−0.2215−0.24330.2054
Standardized Moran’s I0.07806.10294.1226−2.3114−2.54012.1637
Percent of spatial types49.9016.0120.535.308.25-
The disagreements in the spatial autocorrelation characteristics in Table 1 and Table 2 suggested that the maximum of a raw spatial correlogram can provide a suitable distance for detecting the local highlights of local spatial pattern, and it also agreed with the law that the closer the distance the more similarity. In contrast, the maximum of a standardized spatial correlogram focused on the major structure, with a smoothing effect on the details.
For the four heavy metals, their standardized spatial correlograms had more than one distinct waveform (Figure 3). Figure 3 provides the standardized Moran’s I values of Cr, which were positive at a distance from 1 km to 57 km, 79 km to 140 km. This indicated spatial clusters of similar Cr concentrations at these distance ranges. The standardized Moran’s I values of Cr were negative at a distance from 60 km to 72 km, which implied clusters of dissimilar Cr concentrations; that is, spatial outliers. Similarly, Zn and Hg all showed spatial clusters and spatial outliers over an entire region.
For the four metals, the amplitudes of spatial clusters were larger than for spatial outliers, indicating that positive spatial autocorrelation dominated at the global level. The same characteristics of the raw and standardized spatial correlograms were the distances where the 0 value first appeared, which were 57 km, 75 km, 57 km, and 55 km for Cr, Ni, Zn, and Hg, respectively (Figure 2, Figure 3). Thus, these were the maximal spatial positive correlation ranges of the corresponding heavy metals.

3.2. Spatial Pattern Analysis of Heavy Metals with Geostatistics

Table 3 lists the attributes of the semivariograms for the four heavy metals, their spatial correlation distances were identified as 16.5 km, 20 km, 20 km, and 55 km through trial and error, respectively. The semivariograms of Cr, Ni, Zn, and Hg were well fitted with the exponential model, as indicated by regression coefficients (R2) greater than 0.9, which indicated that Cr, Ni, Zn, and Hg had stronger spatial structure. The nugget/sill ratio of Cr, Ni, Zn, and Hg ranged from 34.2% to 48.9%, suggesting moderate spatial dependence, which indicates that the spatial variability of Cr, Ni, Zn, and Hg may be affected by intrinsic (soil formation factors, such as soil parent materials) and extrinsic factors (soil management practices, such as fertilization and pesticides). The spatial dependence of Cr is the strongest, followed by Ni, Zn, and then Hg. There was agreement in the spatial dependence of Cr and Ni metals between the geostatistics analysis and the spatial autocorrelation analysis, but disagreement for Zn and Hg. Comparing the significant spatial pattern of different variables, spatial autocorrelation analysis is preferred over geostatistics in spatial dependence because semivariance cannot test the significance of spatial dependence.
The ranges (spatial correlation distances) of Cr, Zn, and Hg are around 60 km, which indicated the spatial distribution of these three heavy metals may be similar, and the source may be the same. The spatial variability of Hg was significant, which indicated that the concentrations of Hg in soils were mainly affected by random factors (human activities).
The spatial correlation distances of Cr, Ni, Zn, and Hg from geostatistics were 59.55 km, 94.50 km, 65.79 km and 65.10 km, respectively (Table 3). These were larger than the corresponding ranges obtained from spatial autocorrelation analysis. The difference arose in part from the fact that geostatistics includes a mixture of positive and negative correlation [17].
Table 3. Semivariogram models for Cr, Ni, Zn, and Hg and their parameters (range, km).
Table 3. Semivariogram models for Cr, Ni, Zn, and Hg and their parameters (range, km).
Heavy metalsModelNugget (C0)Sill (C0 + C)Range (A0)Nugget/sill (C0/(C0 + C))/%R2RSS
CrExponential0.02510.073359.5534.20.9801.21 × 10−5
NiExponential0.05960.142394.5041.90.9723.52 × 10−5
ZnExponential0.03770.080165.7947.10.9303.80 × 10−5
HgExponential0.50101.025065.1048.90.9695.20 × 10−3

3.3. The Comparative Analysis of Estimation Accuracy for Global Heavy Metals

According to the spatial correlograms, four representative distances were selected: Ijerph 09 00995 i050 and Ijerph 09 00995 i051 represented the distance where the Moran’s I and the standardized Moran’s I, Ijerph 09 00995 i012 reached maximum, respectively. Ijerph 09 00995 i052 was the maximal positive correlation distance for Cr, Ni, Zn, and Hg, and Ijerph 09 00995 i053 was the sum distance of the positive and negative correlation. Table 4 lists these characteristic distances. There was no Ijerph 09 00995 i053 for Ni because no distinct negative correlation range was found between the two positive correlation ranges.
Table 4. The four characteristic distances of Cr, Ni, Zn, and Hg (km).
Table 4. The four characteristic distances of Cr, Ni, Zn, and Hg (km).
Heavy metalsDistance Ijerph 09 00995 i054Distance Ijerph 09 00995 i055Distance Ijerph 09 00995 i056Distance Ijerph 09 00995 i053
Cr6165776
Ni41075-
Zn475778
Hg4115591
A variogram plot can be acquired by calculating variograms at different lags. For irregular sampling, the active lag distance is often represented by a distance band. Generally, the distance band was adjusted repeatedly for higher match between the theoretical model and the experimental semivariance. In this study, in order to effectively and quickly find the suitable active lag distance, we tried to use the distances parameters extracted from the Moran’s I analysis as an auxiliary tool. Therefore, the four characteristic distances were tested as the active lag distance to fit the semivariance and produce spatial interpolation, and these were labeled by model Ijerph 09 00995 i050, Ijerph 09 00995 i057, Ijerph 09 00995 i056, Ijerph 09 00995 i053, respectively. Because the Moran’s I analysis based on Ijerph 09 00995 i050 and Ijerph 09 00995 i057 can highlight the local (often the spatial outliers) and major spatial structure and the complexity of spatial variables, the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057 was used to describe the spatial variability of heavy metals at local and global scales. It is worth noting that model Ijerph 09 00995 i058 was the traditional geostatistics model in the following.
Take Cr as a case. Semivariogram γ(h) for Cr was calculated and the scatter plot of γ(h) vs. Ijerph 09 00995 i025 was generated. The different theoretical semivariance models were used to try with the best fitting value. The smallest nugget values of mean prediction error (ME) close to 0 and root-mean-square standardized prediction errors (RMSSE) close to 1 were selected. Figure 4 presents the scatter plots and fitted model based on the traditional geoststistics model and model Ijerph 09 00995 i059. According to the scatter plots, the values ME and RMSSE were calculated. The ME values of traditional geostatistics model and model Ijerph 09 00995 i059 were 0.234, 0.1015, and the RMSSE were 1.119 and 1.075, respectively. Therefore, the theoretical semivariance of model Ijerph 09 00995 i059 was better than traditional geostatistics model. For Ni, Zn, and Hg, The selection of the semivariance were similar to Cr.
Figure 4. The scatter plots and fitted model based on the traditional geoststistics model (a) and model Ijerph 09 00995 i059 (b) (distance, m).
Figure 4. The scatter plots and fitted model based on the traditional geoststistics model (a) and model Ijerph 09 00995 i059 (b) (distance, m).
Ijerph 09 00995 g004
The spatial interpolation maps of the four heavy metals were conducted using the ordinary kriging method based on model Ijerph 09 00995 i050, Ijerph 09 00995 i057, Ijerph 09 00995 i056, Ijerph 09 00995 i060, Ijerph 09 00995 i058, and the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057. The cross-validation method was used to validate the parameters of these models. Evaluation indices from cross-validation of spatial maps of the four heavy metals are given in Table 5. For Cr, Ni, Zn, and Hg, except for the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057, the MAE, RMSE, and MAE values of model Ijerph 09 00995 i050 and model Ijerph 09 00995 i061 were the smallest, which indicated that semivariance parameters obtained from fitting of experimental semivariogram values were reasonable and spatial prediction using semivariogram parameters is better than other models. The values of model Ijerph 09 00995 i052 and Ijerph 09 00995 i053 were close and the largest, implying unsuitable semivariogram parameters. The values of model Ijerph 09 00995 i058 were close to the values of model Ijerph 09 00995 i061, which were the third smallest. This suggested that the dominant spatial pattern can be captured through traditional geostatistics. Compared with model Ijerph 09 00995 i050, the values of the MAE, RMSE, and MAE of the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057 were smaller, which indicated that the nested model can further improve the global prediction. It may be caused by the model Ijerph 09 00995 i059 being able to highlight local spatial patterns while the model Ijerph 09 00995 i061 focuses on the global dominant spatial pattern. Thus, their coupling may allow fusion of the details and the dominating spatial patterns of heavy metals. Table 5 shown that the best model was the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057, followed by model Ijerph 09 00995 i050, Ijerph 09 00995 i057, Ijerph 09 00995 i062. Therefore, the Moran’s I analysis can help to provide the better semivariogram parameters and improve estimation accuracy.
Table 6 lists the statistics results of measured and predicted heavy metals concentrations. After spatial interpolation, the mean concentrations of Cr, Ni, Zn, and Hg for different models were similar to those of the measured values. The minimum values were larger than those of the measured value, while the maximum values, standard deviations, and CVs were smaller than those of the measured value. The result indicated that all four interpolation methods caused a degree of compression due to the smoothing effect of kringing interpolation. When the model Ijerph 09 00995 i059 and the nested model Ijerph 09 00995 i054 plus Ijerph 09 00995 i055 were used, the data ranges were wider, and the standard deviations and CVs were larger than when other models were used. This indicated that the model Ijerph 09 00995 i059 and the nested model Ijerph 09 00995 i054 plus Ijerph 09 00995 i055 can reduce the smoothing effect of kriging and improve the estimatiom accuracy, especially the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057.
Table 5. Evaluation indices of the interpolation maps of heavy metals.
Table 5. Evaluation indices of the interpolation maps of heavy metals.
Evaluation indicesCrNiZnHg
Model Ijerph 09 00995 i054MAE6.371.458.660.0935
RMSE10.563.2012.560.2170
MSE111.5710.27157.660.0471
Model Ijerph 09 00995 i055MAE6.473.559.920.1195
RMSE10.817.2914.340.2683
MSE116.7653.19205.530.0720
Model Ijerph 09 00995 i063MAE7.755.4412.160.1147
RMSE12.739.9117.420.2540
MSE161.9798.20303.610.0645
Model Ijerph 09 00995 i060MAE7.66-12.530.1351
RMSE12.59-17.880.2817
MSE158.51-319.550.0794
Model Ijerph 09 00995 i058MAE7.064.8310.860.1160
RMSE11.699.1615.660.2557
MSE136.7383.86245.160.0654
Model Ijerph 09 00995 i054 plus Ijerph 09 00995 i055MAE5.311.177.670.0918
RMSE8.902.5911.090.2291
MSE79.196.73123.030.0525
Table 6. Statistics results of measured and predicted heavy metals concentrations (mg·kg−1).
Table 6. Statistics results of measured and predicted heavy metals concentrations (mg·kg−1).
Heavy metals MeanMinimumMaximumRangeStandard DdeviationCV (%)
CrMeasured value60.7531.60300.00268.4020.4933.73
Model Ijerph 09 00995 i05960.5739.45156.59117.1414.2823.57
Model Ijerph 09 00995 i06160.5639.93142.62102.6913.9022.95
Model Ijerph 09 00995 i05860.5040.69136.4195.7213.5022.32
Model Ijerph 09 00995 i054 plus Ijerph 09 00995 i05560.8238.19162.12123.9415.1824.95
NiMeasured value28.498.87203.38194.5111.2539.49
Model Ijerph 09 00995 i05928.4210.16139.08128.918.6430.42
Model Ijerph 09 00995 i06128.3913.4559.3445.895.8820.72
Model Ijerph 09 00995 i05828.4415.0444.7229.684.8617.10
Model Ijerph 09 00995 i054 plus Ijerph 09 00995 i05528.5510.10147.14137.049.1532.07
ZnMeasured value76.2728.50221.62193.1221.0327.57
Model Ijerph 09 00995 i05976.2645.44144.8299.3812.3116.14
Model Ijerph 09 00995 i06176.2047.76128.9381.1711.1814.67
Model Ijerph 09 00995 i05876.1449.99117.2567.2710.3713.62
Model Ijerph 09 00995 i054 plus Ijerph 09 00995 i05576.5543.75159.64115.8813.4317.55
HgMeasured value0.21750.00054.29004.28950.3210147.59
Model Ijerph 09 00995 i0590.21130.02191.62561.60370.160275.80
Model Ijerph 09 00995 i0610.20350.05090.88370.83280.117757.87
Model Ijerph 09 00995 i0580.20720.04851.13821.08970.132764.04
Model Ijerph 09 00995 i054 plus Ijerph 09 00995 i0550.21290.03141.15471.12330.151070.90
The above results shown that the characteristic distances provided by Moran’s I analysis are feasible for improving the spatial estimation accuracy, but the mathematical proof of the methodology was not explored here. In addition, ordinary kriging, the most basic and commom spatial interpolation method, was used to, other kriging model such as indicator kriging (which makes no assumption of normality) with Moran’s I analysis can be further examined for the possible.

3.4. The Impact on Pollution Status of Heavy Metals

In order to understand the impact of the spatial interpolation on the pollution status of heavy metals, a single-factor method was used to assess the pollution status, based on the critical value of “Chinese Environmental Quality Standard for Soils” (GB 15618-1995). The pollution status was classified, according to a single-factor pollution index, as unpolluted or polluted [26]. After spatial interpolation using model Ijerph 09 00995 i050, Ijerph 09 00995 i057, the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057, and the traditional geostatistics model Ijerph 09 00995 i058, the pollution status of both measured and predicted values of the samples were assessed. In Table 7, the rows are the pollution status based on predicted values, while the columns are the pollution status based on measured values. Thus, in each block composed by heavy metal and prediction model, data in the left diagonal present the agreement in pollution status between the measured and predicted values, while those in the right diagonal present the disagreement in pollution status between the measured and predicted values.
Table 7. The sample agreements in pollution status between ground measure and interpolation (%).
Table 7. The sample agreements in pollution status between ground measure and interpolation (%).
CrNiZnHg
PollutedUnpollutedPollutedUnpollutedPollutedUnpollutedPollutedUnpolluted
Model d1Polluted0.10 2.36 3.050.49
Unpolluted0.5999.311.5796.070.1099.903.2493.22
Model d2Polluted 1.280.29 1.870.39
Unpolluted0.6999.312.6595.780.1099.904.4293.32
Model d5Polluted 0.490.29 1.960.39
Unpolluted0.6999.313.4495.780.1099.904.3293.32
Model d1 plus d2Polluted0.10 3.05 3.240.59
Unpolluted0.5999.310.8896.070.1099.903.0593.12
For Zn, 0.1% of measured samples were in pollution status, and these became unpolluted after spatial interpolation using all models, although the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057 improved the Zn global prediction accuracy (Table 7). For Cr, 0.69% of measured values were in pollution status, and these samples were underestimated and became unpolluted after interpolation using models Ijerph 09 00995 i057 and traditional geostatistics model Ijerph 09 00995 i062. However, after interpolation using model Ijerph 09 00995 i059 and the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057, there were still 0.10% of samples remained pollution status. Before and after interpolation, the changes of samples pollution status were more complicated for Ni and Hg, as there was not only an underestimation effect, but also an overestimation effect. These two effects were severe after interpolation using models Ijerph 09 00995 i057 and Ijerph 09 00995 i062 (Table 7). For Ni, 3.93% of measured samples were in pollution status. After spatial interpolation, the model Ijerph 09 00995 i050 avoided the overestimation effect and reduced the underestimation effect from 3.44% to 1.57% compared with the model Ijerph 09 00995 i062. The nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057 showed that there were still 3.05% of samples remained pollution status, but only 0.88% of samples were changed from polluted to unpolluted status after interpolation.
For Hg, 6.29% of measured samples were in pollution status. Compared with the model Ijerph 09 00995 i062, the percent of samples remained pollution status increased from 1.96% to 3.05% after interpolation using the model Ijerph 09 00995 i050, and the percent of samples changed from polluted to unpolluted reduced from 4.32% to 3.24%. Thus, it avoided designating 1.08% of polluted samples wrongly as unpolluted. Models Ijerph 09 00995 i050 and Ijerph 09 00995 i062 also had the overestimation effect, the percent of samples changed unpolluted to polluted were 0.39% and 0.49%, respectively, therefore, 0.10% of unpolluted samples were wrongly identified as polluted. Compared with model Ijerph 09 00995 i050, the percent of samples remained pollution status increased by 0.19%, and the percent of samples changed from unpolluted to polluted increased by 0.10% after interpolaton using the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057.
The results of pollution status assessment of samples showed that, compared with the traditional geostatistics model Ijerph 09 00995 i058, the method combining geostatistics with Moran’s I analysis (models Ijerph 09 00995 i050, Ijerph 09 00995 i057) can avoid effectively the underestimation and overestimation effect so that more samples remain original pollution status.
In general, heavy metal polluted soil samples, as a small probability event, would be underestimated by interpolation, which is exemplified by Zn and Cr in this study. If the distance where the Moran’s I reached maximum ( Ijerph 09 00995 i059) was used to optimize the semivariance, this underestimation effect can be slightly reduced. With the increase of the polluted samples and a significant local spatial pattern appeared, the effect of reducing underestimation was more significant, such as Ni and Hg in this study. If the distances where the Moran’s I and the standardized Moran’s I, Ijerph 09 00995 i012 reached maximum ( Ijerph 09 00995 i059 and Ijerph 09 00995 i057) were combined to optimize the semivariance, the underestimation effect can be reduced further. The reason for this may be that the optimality of semivariance using Ijerph 09 00995 i059 can highlight the most significant local spatial patterns and reflects the more abnormal information of the polluted samples. Thus, the nested model would better balance the details and dominating spatial patterns. But, there also was an overestimation effect. Because most of the overestimation errors happened near actual polluted samples, making the potentially high-risk areas, in terms of soil heavy metal pollution, a little overestimation may be better than more underestimation. By comparing all the models Ijerph 09 00995 i050, Ijerph 09 00995 i057, Ijerph 09 00995 i058, the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057, the nested model was defined as the optimality model.
As mentioned previously, spatial autocorrelation analysis was adopted for the optimality of semivariance by simply deleting the spatial outliers. However, in pollution studies, this may cause severe hypercorrectness when the spatial outliers are in fact reasonable. If global dominating spatial patterns are the focus, then standardized spatial correlogram can provide the spatial correlation distance for the optimality of semivariance. In contrast, if abnormal situations or details are the key, such as evaluation of the soil heavy metal pollution, then raw spatial correlogram can provide the information to help the optimality of semivariance. Moreover, a nested model that fuses both the details and the dominating spatial patterns can provide an even better prediction.
In the current study, for Hg metal element, there still was a large gap of assessment accuracy between the nested model interpolation and the measured values. A greater improvement in assessment accuracy may occur if zonal geostatistics are interpolated according the spatial distribution of the local Moran’s I spatial pattern types. In addition, as this study primarily focused on the spatial autocorrelation analysis for the optimality of semivariance, the ordinary kriging type may not have been the optimal for the estimation accuracy.

3.5. The Spatial Distribution of Heavy Metals Using the Nested Model Interpolation

Figure 5 shows the spatial distribution of heavy metal concentrations interpolated using the optimality model, that is the nested model Ijerph 09 00995 i059 and Ijerph 09 00995 i057. Heavy metals concentrations are separated into classes according to the background values of soil heavy metals of Beijing and their multiples to highlight the spatial differences of different classes.
Figure 5. Distribution maps of heavy metals based on the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057 (a) Cr, (b) Ni, (c) Zn, (d) Hg.
Figure 5. Distribution maps of heavy metals based on the nested model of Ijerph 09 00995 i050 and Ijerph 09 00995 i057 (a) Cr, (b) Ni, (c) Zn, (d) Hg.
Ijerph 09 00995 g005
Cr concentrations in the study areas were greater than the background value (29.8 mg·kg−1), and areas where Cr concentrations were three times the background value were observed in northeast Beijing (Figure 5(a)). The areas where Ni concentrations were lower than the background value (26.8 mg·kg−1) were more, and areas where Ni concentrations were twice the background value were found in the Miyun, Shunyi, and Daxing districts (Figure 5(b)). Areas where Zn concentrations were lower than the background value (57.5 mg·kg−1) were in the Daxing and Fangshan districts, and the areas that Zn concentrations were twice of the background value were in the Daxing, Huairou, Pinggu, and Fangshan districts (Figure 5(c)). The areas that Hg concentrations were lower than background value (0.08 mg·kg−1) were in Daxing and Shunyi districts, and the areas that Hg concentrations exceeded 0.32 mg·kg−1 (four times of Hg background value) were primarily distributed in the urban fringe around the city, while a few areas with higher levels were also observed in the Miyun, Pinggu and Fangshan districts (Figure 5(d)). The background value was the original content of soil heavy metals which was not or rarely influenced by human activities. However, the spatial distribution of the four heavy metals showed that human activities had caused the content of heavy metals to be more than the background value and to accumulate. Three major areas were identified that should be received more attention. One is the northeast region of Beijing, where mining operations are responsible for the enrichment in heavy metals. The second is the southeast part of Beijing where wastewater irrigation has strongly changed the content of metals in soils, particularly those of Cr and Zn. The last is the urban fringe around the city, where Hg has had a significant increase.

4. Conclusions

Both geostatistics and spatial autocorrelation analysis can evaluate the spatial patterns of heavy metals. However, the two methods have their advantages and disadvantages. Geostatistics can provide a technique of semivariance to quantify the spatial patterns of soil parameters, but the fitting of variogram is influenced by subjective factors, and it will affect the kriging estimation. On the other hand, the Moran’s I analysis just can provide some spatial autocorrelation distances of variable, which have the same meaning as the range calculation from the variogram, so in this paper we tried to use this information to help calculate the semivariance in geostatistics and produce spatial interpolation to improve the accuracy of traditional geostatistics. This is the method combining geostatistics with Moran’s I analysis.
According to spatial correlogram of Moran’s I analysis, four characteristics distances were obtained and used as the active lag distance to calculate the semivariance. The resulted showed that the fitting accuracy of semivariance based on the distances where the Moran’s I and the standardized Moran’s I, Z(I) reached maximum were better than traditional geostatistics, because the Moran’s I analysis based on the two distance can detect the local spatial and major spatial pattern of heavy metals, respectively. Then, the spatial interpolation was produced based on the two distances. By comparing the values MAE, RMSE, and MSE, and the pollution status of measured and predicted values, the results showed that the estimation accuracy of the method combining geostatistics with Moran’s I analysis were better than the traditional geostatistics. In addition, because of the complexity of spatial data, the nested model, which is coupled by the distances where the Moran’s I and the standardized Moran’s I, Z(I) reached a maximum, was used to calculate the semivariance and produce spatial interpolation, the results showed that the nested model was the optimality and improve the fitting and estimation accuracy further. Consequently, the Moran’s I analysis can be used as a useful complement to geostatistics and produce a high quality spatial interpolation maps of heavy metals. Based on the interpolation maps produced using the nested model, areas of high concentrations of heavy metals can be identified, this is very important for risk assessment of environmental and pollution remediation.

Acknowledgements

The research was supported by the National Natural Science Foundation (41130526) and Beijing Municipal Bureau of Finance programs support. The authors are grateful to the Agricultural Environmental Monitoring Station of Beijing for their soil sampling and analysis, as well as the editor and three anonymous referees for their comments on earlier versions of the manuscript.

References

  1. McGrath, D.; Zhang, C.S.; Carton, O.W. Geostatistical analyses and hazard assessment on soil lead in Slivermines areas, Ireland. Environ. Pollut. 2004, 127, 239–248. [Google Scholar]
  2. Franco, C.; Soares, A.; Delgado, J. Geostatistical modelling of heavy metal contamination in the topsoil of Guadiamar river margins (S Spain) using a stochastic simulation technique. Geoderma 2006, 136, 852–864. [Google Scholar]
  3. Guo, X.D.; Fu, B.J.; Ma, K.M.; Chen, L.D.; Wang, J. Spatio-temporal variability of soil nutrients in Zunhua plain, northern China. Phys. Geogr. 2001, 22, 343–360. [Google Scholar]
  4. Sun, B.; Zhou, S.L.; Zhao, Q.G. Evaluation of spatial and temporal changes of soil quality based on geostatistical analysis in the hill region of subtropical China. Geoderma 2003, 115, 85–99. [Google Scholar]
  5. Western, A.W.; Bloschl, G.; Grayson, R.B. Geostatistical characterisation of soil moisture patterns in the Tarrawarra catchment. J. Hydrol. 1998, 205, 20–37. [Google Scholar]
  6. Huang, Z.G.; Li, B.G.; Hu, K.L. Characteristics of the spatio-temporal changes of soil organic matter of sugarcane field in red soil hill areas (in Chinese, with English abstract). Trans. CSAE 2006, 22, 58–63. [Google Scholar]
  7. Zhang, C.S.; McGrath, D. Geostatistical and GIS analyses on soil organic carbon concentrations in grassland of southeastern Ireland from two different periods. Geoderma 2004, 119, 261–275. [Google Scholar]
  8. Huo, X.N.; Zhang, W.W.; Sun, D.F.; Li, H.; Zhou, L.D.; Li, B.G. Spatial pattern analysis of heavy metals in Beijing agricultural soils based on spatial autocorrelation statistics. Int. J. Environ. Res. Public Health 2011, 8, 2074–2089. [Google Scholar] [CrossRef]
  9. Zhang, C.S.; Tao, S.; Yuan, G.P.; Liu, S. Spatial autocorrelation analysis of trace element contents of soil in Tianjin plain area (in Chinese, with English abstract). Acta Pedol. Sin. 1995, 32, 50–57. [Google Scholar]
  10. Anselin, L. Spatial effects in econometric practice in environmental and resource economics. Am. J. Agr. Econ. 2001, 83, 705–710. [Google Scholar]
  11. Overmars, K.P.; de Koning, G.H.J.; Veldkamp, A. Spatial autocorrelation in multi-scale land use models. Ecol. Model. 2003, 164, 257–270. [Google Scholar]
  12. Anselin, L. Local indicators of association—LISA. Geogr. Anal. 1995, 27, 93–115. [Google Scholar]
  13. Cliff, A.D.; Ord, J.K. Spatial Autocorrelation; Pion: London, UK, 1973; Volume 178. [Google Scholar]
  14. Cliff, A.D.; Ord, J.K. Spatial Processes: Models and Applications; Pion: London, UK, 1981. [Google Scholar]
  15. Goodchild, M.F. Spatial Autocorrelation. CATMOG 47; Geobooks: Norwich, UK, 1986; pp. 6–25. [Google Scholar]
  16. Legendre, P.; Fortin, M.J. Spatial pattern and ecological analysis. Plant Ecol. 1989, 80, 107–138. [Google Scholar]
  17. Zhang, C.S.; Zhang, S.; He, J.B. Spatial distribution characteristics of heavy metals in the sediments of Changjiang River system—Spatial autocorrelation and fractal methods (in Chinese, with English abstract). Acta Geogr. Sin. 1998, 53, 87–96. [Google Scholar]
  18. Webster, R.; Oliver, M.A. Geostatistics for Environmental Scientists; Wiley: Chichester, UK, 2001. [Google Scholar]
  19. Isaaks, E.H.; Srivastava, R.M. An Introduction to Applied Geostatistics; Oxford University Press: New York, NY, USA, 1989; pp. 140–398. [Google Scholar]
  20. Wang, Z.Q. Geostatistics and Its Application in Ecology (in Chinese); Science Press: Beijing, China, 1999; pp. 162–192. [Google Scholar]
  21. Willmott, C.J. Some comments on the evaluation of model performance. Bull. Am. Meteorol. Soc. 1982, 63, 1309–1313. [Google Scholar]
  22. Hernandez-Stefanoni, J.L.; Ponce-Hernandez, R. Mapping the spatial variability of plant diversity in a tropical forest: Comparison of spatial interpolation methods. Environ. Monit. Assess. 2006, 117, 307–334. [Google Scholar]
  23. Vicente-Serrano, S.M.; Saz-Sánchez, M.A.; Cuadrat, J.M. Comparative analysis of interpolation methods in the middle Ebro Valley (Spain): Application to annual precipitation and temperature. Clim. Res. 2003, 24, 161–180. [Google Scholar]
  24. Priyabrata, S.; Chopra, U.K.; Debashis, C. Spatial variability of soil properties and its application in predicting surface map of hydraulic parameters in an agricultural farm. Curr. Sci. 2008, 95, 937–945. [Google Scholar]
  25. Huo, X.N.; Li, H.; Sun, D.F.; Zhou, L.D.; Li, B.G. Multi-scale spatial structure of heavy metals in agricultural soils in Beijing. Environ. Monit. Assess. 2010, 164, 605–616. [Google Scholar]
  26. Chen, T.B.; Zheng, Y.M.; Lei, M.; Huang, Z.C.; Wu, H.T.; Chen, H.; Fan, K.K.; Yu, K.; Wu, X.; Tian, Q.Z. Assessment of heavy metal pollution in surface soils of urban parks in Beijing, China. Chemosphere 2005, 60, 542–551. [Google Scholar]

Share and Cite

MDPI and ACS Style

Huo, X.-N.; Li, H.; Sun, D.-F.; Zhou, L.-D.; Li, B.-G. Combining Geostatistics with Moran’s I Analysis for Mapping Soil Heavy Metals in Beijing, China. Int. J. Environ. Res. Public Health 2012, 9, 995-1017. https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph9030995

AMA Style

Huo X-N, Li H, Sun D-F, Zhou L-D, Li B-G. Combining Geostatistics with Moran’s I Analysis for Mapping Soil Heavy Metals in Beijing, China. International Journal of Environmental Research and Public Health. 2012; 9(3):995-1017. https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph9030995

Chicago/Turabian Style

Huo, Xiao-Ni, Hong Li, Dan-Feng Sun, Lian-Di Zhou, and Bao-Guo Li. 2012. "Combining Geostatistics with Moran’s I Analysis for Mapping Soil Heavy Metals in Beijing, China" International Journal of Environmental Research and Public Health 9, no. 3: 995-1017. https://0-doi-org.brum.beds.ac.uk/10.3390/ijerph9030995

Article Metrics

Back to TopTop