Effects of Vegetable Fields on the Spatial Distribution Patterns of Metal(loid)s in Soils Based on GIS and Moran’s I

Metal(loid) pollution in vegetable field soils has become increasingly severe and affects the safety of vegetable crops. Research in China has mainly focused on greenhouse vegetables (GV), while open field vegetables (OV) and the spatial distribution patterns of metal(loid)s in the surrounding soils have rarely been assessed. In the present study, spatial analysis methods combining Geographic Information Systems (GIS) and Moran’s I were applied to analyze the effects of vegetable fields on metal(loid) accumulation in soils. Overall, vegetable fields affected the spatial distribution of metal(loid)s in soils. In long-term vegetable production, the use of large amounts of organic fertilizer led to the bioconcentration of cadmium (Cd) and mercury (Hg), and long-term fertilization resulted in a significant pH decrease and consequent transformation and migration of chromium (Cr), lead (Pb), and arsenic (As). Thus, OV fields with a long history of planting had lower average pH and Cd, and higher average As, Cr, Hg, and Pb than GV fields, reached 0.93%, 10.1%, 5.8%, 3.0%, 80.8%, and 0.43% respectively. Due to the migration and transformation of metal(loid)s in OV soils, these should be further investigated regarding their abilities to reduce the accumulation of metal(loid)s in soils and protect the quality of the cultivated land.


Introduction
The term "metal(loid)s" represents both metals and metalloids. Metal(loid) accumulation can lead to the contamination of surface water, groundwater, organisms, sediments, and oceans. Metal(loid) pollution in agricultural soils has become an urgent issue worldwide. This is of particular concern in China due to its rapid economic growth over the past 40 years, and metal(loid) pollution in soils is now being considered of high-risk to the environment and human health [1]; the consequential environmental problems have also received widespread attention [2]. Some researchers have employed multivariate statistical analysis and geostatistical analysis for qualitative or quantitative research to determine the heavy metal sources in soils; however, it is generally believed that natural sources and human activities are the two major providers of heavy metals [3][4][5]. Natural sources include rock components [6][7][8], soil parent material [9], and atmospheric sediments from soil formation processes, and these are concentrated in soils after weathering and leaching, attaining high geological background values. Human activities mainly comprise industrial activities such as mining, industrial emissions, coal combustion, point source emissions [10][11][12][13][14][15], agricultural production from the long-term and massive application of fertilizers [16,17], and life activities.     From July to December 2017, 375 0-20 cm topsoil samples were collected across Feidong County ( Figure 1) using bamboo shovels. During the sampling, areas with new and locally contaminated soils were avoided, and the geographic coordinates of the sampling sites were recorded. After sampling, the GPS-measured sampling points with coordinate records were converted into points with spatial coordinates using ArcGIS 10.2 (Environmental Systems Research Institute, Inc, RedLands, CA, USA) and projection transformation was applied to generate a sample distribution map with soil metal(loid) information.

Vegetable Field Types
Among the vegetable field types, the greenhouse vegetables are an area covered with a greenhouse facility. The location and area of greenhouse facilities in 2019 were extracted by the four spatial attributes of area, length, rectangular_fit, and major_Length with the Rule-based Classification step in the ENVI software from the Google Earth Level 17 data, as shown on the left and middle in Figure 3. The overall accuracy was 77.439% and the kappa coefficient was 0.618, then the extraction results were modified manually. The open field vegetables are the vegetable fields without greenhouse facilities minus the greenhouse vegetable areas, to reduce the impact of different types of vegetable fields, as shown on the right in Figure 3. Distribution of the soil types in Feidong County (data from Anhui Provincial Agriculture Committee).
From July to December 2017, 375 0-20 cm topsoil samples were collected across Feidong County ( Figure 1) using bamboo shovels. During the sampling, areas with new and locally contaminated soils were avoided, and the geographic coordinates of the sampling sites were recorded. After sampling, the GPS-measured sampling points with coordinate records were converted into points with spatial coordinates using ArcGIS 10.2 (Environmental Systems Research Institute, Inc, RedLands, CA, USA) and projection transformation was applied to generate a sample distribution map with soil metal(loid) information.

Vegetable Field Types
Among the vegetable field types, the greenhouse vegetables are an area covered with a greenhouse facility. The location and area of greenhouse facilities in 2019 were extracted by the four spatial attributes of area, length, rectangular_fit, and major_Length with the Rule-based Classification step in the ENVI software from the Google Earth Level 17 data, as shown on the left and middle in Figure 3. The overall accuracy was 77.439% and the kappa coefficient was 0.618, then the extraction results were modified manually. The open field vegetables are the vegetable fields without greenhouse facilities minus the greenhouse vegetable areas, to reduce the impact of different types of vegetable fields, as shown on the right in Figure 3. The obtained spatial distribution was intersected with the spatial distribution of the total vegetable fields from Feidong's agricultural survey data in 2016 to obtain the spatial distributions of the OV and GV fields ( Figure 4).

Soil Sampling and Analysis
The soil tests were completed at the Anhui Institute of Geological Experiment. After the samples were digested by HCl-HNO3-HClO4-HF, the concentrations of Cd, Pb, and Cr were determined by inductively coupled plasma mass spectrometry (ICP-MS 7700, Agilent Technologies, Santa Clara, CA, USA); the concentrations of Hg and As were determined by atomic fluorescence spectrometry (AFS-8220, JiTian Technologies, Beijing, China) after digestion by aqua regia; the soil pH was measured using a pH meter (SG8, METTLER TOLEDO, Zurich, Switzerland). The laboratory used national standard materials and repeated sample tests for quality monitoring. The accuracy, precision, and reported percent of the various indicators were controlled at 0.10-0.12, 10-20%, and over 98%, respectively. The obtained spatial distribution was intersected with the spatial distribution of the total vegetable fields from Feidong's agricultural survey data in 2016 to obtain the spatial distributions of the OV and GV fields ( Figure 4). The obtained spatial distribution was intersected with the spatial distribution of the total vegetable fields from Feidong's agricultural survey data in 2016 to obtain the spatial distributions of the OV and GV fields ( Figure 4).

Soil Sampling and Analysis
The soil tests were completed at the Anhui Institute of Geological Experiment. After the samples were digested by HCl-HNO3-HClO4-HF, the concentrations of Cd, Pb, and Cr were determined by inductively coupled plasma mass spectrometry (ICP-MS 7700, Agilent Technologies, Santa Clara, CA, USA); the concentrations of Hg and As were determined by atomic fluorescence spectrometry (AFS-8220, JiTian Technologies, Beijing, China) after digestion by aqua regia; the soil pH was measured using a pH meter (SG8, METTLER TOLEDO, Zurich, Switzerland). The laboratory used national standard materials and repeated sample tests for quality monitoring. The accuracy, precision, and reported percent of the various indicators were controlled at 0.10-0.12, 10-20%, and over 98%, respectively.

Soil Sampling and Analysis
The soil tests were completed at the Anhui Institute of Geological Experiment. After the samples were digested by HCl-HNO 3 -HClO 4 -HF, the concentrations of Cd, Pb, and Cr were determined by inductively coupled plasma mass spectrometry (ICP-MS 7700, Agilent Technologies, Santa Clara, CA, USA); the concentrations of Hg and As were determined by atomic fluorescence spectrometry (AFS-8220, JiTian Technologies, Beijing, China) after digestion by aqua regia; the soil pH was measured using a pH meter (SG8, METTLER TOLEDO, Zurich, Switzerland). The laboratory used national standard materials and repeated sample tests for quality monitoring. The accuracy, precision, and reported percent of the various indicators were controlled at 0.10-0.12, 10-20%, and over 98%, respectively.

Geostatistical Interpolation and Kernel Density Simulation
The ARCMAP10.2 platform (Environmental Systems Research Institute, Redlands, CA, USA), designed by the Environmental Systems Research Institute (ESRI), was used for the geostatistical interpolation of metal(loid)s and the kernel density simulation of vegetable fields. The equations of the two techniques used were: (1) Ordinary Kriging According to the ordinary kriging interpolation technique, if the attribute value of variable Z(x) for study area a at point x i ∈ A (i = 1, 2, . . . , n) was Z(x i ), then the kriging interpolation result of the attribute value Z * (x 0 ) at the to-be-interpolated point x 0 ∈ A was the weighted sum of the attribute values Z(x i )(i = 1, 2, . . . , n) of the most neighboring known sample points, namely: where λ i is the weight assigned to each Z(x i ) value, with their sum being 1, and n is the amount of the most neighboring sampled data points used for the estimation.
(2) Kernel Density Kernel density estimation is a particularly useful method of density estimation. A set of data can be used to continuously replace discrete histograms to create a smooth curve. The universal equation of kernel density estimation is mathematically expressed as follows: where k(x) is the kernel function, which is usually a smooth symmetric function, such as a Gaussian function, and h is the smoothing bandwidth if it is greater than 0, which controls the amount of smoothing. Kernel density estimation smooths each data point X i into a small density bump and then adds all of the small bumps to obtain a final density estimate.

Statistical Analyses
Statistical analyses were performed using SPSS 19.0 (International Business Machines Corporation, Armonk, NY, USA) [28]. Pearson's correlation analysis was used to measure the linear relationship between the two variables by detecting if two phenomena (statistics) are correlated. In this study, Pearson's correlations were performed to verify if the distributions of the soil metal(loid)s were closely related to the locations of the vegetable fields.

Spatial Autocorrelation Analysis
Moran's I is the most commonly used method for calculating spatial autocorrelation, both global and local spatial autocorrelations [1,29,30].
The global Moran's I is derived from the covariance relationship between correlation coefficients. The value of the covariate also represents the correlation between the two sets of values. The global Moran's I is expressed as follows: where Wij is the spatial adjacency weight matrix of each spatial unit i and the spatial unit j (j = {1, 2, 3, . . . , n } ) in the study area (1 indicates that i is adjacent to j, while 0 indicates that i and j are not adjacent), x i is the value of each variable in a variable set, with their average being x, and y i is the value of each variable in another variable set, with their average being y.
The local Moran's I is a local measure of spatial autocorrelation, which is used to identify the locations of spatial clusters and spatial outliers. It is calculated as follows: The definition of each variable is similar to that provided for the above equation. The Moran's I value calculated according Equation (4) is between −1 and 1. If the value is greater than 0, the correlation is positive. If the value is less than 0, the correlation is negative. A larger value indicates a greater spatial distribution correlation, i.e., the spatial distribution is clustered. A smaller value indicates a small spatial distribution correlation. When the value tends to be 0, the spatial distribution is random.
The global and local Moran's I were analyzed in GeoDa [31], which is a user-friendly software with a wide set of spatial analysis methods, by which the global Moran's I value and its significance can be obtained, as well as the local spatial autocorrelation classification results of local Moran's I statistical analysis.

Technical Path
A geographic polygon database was created in ArcMap (Environmental Systems Research Institute, Inc, RedLands, CA, USA). The data for the metal(loid)s in soils, including the soil metal(loid) concentrations and latitude/longitude information, were imported into the database. The Kernel Density tool in ArcMap was used to obtain the kernel density estimate of vegetable field data in the study area. A 1 × 1 km grid was created for the entire study area. The join association method in ArcGIS was then used to correlate the metal(loid) concentrations with the kernel density of vegetable fields based on spatial locations, and the correlated grid data were finally exported. Then, a spatial weight matrix was constructed in GeoDa (GeoDa Center for Geospatial Analysis in the University of Chicago, Chicago, IL, USA), and the bivariate local Moran's I was used to generate a cluster map. The generated type results were saved to the data file for subsequent operations and analyses in ArcGIS.

Variations and Correlations of Soil Heavy Metals
The analysis of the raw sampling data (Table 2) indicated that the coefficient of variation of soil Hg was the highest (85.597%), followed by soil Cd (41.804%), and soil As (35.081%), which may be related to human activities. The variation intensity of soil pH was the lowest. The coefficient of variation can only qualitatively reflect the overall level of soil attributes and their changing trends but cannot reflect changes in their spatial characteristics. Neither can it quantitatively describe how soil attributes change with sampling sites or determine whether structural or random factors are the most influential on variation. Semivariogram features are needed to explain the above characteristics.
The nugget effect, one of the features of the semivariogram, is the ratio of the nugget value to the sill value, i.e., C0/(C0 + C), which indicates the proportion of the spatial variability due to randomness in the total variation of the system. It falls into three grades according to the spatial correlation degrees of the region-specific variables. If the nugget coefficient is <25%, the spatial correlation is strong. If the nugget coefficient is 25-75%, the spatial correlation is moderate. If the nugget effect is >75%, the spatial correlation is weak. It can be seen from Table 3 that the soil As showed a strong spatial correlation, while the other soil attributes showed a moderate spatial correlation.

Spatial Distribution of Metal(loid)s-Kriging Model Interpolation
The concentrations of the five metal(loid)s examined here (Table 2) were compared with the soil pollution risk values of agricultural lands published in China's Environmental Quality Standard for Soils [32]. The average concentrations of the five metal(loid)s did not exceed Grade II soil quality. However, the concentrations of the five metal(loid)s exceeded the soil background value in some areas, indicating that the farmland soils in the study area may have been contaminated by metal(loid)s to different degrees. Combining the results of the spatial distribution of vegetable fields with that of metal(loid)s, it can be depicted that the locations of the vegetable fields were closely related to soil pH and heavy metal spatial distributions.
For the purpose of understanding the spatial distribution regularities of the different metal(loid)s, the kriging interpolation method of the geostatistical module in ArcGIS was used to assess the metal(loid) concentrations in each type of soil in the study area. Based on this, a spatial distribution map of the heavy metal concentrations in the soil was created ( Figure 4). The spatial distributions of Cr, Pb, Cd, As, and Hg, as well as the pH of sampled areas showed an increasing or decreasing trend, and all these attributes had obvious high-value areas. All metal(loid)s showed high values in the south-eastern area, possibly due to phosphate mining. The concentrations of Cr, Pb, and Hg were generally low in the study area. While Cr was clustered in several small sites of the north-eastern area, possibly due to local paint factories, the concentration of As was high in the north-western area possibly due to building material and paint factories. High Hg concentrations were found in the south-eastern area, mainly due to mining plants, and Cd was mainly distributed in the southern area due to cities and mines, and in the north-western area due to agriculture, especially vegetable fields. Lead and Cr were mainly distributed in the western area, possibly due to agricultural production rather than industry and mining. Since long-term vegetable planting caused soil acidification, high pH values ( Figure 5) were found to be negatively correlated with the spatial distribution of vegetable fields ( Figure 4).

The Vegetable Field Kernel Density
The probability density of different vegetable field types in the study area was obtained by kernel density analysis method. The centroid of each polygon was extracted as the point data for the kernel density calculation. The default search radius (bandwidth) was calculated based on the spatial configuration and number of input points. The bandwidth was calculated by dividing the minimum value of the width or the height of the output range in the study area by 30. For example, the east-west width of the study area is 47,670 m, and the default bandwidth is 1589. This approach to calculate the default radius generally avoids the "ring around the points" phenomenon that often occurs with sparse datasets. This approach corrects for spatial outliers-input points that are very far away from the rest-so they will not make the search radius unreasonably large. The results show that both vegetable fields have spatial agglomeration characteristics, mostly around villages and towns ( Figure 6).

Correlation Analysis
The effects of the vegetable fields on the geochemical characteristics of soil metal(loid)s and pH were analyzed, but the normality test did not pass for the Kolmogorov-Smirnova test and the Shapiro-Wilk test in Table 4. Thus, the original data were transformed by the bloom method and then were assessed by examining Pearson's correlations between metal(loid) interpolation data and vegetable field kernel density in the 2276 1 × 1 km grids across the study area. The correlation coefficients are listed in Table 5. The p-values lower than 0.01 evidenced the significant correlation between soil Cr, Pb, Cd, As, Hg concentrations, and pH with the density of OV and GV fields. Since the planting times of GV fields varied greatly, there was no obvious regularity of their effects on metal(loid) concentrations. Open field vegetables showed high significant correlations with soil attributes due to their long-term planting history, and these were positive for Cd and Hg, and negative for Cr, Pb, As, and pH. This was primarily because in OV fields, the massive use of organic fertilizers resulted in bioconcentrations of Cd and Hg, and long-term fertilization significantly reduced the pH, leading to the transformation and migration of Cr, Pb, and As.

Heavy Metal Concentrations on the Different Types of Soil and Vegetable Fields
To analyze whether heavy metal accumulation in the soil was related to the soil type and vegetable field type, the ArcGIS statistical data tool (Zonal Statistics as Table) was used to superimpose the soil type, vegetable field type, and heavy metal grid data to obtain the number, mean, standard deviation, etc., of pixels in the type-specific statistically framed areas. Table 6 and Figure 7 show the average concentrations of metal(loid)s in the different types of soils. The average concentrations of all metal(loid)s were low and did not reach the soil pollution background values of Anhui. However, the heavy metal concentrations in the different types of soils varied significantly, which was related to the soil background values. In the limestone and skeletal soils, the concentrations of the four metal(loid)s other than Cd were high. In paddy soils, the concentration of Cd was high, while that of the other four metal(loid)s were low, which may be related to soil fertilization characteristics. number, mean, standard deviation, etc., of pixels in the type-specific statistically framed areas. Table 6 and Figure 7 show the average concentrations of metal(loid)s in the different types of soils. The average concentrations of all metal(loid)s were low and did not reach the soil pollution background values of Anhui. However, the heavy metal concentrations in the different types of soils varied significantly, which was related to the soil background values. In the limestone and skeletal soils, the concentrations of the four metal(loid)s other than Cd were high. In paddy soils, the concentration of Cd was high, while that of the other four metal(loid)s were low, which may be related to soil fertilization characteristics.   Table 7 and Figure 8 show the average concentrations of metal(loid)s in the different vegetable fields. The average pH and Cd concentrations of OV fields were lower than those of the GV fields, whereas the average As, Cr, Hg, and Pb concentrations of the OV fields were higher than those of the GV fields. This was possibly because the long planting history of OV   Figure 8 show the average concentrations of metal(loid)s in the different vegetable fields. The average pH and Cd concentrations of OV fields were lower than those of the GV fields, whereas the average As, Cr, Hg, and Pb concentrations of the OV fields were higher than those of the GV fields. This was possibly because the long planting history of OV fields led to soil acidification (low pH), and because OV fields are susceptible to atmospheric deposition of metal(loid)s. Therefore, the As, Cr, Hg, and Pb concentrations of the OV fields were higher. The average concentration of Cd in the GV fields was high, possibly because the amount of applied organic fertilizer in such fields far exceeded that of the OV fields and resulted in Cd accumulation in the soils. fields led to soil acidification (low pH), and because OV fields are susceptible to atmospheric deposition of metal(loid)s. Therefore, the As, Cr, Hg, and Pb concentrations of the OV fields were higher. The average concentration of Cd in the GV fields was high, possibly because the amount of applied organic fertilizer in such fields far exceeded that of the OV fields and resulted in Cd accumulation in the soils.  The pH of the GV fields increased with the planting area (Table 8), which was related to the latest compensation policy for greenhouse vegetable cultivation in Feidong, as described below. The compensation standards vary with planting areas and the larger the planting area, the higher the compensation. Compensation is given according to five levels of planting areas, i.e., 1 (0-1.3 ha), 2 (1.3-3.3 ha), 3 (3.3-5 ha), 4 (5-6.6 ha), and 5 (≥6.6 ha). Therefore, some large-scale GV fields have emerged in recent years, but their soils are not severely damaged, and their pH is high due to the short planting time. Level 1 GV fields were characterized by small planting areas and long planting times, and vegetables were mostly for the collective consumption of a family or village, so the soil pH was the lowest, posing a certain risk of migration. Level 2 and 3 GV fields used more fertilizer as they produced vegetables for sale, resulting in higher metal(loid) concentrations in the soils. However, the average metal(loid) concentrations varied with the use of different organic fertilizers. For Level 4 GV fields, the planting area was large, and the planting time was short; the use of organic fertilizers increased the concentrations of As and Cd, and reduced Pb accumulation in soils. The pH of the GV fields increased with the planting area (Table 8), which was related to the latest compensation policy for greenhouse vegetable cultivation in Feidong, as described below. The compensation standards vary with planting areas and the larger the planting area, the higher the compensation. Compensation is given according to five levels of planting areas, i.e., 1 (0-1.3 ha), 2 (1.3-3.3 ha), 3 (3.3-5 ha), 4 (5-6.6 ha), and 5 (≥6.6 ha). Therefore, some large-scale GV fields have emerged in recent years, but their soils are not severely damaged, and their pH is high due to the short planting time. Level 1 GV fields were characterized by small planting areas and long planting times, and vegetables were mostly for the collective consumption of a family or village, so the soil pH was the lowest, posing a certain risk of migration. Level 2 and 3 GV fields used more fertilizer as they produced vegetables for sale, resulting in higher metal(loid) concentrations in the soils. However, the average metal(loid) concentrations varied with the use of different organic fertilizers. For Level 4 GV fields, the planting area was large, and the planting time was short; the use of organic fertilizers increased the concentrations of As and Cd, and reduced Pb accumulation in soils. No regularity was found in OV fields, which was possibly related to the size of the investigated area and to differences in the fertilization habits. Level 1 OV fields were characterized by small planting areas and had the highest pH. Their Cd and Hg concentrations were also the highest, probably because the small vegetable fields, located on the roadsides for the convenience of vegetable growers, were susceptible to automobile exhaust and atmospheric conditions. Level 5 OV fields were large, and had the lowest pH, as well as the lowest Pb and Hg concentrations. Level 3 OV fields had the highest As concentrations and level 4 OV fields had the highest Pb and Cr concentrations (Table 9). Open vegetable fields had lower soil pH and higher soil As, Cr, and Hg concentrations than GV fields in the same area. Only OV fields of levels 2, 3, and 4 had lower Cd and Pb than GV fields in the same area. Heavy metal pollution in the OV fields was more severe than that in the GV fields. Moreover, OV fields were subject to the external environment, resulting in the loss of metal(loid)s ( Figure 9). No regularity was found in OV fields, which was possibly related to the size of the investigated area and to differences in the fertilization habits. Level 1 OV fields were characterized by small planting areas and had the highest pH. Their Cd and Hg concentrations were also the highest, probably because the small vegetable fields, located on the roadsides for the convenience of vegetable growers, were susceptible to automobile exhaust and atmospheric conditions. Level 5 OV fields were large, and had the lowest pH, as well as the lowest Pb and Hg concentrations. Level 3 OV fields had the highest As concentrations and level 4 OV fields had the highest Pb and Cr concentrations (Table 9). Open vegetable fields had lower soil pH and higher soil As, Cr, and Hg concentrations than GV fields in the same area. Only OV fields of levels 2, 3, and 4 had lower Cd and Pb than GV fields in the same area. Heavy metal pollution in the OV fields was more severe than that in the GV fields. Moreover, OV fields were subject to the external environment, resulting in the loss of metal(loid)s (Figure 9).

Spatial Autocorrelation Comparison between OV and GV
The global Moran's I describes the overall distribution of a phenomenon and determines whether this phenomenon has a spatial clustering characteristic. Generally, at the 5% significance level, if Z(I) is greater than 1.96, it indicates that the distribution of the phenomenon has significant spatial autocorrelation. If Z(I) is less than −1.96, it indicates that the distribution of the phenomenon within the study area has negative spatial autocorrelation. The global Moran's I values found here revealed that the five metal(loid)s and pH had significant spatial autocorrelations (all p < 0.01). The global Moran's I values were ranked in the following sequence: Pb > As > pH > Cr > Hg > Cd (Table 10). The bivariate local Moran's I is a measure of the spatial autocorrelation between two variables. Our results showed that the five metal(loid)s and pH had significant spatial autocorrelations (all p < 0.01), as presented in Tables 11 and 12. When the absolute bivariate local Moran's I of each soil attribute were compared, OV fields showed greater correlations than GV fields, possibly because the latter did not witness an apparent regularity in time and space. Open vegetable fields had a positive spatial correlation with soil Hg and Cd concentrations, and a negative spatial correlation with soil Pb, As, and Cr concentrations, possibly because the reduction in pH accelerated the loss of the latter three metal(loid)s after activation. The modifiable area unit problem (MAUP), a term introduced by Openshaw and Taylor's classic paper [33], has long been recognized as a potentially troublesome feature of aggregated data. The deviation between the statistical results and the spatial analysis results caused by the change of the aggregated spatial units is often described by scale effect and zoning effect. In order to clarify the MAUP effect in this paper, three scales of geographical units were selected to analyze the modifiable areal unit problem sensitivities. The results of bivariate local Moran's I value coefficients were calculated, as shown in Table 13, which was sensitive and depended on the specific geographic unit used in the study. The spatial autocorrelation coefficients of different geographic units were different, but there was a scale (1000 × 1000 m grid) that could explain the spatial sensitivity reasonably. The other two scales, the positive correlation between pH and OV was obviously inconsistent with the survey results, and the spatial relationship between other metal(loid)s and OV was difficult to explain. Therefore, the follow-up analysis only analyzed the medium-scale (1 × 1 km). Although the global Moran's I described the overall distribution of metal(loid)s and their spatial clustering, it did not specify in which regions the metal(loid)s were clustered. On the other hand, the bivariate local Moran's I distribution map produced ( Figure 10) evidenced four types of local spatial correlations of variables in each region and surrounding regions: OV fields were highly clustered and witnessed no accumulation of metal(loid)s. (4) Areas with low Moran's I for both soil metal(loid)s and OV fields (low-low area): OV fields were lowly clustered and witnessed a low accumulation of metal(loid)s. (5) As shown in Figure 10, the results of spatial autocorrelation distribution of variables were quite different due to different geographical units. The spatial autocorrelation of smaller geographical units was much stronger than that of larger geographical units. Table 13 also shows that the coefficient of bivariate local Moran's I values decreased with the increase of larger spatial units. The sensitivity analysis of the modifiable area unit problem illustrates the uncertainty caused by the scale effect, and different scales may lead to some information loss or bias.
As the bivariate local Moran's I values of the GV fields were low, this study only analyzed the bivariate local Moran's I distribution map of the OV fields ( Figure 11). Their spatial distribution characteristics showed that OV fields had different effects on the spatial distributions of the six soil attributes. Concentrations of Pb and Cr and pH were highest in the soils of high-low areas, among which the pH was the most significant. Open vegetable fields were lowly clustered, and their soil pH was high, probably because vegetable planting in OV fields removed excessive base elements, thereby leading to soil acidification. Concentrations of Pb and Cr showed a negative spatial correlation, probably because the humus resulting from organic fertilizer decomposition caused soil acidification to increase the activity of Pb and Cr in soils, thereby resulting in their loss via mobilization. Concentrations of Hg and Cd were highest in the soils of low-low areas and showed no obvious negative spatial correlations. Farmers in Feidong County used to apply organic fertilizers such as chicken manure and cow dung to grow vegetables, generally 18,750 kg/hectares for common vegetables and tomatoes and 15,000 kg/hectares for lettuce, cabbage, Hangzhou pepper, and pepper. Such organic matter had a certain reduction ability that enabled Hg and Cd in soil solutions to form sulfides and precipitate, so their migration was less important than that of Pb and Cr. The five metal(loid)s presented low values in soils of high-high areas, i.e., the cumulative effects of OV fields and metal(loid)s were not obvious, possibly because the bioconcentration of metal(loid)s in the soils was affected by many factors, such as the application of fertilizers and pesticides, agricultural planting structure, sewage irrigation, soil properties, different climatic conditions, and atmospheric deposition. Thus, OV fields alone cannot explain all the spatial variability of metal(loid) bioconcentration.  Pb and Cr. The five metal(loid)s presented low values in soils of high-high areas, i.e., the cumulative effects of OV fields and metal(loid)s were not obvious, possibly because the bioconcentration of metal(loid)s in the soils was affected by many factors, such as the application of fertilizers and pesticides, agricultural planting structure, sewage irrigation, soil properties, different climatic conditions, and atmospheric deposition. Thus, OV fields alone cannot explain all the spatial variability of metal(loid) bioconcentration. Figure 11. Bivariate local Moran's I distribution maps for metal(loid)s and pH by 1 × 1 km Unit. Figure 11. Bivariate local Moran's I distribution maps for metal(loid)s and pH by 1 × 1 km Unit.

Conclusions
By combining GIS and Moran's I spatial analysis to generate spatial autocorrelation distribution maps, the spatial clustering (positive autocorrelation) and spatial anomaly (negative autocorrelation) of vegetable fields and soil heavy metal concentrations were identified. It can be seen from GIS spatial analysis that the spatial distributions of Cr, Pb, Cd, As, and Hg content and pH in the study area showed a certain increasing or decreasing trend, and all showed spatial clustering to some extent. It can be concluded that even though the background values of metal(loid)s in the soils of Feidong are low, anthropogenic activities have contributed to an alarming increase of metal(loid) concentrations. On this basis, spatial correlation analysis was performed to assess the effects of vegetable fields on metal(loid) accumulation in soils, which is important for protecting soil quality and improving crop quality. Long-term vegetable production in the study area involves the massive use of organic fertilizers, which affects the spatial distribution of metal(loid)s in soils, resulting in the bioconcentration of Cd and Hg. Long-term fertilization significantly reduced pH, thereby causing transformation and migration of Cr, Pb, and As. The soil pH of the GV fields increased with the planting area, which was related to the vegetable compensation policy in the study area. The planting history is short, but the accumulation rate of Cd is faster, which has a certain migration risk. Open vegetable fields, with a long planting history, had lower than average soil pH and Cd concentration than GV fields, and higher average soil As, Cr, Hg, and Pb concentrations than GV fields. Metal(loid)s in soil samples caused potential hazards through non-point source pollution transfer. Therefore, special attention should be paid to the type and quantity of fertilizer in order to reduce the transformation and loss of metal(loid)s in soils.
In the future, the three-dimensional spatial variability of metal(loid)s in soil will be focused on. Through the study of metal(loid) migration and transformation in vertical and horizontal directions, we can clarify the environmental impact of metal(loid)s in OV and GV to the groundwater and surrounding ditches respectively.