Comparative Study on Different Interpolation Methods and Source Analysis of Soil Toxic Element Pollution in Cangxi County, Guangyuan City, China

: Spatial interpolation is a crucial aspect of soil toxic element pollution research, serving as a vital foundation for pollution assessment, treatment, and sustainability efforts. The selection and adjustment of interpolation methods directly influences the accuracy of spatial distribution maps and data results, thereby indirectly impacting related research. This paper conducts a comparative study of different interpolation methods and analyses the sources of soil toxic elements in the study area of Cangxi County, aiming to provide a scientific foundation for future soil management, remediation, and enhanced local sustainability. The spatial correlation of As, Cd, Hg, Mn, Pb, and Mo in 228 surface soil samples in the study area of Cangxi County is analyzed. The interpolation results, spatial distribution of OK (ordinary Kriging), IDW (inverse distance weighting), RBF (radial basis function) and the changes of pollution area after interpolation are compared. The smoothing effect is assessed based on the comparison results, interpolation accuracy, and impact on pollution assessment of OK, IDW, and RBF. The interpolation method most suitable for each metal in the study area is selected. It can be concluded that the optimal interpolation method for As, Hg, and Mn is IDW; for Cd and Mo, it is RBF; and for Pb, it is OK. After the correlation analysis of toxic elements in the soil of the study area, the PMF (positive matrix factorization) model and hotspot analysis is applied to analyzing the source of toxic elements. The analysis indicates that the predominant sources of pollution are anthropogenic, categorized into industrial activities (30.8%), atmospheric deposition caused by coal combustion and traffic exhaust (21.5%) and agricultural activities (19.5%). Natural sources, such as soil parent material, contribute to 28.2% of the pollution on average.


Introduction
Soil plays a critical role in sustainability, impacting agricultural, economic, and human health aspects.Moreover, soil pollution by toxic elements is a significant environmental concern worldwide [1], garnering considerable attention [2].The latest soil pollution survey in China revealed that 16.1% of samples surpassed the environmental quality standards set by the Ministry of Environmental Protection, with toxic elements constituting 82.4% of the pollutants [3].Pollution from toxic elements differs from organic pollution; it is nonbiodegradable [4], and can readily accumulate in the bodies of humans and animals through the food chain, posing health risks [5,6].Soil pollution by toxic elements negatively affects local sustainability, contaminating farm crops, undermining the agricultural economy, and increasing health risks for local residents.

Overview of the Study Area
The study area (31 • 51 ′ N~32 • 55 ′ N, 105 • 57 ′ E~106 • 02 ′ E) is at the junction of Sanchuan Town, Baihe Town, and Shimen Town in Cangxi County, Guangyuan City, and Sichuan Province, with an acreage of 14.88 km 2 (Figure 1).The study area is a purple hilly area, the terrain is tilted from northwest to southeast, and the slope is gentle.The study area soil is dominated by purple soil, which accounts for 70% of the total area, while paddy soil accounts for 30%.The main crop is kiwifruit.Cangxi County belongs to the subtropical humid monsoon climate zone.The average annual temperature is 16.9 • C, the average temperature in January is 6 • C, and in July it is 27 • C. The extreme minimum temperature is −4.6 • C, while the maximum temperature is 39.3 • C. The diurnal amplitude is 3~7 • C. The annual frost-free period lasts for 288 days, and the average annual rainfall is about 1100 mm.

Overview of the Study Area
The study area (31°51′ N~32°55′ N, 105°57′ E~106°02′ E) is at the junction of Sanchuan Town, Baihe Town, and Shimen Town in Cangxi County, Guangyuan City, and Sichuan Province, with an acreage of 14.88 km 2 (Figure 1).The study area is a purple hilly area, the terrain is tilted from northwest to southeast, and the slope is gentle.The study area soil is dominated by purple soil, which accounts for 70% of the total area, while paddy soil accounts for 30%.The main crop is kiwifruit.Cangxi County belongs to the subtropical humid monsoon climate zone.The average annual temperature is 16.9 °C, the average temperature in January is 6 °C, and in July it is 27 °C.The extreme minimum temperature is −4.6 °C, while the maximum temperature is 39.3 °C.The diurnal amplitude is 3~7 °C.The annual frost-free period lasts for 288 days, and the average annual rainfall is about 1100 mm.

Sample Collection and Analysis
A total of 228 topsoil (cultivated layer) samples at an average density of 16 pieces/km 2 were collected.GPS was used to navigate to the pre-set sampling point and measure its geographical coordinates.Each sampling point was marked, numbered, and documented with corresponding landscape photographs.The sampling tool used was a stainless-steel dredge shovel.The sampling medium was the soil column of 0~20 cm on the surface, and the soil column combination of 3~5 points (radial) was collected within the radius of 50 m around the sampling point as a sample.Efforts were made to minimize soil disruption to ensure that the samples were not secondary polluted during collection.The samples were stored in designated sample bags and isolated using a polyester bag.
In this study, the soil sample analysis was undertaken by the Chengdu Mineral Resources Supervision and Inspection Center of the Ministry of Land and Resources.The analysis unit strictly adheres to the "Technical standard of geological survey of China Geological Survey (DD2005-01~DD2005-03)" published by China Geological Survey in October 2005 and rigorously monitors the analysis quality of various samples through the use of standard samples, internal laboratory inspections, cryptographic inspections, and external sampling inspections.The detection center selected the full amount of soil sample elements.After being received by the detection center, air-dried, and gravel and plant debris were removed, the sample was divided into two parts.One part was ground with agate mortar, passed through a soil sieve with a diameter of 2 mm, and then ground all through a 100-mesh sieve.Different spectroscopic analysis methods were used to determine the toxic element concentration of samples.The spectroscopic analysis methods,

Sample Collection and Analysis
A total of 228 topsoil (cultivated layer) samples at an average density of 16 pieces/km 2 were collected.GPS was used to navigate to the pre-set sampling point and measure its geographical coordinates.Each sampling point was marked, numbered, and documented with corresponding landscape photographs.The sampling tool used was a stainless-steel dredge shovel.The sampling medium was the soil column of 0~20 cm on the surface, and the soil column combination of 3~5 points (radial) was collected within the radius of 50 m around the sampling point as a sample.Efforts were made to minimize soil disruption to ensure that the samples were not secondary polluted during collection.The samples were stored in designated sample bags and isolated using a polyester bag.
In this study, the soil sample analysis was undertaken by the Chengdu Mineral Resources Supervision and Inspection Center of the Ministry of Land and Resources.The analysis unit strictly adheres to the "Technical standard of geological survey of China Geological Survey (DD2005-01~DD2005-03)" published by China Geological Survey in October 2005 and rigorously monitors the analysis quality of various samples through the use of standard samples, internal laboratory inspections, cryptographic inspections, and external sampling inspections.The detection center selected the full amount of soil sample elements.After being received by the detection center, air-dried, and gravel and plant debris were removed, the sample was divided into two parts.One part was ground with agate mortar, passed through a soil sieve with a diameter of 2 mm, and then ground all through a 100-mesh sieve.Different spectroscopic analysis methods were used to determine the toxic element concentration of samples.The spectroscopic analysis methods, analytical equipment and method detection limits (MDL) of each toxic element were mea-sured are shown in Table 1.The standards for the selection of determination methods are based on the following: "Soil and sediment-Determination of mercury, arsenic, selenium, bismuth, antimony-Microwave dissolution/Atomic Fluorescence Spectrometry (HJ-680 2013)" published by Ministry of Environmental Protection (China) on 21 November 2013, for As and Hg determination using AFS; "Soil and sediment-Determination of inorganic element-Wavelength dispersive X-ray fluorescence spectrometry (HJ-780 2013)" published by Ministry of Environmental Protection (China) on 21 November 2013, for Mn and Pb determination using XRF; and "Soil and sediment-Determination of aqua regia extracts of 12 metal elements-Inductively coupled plasma mass spectrometry (HJ-803 2016)" published by Ministry of Environmental Protection (China) on 24 June 2016, for Cd and Mo determination using ICP-MS.Another part of the sample was for pH determination and direct sunlight, acid, alkali, other gases, and dust pollution were avoided.At the same time, another spare part of the sample was extracted by quartering method, the debris other than soil was removed, and the test screen with 0.84 mm aperture was used to filter.Potentiometry was used to determine pH value.

Analysis Methods
Microsoft Excel 2013 and SPSS (Statistical Package for the Social Sciences) 27.0 was used to analyze the data.ArcGIS was used to proceed geological statistical analysis.OK, IDW, and RBF interpolation were also conducted with ArcGIS.PMF was analyzed using EPA PMF 5.0.The methods and formulas employed in the data analysis are as follows: The coefficient of variation (CV) is an index that assesses the dispersion degree of a single factor in spatial distribution [28,29].The higher the CV, the stronger the interference from human activities or the greater the degree of pollution [30][31][32][33].The coefficient of variation can be categorized into weak variation (CV < 15%), moderate variation (15% ≤ CV ≤ 36%) and strong variation (CV > 36%) [34].
Kriging, widely used in geostatistics, is the most usual interpolation method [35].It is an unbiased, optimal estimation method for regionalized variables in a confined area, based on variogram theory and structural analysis [36,37].Its main advantages include not only predictive results but also error estimation and the spatial autocorrelation of soil characteristics, interpolated from known to estimated points.It exhibits excellent intrinsic correlation attributes and accuracy, beneficial for evaluating the uncertainty of prediction results [35,37].Ordinary Kriging (OK) was used in this study.
Inverse distance weighting (IDW) is an interpolation method related to spatial distance [22].It operates on the principle of similarity that the closer two objects are, the more similar they are likely to be; conversely, the greater the distance between them, the less similar they are.IDW uses the distance between the interpolation point and the sample point as a weight, assigning greater weight to sample points that are closer to the interpolation point [16].
The radial basis function method (RBF) is a variation function model.It calculates a set of weight coefficients of the nodes to be estimated by the base function and adjusts the smoothing factor in the base function to control the smoothness of the interpolation surface and the estimation accuracy, so as to achieve smooth interpolation.This method is suitable for interpolation of surfaces with gentle surface changes, which is susceptible to the influence of maximum and minimum values, and can predict the values of unknown points higher or lower than the sample points [35].
To compare interpolation methods, it is essential to minimize the original data's error, and an effective strategy for error reduction is to optimize parameters during the interpolation process for optimal results [38].Analyzing the interpolation outcomes for various toxic elements with different parameters in ArcGIS reveals that using a combination of parameters which yields the least inaccuracy is beneficial.By adjusting the range of each toxic element in the interpolation process to ensure it is sufficiently large, and by keeping the nugget coefficient (C o /(C o + C)) small, interpolation can achieve a broader spatial autocorrelation range and stronger spatial autocorrelation, thereby enhancing accuracy [39].
In cross validation, performing linear regression analysis helps achieve the algorithm's goal of unbiased mean value prediction.The measured values of sampling points were used as independent variables, and the predicted values of cross-validation were used as dependent variables.Before the intersection of the linear model and the 1:1 straight line, the predicted value is greater than the measured value, while the result is opposite after the intersection [40,41].
The commonly used evaluation indexes for cross-validation include RMSE (Root Mean Square Error) and ME (Mean Error), among others.RMSE serves as a metric for assessing the accuracy of predictions.A smaller RMSE value indicates a higher level of precision in the interpolation.ME is the measure of prediction accuracy.ME serves as a metric for assessing the impartiality of predictions, with values closer to 0 indicating a higher degree of unbiasedness in the interpolation [39].Their expressions are as follows: where V i is the error between the measured value and the predicted value.IP (Imprecision) is the variation of prediction error.The smaller the IP is, the more accurate the interpolation result is [30].Its expression is as follows: RI (the potential ecological risk index) is widely used to assess the degree of soil toxic element pollution.It introduces the toxicity coefficient of toxic elements, links the environmental ecology with the toxicological effects of toxic elements, and considers various factors for analysis [31].Its expression is as follows: where C i f is the i-th toxic element's enrichment coefficient; C i is the i-th toxic element's measured concentration; and C i n is the i-th toxic element's assessment standard, which is the i-th toxic element's soil background value of Chengdu.T i r is the i-th toxic element's toxicity coefficient.Toxicity coefficients (Tr) of the examined elements are shown in Table 2. RI is the potential ecological risk index of multi-element environment in soil; E i r is the i-th toxic element's potential ecological risk index [42,43].The potential ecological risk index E i r , RI classification and its corresponding pollution degree are shown in Table 3.The positive matrix factor analysis (PMF) is a method for quantitative analysis of multivariate factors of pollution sources by using sample composition or fingerprints through the mathematical method of receptor model [45].The operation principle of PMF model is as follows: where i is the number of samples; j is the number of elements; X ij is the j-th toxic element concentration in the i-th sample; g ik is the relative contribution of source k to sample i; f ki is the content of element j on the source; and e ij is the residual matrix.
The optimal matrices G and F are obtained by decomposing the original matrix X of the PMF model, so that the objective function Q is minimized.The calculation of the objective function Q is as follows: where U ij is the uncertainty of the content of the j-th element in the i-th sample.The model can weight each individual data point and give each data point the appropriate uncertainty.When the concentration of the element is lower than or equal to the corresponding method detection limit (MDL), the uncertainty is calculated as Formula (9), otherwise the uncertainty is calculated as Formula (10).
where δ is the relative standard deviation; and ∂ is the concentration of chemical elements [46].The method detection limit (MDL) of toxic element elements involved in this study is shown in Table 1.

Descriptive Statistics of Soil Toxic Element Content
The experimental results of this study are presented in Table 4.The soil pH value was between 5.26 and 8.56.The CV of Hg in the study area is 46.545%, suggesting a strong level of variation, which implies that Hg pollution is more likely related to human activities.The CV of Pb is 12.523%, indicating a weak variation and suggesting that Pb pollution is less likely to be related to human activities.Other toxic elements are moderately variable and related to human activities.

Assessment of Toxic Element Pollution in Soil
The results of the potential ecological risk index method are shown in Table 5.From the average value of E i r , the potential ecological risk of the surface layer in the study area, from light to moderate, is as follows: Mn (0.581) < Pb (5.382) < As (7.489) < Mo (13.848) < Hg (36.827) < Cd (51.454).According to the classification criteria in Table 3, Cd in the surface soil reaches moderate pollution, and As, Hg, Mn, Pb, and Mo are all lightly polluted in the study area.From the RI value, the soil of the study area is generally at a light pollution level.To further analyze the pollution status of each toxic element, the frequency distribution of E i r and RI was calculated and counted.The results are shown in Table 6.For all soil samples, As, Mn, Pb, and Mo are entirely at a light potential ecological risk level; For Cd, 21.93% of the sampling points are at light pollution level, 75.44% of the sampling points are at moderate pollution level, and 2.63% of the sampling points are at strong pollution level.For Hg, 58.77% of the sampling points are at light pollution level, 39.04% of the sampling sites are at moderate pollution level, and 2.19% of the sampling sites are at strong pollution level.From RI, 92.54% of the sampling points are at light pollution level, and 7.46% of the sampling points are at moderate pollution level.

Parameter Optimization of Different Interpolation Methods for Soil Toxic Elements
Through experiments, the optimal parameters of the variogram theoretical model of OK method are shown in Table 7.The optimal theoretical model for As, Mn, Pb, and Mo is exponential.Spherical is the optimal theoretical model for Cd.Gaussian is optimal for Hg.The optimal parameters of IDW and RBF are presented in Table 8.Through experiments, it is found that when the power value in the interpolation process of each element is 1, the inaccuracy is relatively minimal.In the RBF experiment, As, Hg, Mn, and Pb were interpolated using the inverse multiquadric function, while the Cd and Mo elements were interpolated using the spline with tension function.

Accuracy Comparison of Different Interpolation Methods for Soil Toxic Elements
The cross-validation method was used to assess the accuracy of the interpolation method.The regression line intersected with the 1:1 straight line, as shown in Figure 2. The correlation coefficient between the predicted and the measured values of As and Hg by IDW is the largest, followed by RBF and OK.The correlation coefficient of Cd, Pb, and Mo by OK is the smallest, significantly lower than the similar coefficients with IDW and RBF.The correlation coefficient of Mn is highest using IDW, and with RBF it is slightly lower than with OK.The correlation coefficient using OK ranges from 0.1170 to 0.3997.The correlation coefficient using IDW ranges from 0.8030 to 0.8388.The correlation coefficient by RBF ranges from 0.0992 to 0.8388.In general, the fitting degree of OK is significantly lower than that of IDW and RBF, and its correlation coefficient is also generally smaller than that of the other two methods.
The RMSE, ME, and IP values were computed to assess the accuracy of the three interpolation methods.The results are presented in Table 9    The RMSE, ME, and IP values were computed to assess the accuracy of the three interpolation methods.The results are presented in Table 9.For As interpolation, OK yields an ME closest to 0 and has the smallest IP, whereas RBF achieves the smallest RMSE.For Cd interpolation, RBF's ME is closest to 0, with OK and IDW showing similar results, and relatively low RMSE and IP values.When OK is used to interpolate Hg, it results in an ME closest to 0, with its RMSE and IP values being slightly smaller than those achieved by RBF and IDW.For Mn interpolation, IDW's ME is closest to 0, with the smallest RMSE and IP observed when using OK method.Pb interpolation conducted by OK is remarkably accurate, with ME near 0 and the lowest RMSE and IP values.OK's ME for Mo interpolation approaches 0, while RBF and IDW present similar and relatively low RMSE and IP.On average, the ME generated by IDW is relatively close to 0, and the RMSE and IP produced by OK are comparatively small.

Comparison of Interpolation Results' Statistics and Spatial Distribution of Different Interpolation Methods for Soil Toxic Elements
The maximum, minimum, mean value, SD (standard deviation) and CV (coefficient of variation) of the results by the three interpolation methods were statistically analyzed and compared with the relevant data of the sampling points.The results are presented in

Comparison of Interpolation Results' Statistics and Spatial Distribution of Different Interpolation Methods for Soil Toxic Elements
The maximum, minimum, mean value, SD (standard deviation) and CV (coefficient of variation) of the results by the three interpolation methods were statistically analyzed and compared with the relevant data of the sampling points.The results are presented in Table 10.The smoothing effect of the three interpolation methods impacts the measured concentration range of toxic elements and the intensity of the smoothing effect can be explained by the impact amplitude.In general, the smoothing effects of IDW and RBF are significantly weaker than that of OK.The concentration range of toxic elements showed different degrees of reduction after interpolation, OK had the largest reduction, followed by RBF, while IDW had the smallest reduction.The mean values of the six toxic element elements after different interpolation methods did not change significantly.The average value increased for all elements except Cd after IDW interpolation, while the average values of the other elements decreased by varying degrees after different interpolation methods were applied.Except for Cd and Mo interpolated by IDW, which are closest to the measured values, the mean values of other elements obtained by RBF are the closest to the sampling points, followed by IDW, with OK showing the largest discrepancy.After interpolation, both SD and CV significantly reduced, a result of the smoothing effect.The SD and CV of Cd and Hg are the closest to the sampling point after OK interpolation.For As, Mn, and Pb, the smallest SD and CV values occur after IDW interpolation, while those for Mo are closest to the sampling point after RBF interpolation.
The content distribution maps of different toxic element elements were generated using various spatial interpolation methods in ArcGIS, and the results are presented in Figure 3.The soil content of As is higher in the northern part of the study area and at the junction of Baihe Township and Shimen Township in the central part.When comparing the three methods for the prediction of low concentration (3.100~8.762mg•kg −1 ) of the As element, it is found that OK is rougher, while IDW provides more detail for a higher concentration range (9.657~20.200mg•kg −1 ).The Cd soil content in the study area is higher in the northern Sanchuan Town and the central and southern Baihe Township.All three interpolation methods have a "bull's eye" phenomenon near the sampling points with higher content, potentially indicating point source pollution and a small contaminated area due to concentrations near the relevant sampling points exceeding standards.The prediction of OK for higher concentrations (0.235~0.589 mg•kg −1 ) is significantly rougher, whereas the prediction of RBF for both higher and lower concentrations is relatively more accurate.The Hg content is higher at the junction of Baihe Township and Shimen Township in the central part of the study area.The results of RBF interpolation are significantly rougher than the other two methods, while OK provides better predictions than IDW.The Mn content is mainly concentrated in the soil at the junction of Baihe Township and Shimen Township in the south of the study area, in Baihe Township to the east, and in Sanchuan Township to the north.For predicting Mn, IDW performs significantly better than the other two methods for both higher and lower concentration ranges.The soil Pb content at the junction of Baihe Township and Shimen Township in the middle of the study area and Baihe Township in the north is relatively high.For the prediction of Pb, the three methods yield similar results, with the OK method slightly outperforming the other two from a detailed perspective.The soil content of Mo is mainly concentrated along the border between Baihe Township and Shimen Township in the middle of the study area and in Baihe Township to the north.Although the prediction of OK for lower concentration (0.272~0.467 mg•kg −1 ) is more detailed, compared to the other two methods, its prediction for medium and high concentration (0.534~1.240 mg•kg −1 ) is less accurate, with RBF being the best for predicting medium and high concentrations.
the northern Sanchuan Town and the central and southern Baihe Township.All three interpolation methods have a "bull's eye" phenomenon near the sampling points with higher content, potentially indicating point source pollution and a small contaminated area due to concentrations near the relevant sampling points exceeding standards.The prediction of OK for higher concentrations (0.235~0.589 mg•kg −1 ) is significantly rougher, whereas the prediction of RBF for both higher and lower concentrations is relatively more accurate.The Hg content is higher at the junction of Baihe Township and Shimen Township in the central part of the study area.The results of RBF interpolation are significantly rougher than the other two methods, while OK provides better predictions than IDW.The Mn content is mainly concentrated in the soil at the junction of Baihe Township and Shimen Township in the south of the study area, in Baihe Township to the east, and in Sanchuan Township to the north.For predicting Mn, IDW performs significantly better than the other two methods for both higher and lower concentration ranges.The soil Pb content at the junction of Baihe Township and Shimen Township in the middle of the study area and Baihe Township in the north is relatively high.For the prediction of Pb, the three methods yield similar results, with the OK method slightly outperforming the other two from a detailed perspective.The soil content of Mo is mainly concentrated along the border between Baihe Township and Shimen Township in the middle of the study area and in Baihe Township to the north.Although the prediction of OK for lower concentration (0.272~0.467 mg•kg −1 ) is more detailed, compared to the other two methods, its prediction for medium and high concentration (0.534~1.240 mg•kg −1 ) is less accurate, with RBF being the best for predicting medium and high concentrations.

Effects of Different Interpolation Methods on the Assessment Results of Soil Toxic Element Pollution
The data obtained through various interpolation methods were evaluated using the potential ecological risk index method, and the proportion of polluted area was calculated.Concurrently, these results were compared and analyzed alongside those in Table 6.The results are presented in Table 11, and the distribution map of polluted area formed in ArcGIS is shown in Figure 4.All sampling points for As, Mn, Pb, and Mo elements are classified at a light pollution level, and the three interpolation methods do not significantly influence the pollution assessment of these elements.
For Cd, the three interpolation methods all increased the acreage of moderate pollution area (40 ≤  < 80) and reduced the proportion of light pollution area ( < 40) and strong pollution area (80 ≤  < 160).The order of decrement in the proportion of light pollution is IDW (21.103%) > RBF (17.255%) > OK (12.456%).It can be observed in the spatial distribution map that the OK method eliminates the area of strong pollution due to its smoothing effect, which significantly impacts the pollution assessment.For the acreage of light pollution areas, OK is comparable to RBF and the predicted pollution areas interpolated by the two methods is larger than IDW.

Correlation and Source Analysis of Toxic Elements in Soil
The results of the correlation analysis, including the correlation coefficient matrix and the Pearson correlation coefficient for toxic elements, are presented in Figure 5.The results showed that the correlations between As and Mo (r = 0.66), As and Pb (r = 0.The positive matrix factorization (PMF) model was used to analyze the pollution sources of soil toxic elements in the study area.The PMF's input data comprised concentration data and corresponding uncertainty data.The studied elements were classified as "strong" based on their signal-to-noise ratio (S/N).Upon inputting four source factors into the program, the Q value stabilized, indicating the model's applicability [50].Concurrently, the fitting degree (R 2 ) for As, Cd, Hg, Mn, Pb, and Mo reached 0.862,0.700,0.924, 0.977, 0.997, and 0.632, respectively, indicating its applicability and reliability of the results [50].The contribution rates and average contribution rates of each factor are shown in Table 12 and Figure 6.It can be found that the major toxic elements contributed by factor 1 are Cd and Hg, with contribution rates of 45.8% and 71.9%, respectively.The main toxic elements contributed by factor 2 are Mn and As, with contribution rates reaching 49.9% and 37.6%.The major toxic elements contributed by factor 3 are As, Pb, and Mo, with contribution rates of 50.6%, 40.0%, and 50.5%, respectively.The major toxic element contributed by factor 4 is Cd, with a contribution rate of 54.2%.Cd < 160).The order of decrement in the proportion of light pollution is IDW (21.103%) > RBF (17.255%) > OK (12.456%).It can be observed in the spatial distribution map that the OK method eliminates the area of strong pollution due to its smoothing effect, which significantly impacts the pollution assessment.For the acreage of light pollution areas, OK is comparable to RBF and the predicted pollution areas interpolated by the two methods is larger than IDW.
For Hg, all three methods increased the proportion of light pollution area (E i Hg < 40) and reduced the proportion of moderate pollution area (40 ≤ E i Hg < 80) and strong pollution area (80 ≤ E i Hg < 160).The order of amplification for light pollution is OK (20.936%) > IDW (19.104%) > RBF (18.286%), and the order of decrement for moderate pollution was OK (18.746%) > IDW (16.922%) > RBF (16.089%).OK also does not present the prediction of strong pollution area.For the light and moderate pollution area, the predictions of the three methods are similar.
In general, all three methods exhibit the phenomenon of increasing the acreage of the pollution area with larger proportion, and the increases in the three methods are similar.This phenomenon of data concentration and accuracy decline, caused by the decrease in data discreteness, is also a manifestation of the smoothing effect.The greater the increase or decrease, the more pronounced the smoothing effect.

Correlation and Source Analysis of Toxic Elements in Soil
The results of the correlation analysis, including the correlation coefficient matrix and the Pearson correlation coefficient for toxic elements, are presented in Figure 5 The positive matrix factorization (PMF) model was used to analyze the pollution sources of soil toxic elements in the study area.The PMF's input data comprised concentration data and corresponding uncertainty data.The studied elements were classified as "strong" based on their signal-to-noise ratio (S/N).Upon inputting four source factors into the program, the Q value stabilized, indicating the model's applicability [50].Concurrently, the fitting degree (R 2 ) for As, Cd, Hg, Mn, Pb, and Mo reached 0.862,0.700,0.924, 0.977, 0.997, and 0.632, respectively, indicating its applicability and reliability of the results [50].The contribution rates and average contribution rates of each factor are shown in Table 12 and Figure 6.It can be found that the major toxic elements contributed by factor 1 are Cd and Hg, with contribution rates of 45.8% and 71.9%, respectively.The main toxic elements contributed by factor 2 are Mn and As, with contribution rates reaching 49.9% and 37.6%.The major toxic elements contributed by factor 3 are As, Pb, and Mo, with contribution rates of 50.6%, 40.0%, and 50.5%, respectively.The major toxic element contributed by factor 4 is Cd, with a contribution rate of 54.2%.

Correlation and Source Analysis of Toxic Elements in Soil
The results of the correlation analysis, including the correlation coefficient matrix and the Pearson correlation coefficient for toxic elements, are presented in Figure 5.The results showed that the correlations between As and Mo (r = 0.66), As and Pb (r = 0.54), Mo and Pb (r = 0.54) are high.The correlations between As and Hg (r = 0.21), As and Mn (r = 0.015), Cd and Hg (r = 0.28), Cd and Mn (r = 0.13), Cd and Pb (r = 0.078), Hg and Mo (r = 0.44), Hg and Pb (r = 0.43) are general.There is a negative correlation between As and Cd (r = −0.25),Cd and Mo (r = −0.23),Hg and Mn (r = −0.16),Mn and Mo (r = −0.15),Mn and Pb (r = −0.23).The positive matrix factorization (PMF) model was used to analyze the pollution sources of soil toxic elements in the study area.The PMF's input data comprised concentration data and corresponding uncertainty data.The studied elements were classified as "strong" based on their signal-to-noise ratio (S/N).Upon inputting four source factors into the program, the Q value stabilized, indicating the model's applicability [50].Concurrently, the fitting degree (R 2 ) for As, Cd, Hg, Mn, Pb, and Mo reached 0.862,0.700,0.924, 0.977, 0.997, and 0.632, respectively, indicating its applicability and reliability of the results [50].The contribution rates and average contribution rates of each factor are shown in Table 12 and Figure 6.It can be found that the major toxic elements contributed by factor 1 are Cd and Hg, with contribution rates of 45.8% and 71.9%, respectively.The main toxic elements contributed by factor 2 are Mn and As, with contribution rates reaching 49.9% and 37.6%.The major toxic elements contributed by factor 3 are As, Pb, and Mo, with contribution rates of 50.6%, 40.0%, and 50.5%, respectively.The major toxic element contributed by factor 4 is Cd, with a contribution rate of 54.2%.

Spatial Distribution of Hotspots of Toxic Elements in Soil
Based on the source analysis of the soil toxic elements, the contribution rate o source factor was calculated to further explore the relationship between the main factor of each toxic element and its spatial distribution.Using the contribution val the main source factors from the PMF model for the toxic elements of 228 sampling p the spatial distribution of hotspots of soil toxic elements was analyzed.The resu presented in Figure 7.It can be observed that the spatial distributions of As, Mo, a hotspots are very similar, suggesting a possible common source for these three elem The high hotspots (99%, 95%, and 90% confidence intervals) are concentrated in th struction land and cultivated land in the central and northern parts of the study are

Spatial Distribution of Hotspots of Toxic Elements in Soil
Based on the source analysis of the soil toxic elements, the contribution rate of each source factor was calculated to further explore the relationship between the main source factor of each toxic element and its spatial distribution.Using the contribution values of the main source factors from the PMF model for the toxic elements of 228 sampling points, the spatial distribution of hotspots of soil toxic elements was analyzed.The results are presented in Figure 7.It can be observed that the spatial distributions of As, Mo, and Pb hotspots are very similar, suggesting a possible common source for these three elements.The high hotspots (99%, 95%, and 90% confidence intervals) are concentrated in the construction land and cultivated land in the central and northern parts of the study area (Figure 7a,e,f).The high hotspots (99%, 95%, and 90% confidence intervals) of Cd are mainly concentrated in cultivated land and forest land in the eastern and southern parts of the study area (Figure 7b).The high hotspots (99%, 95%, and 90% confidence intervals) of Hg are mainly concentrated in the construction land in the middle of the study area, and a small number of high hotspots are distributed in the northeast and northwest of the study area (Figure 7c).The high hotspots (99%, 95%, and 90% confidence intervals) of Mn are less, mainly distributed in woodland and cultivated land in the northern part of the study area (Figure

Comparative Analysis of Toxic Element Concentration and Pollution Status in Soil
Considering the Chinese and Chengdu soil background values as benchmarks, the soil toxic element pollution in the study area was relatively moderate.Only the mean concentration of Cd simultaneously exceeded the Chinese soil background value by 2.299

Comparative Analysis of Toxic Element Concentration and Pollution Status in Soil
Considering the Chinese and Chengdu soil background values as benchmarks, the soil toxic element pollution in the study area was relatively moderate.Only the mean concentration of Cd simultaneously exceeded the Chinese soil background value by 2.299 times and the Chengdu soil background value by 1.715 times, while the mean concentration of Pb exceeded 7.635% of the Chengdu soil background value.Compared to other parts of China (as shown in Table 13), our study area has a relatively high concentration of As.From the perspective of pollution assessment, the mean E i r of As is relatively low and does not reach the level of moderate pollution (Table 3).The mean concentration of Cd in this study area is higher than that of many other study areas, with an E i r higher than those recorded in Lanzhou [51], and Shuozhou [52], suggesting a need for long-term local monitoring.Compared to other selected study areas, the mean concentration of Hg in this study area is relatively high.Although the average E i r dose not reach the moderate pollution level (Table 3), it still needs long-term monitoring.The average soil concentration and E i r of Mn in this study area are similar to those in other study areas [51,53,54], indicating almost no ecological risk (Table 3).The average concentration of Pb in this study area is relatively high, with the exception of the Luoyang study area [55], which records a significantly higher average concentration.However, its average potential ecological risk index is similar to those of other study areas, suggesting a mild pollution level (Table 3).The data on Mo are relatively limited, with this study area recording a low concentration indicative of light pollution (Table 3).In general, the RI in this study area is lower than that in many other study areas, suggesting a light pollution level (Table 3).This implies that the soil toxic element pollution in this study area is not severe.

The Comparison of Spatial Interpolation Methods
Diverse impacts on the spatial dispersion of soil toxic elements are produced by different interpolation methods.As shown in Table 9, with optimal parameters, the interpolation accuracy of OK, IDW, and RBF for each toxic element can be observed.Generally, the ME, RMSE, and IP values of the same toxic element after interpolation using different methods are similar, whereas the accuracy of the same interpolation method can vary greatly with different toxic elements.Overall, the accuracy of OK is relatively high.For the prediction of As, Mn, and Pb, OK has obvious advantages in accuracy, because it produces the smallest IP.While the ME, RMSE, and IP values for Cd and Hg obtained using the OK, IDW, and RBF methods are similar, the IDW method exhibits relatively higher accuracy despite the slight differences.The interpolation of Mo using the RBF method results in an ME value closer to 0 and smaller RMSE and IP values, making it a relatively optimal choice.
The statistics of the interpolation results across the three methods reveal a common trend: the range, standard deviation (SD), and coefficient of variation (CV) decrease, while the mean value changes compared to the original sampling point data.This can be attributed to the inherent smoothing effect of these interpolation methods.The influence of smoothing effect on interpolation results can be evaluated by the variation range of each value.As indicated in Table 10, the OK method has the most pronounced smoothing effect, with IDW and RBF showing similar effects.This may be because both IDW and RBF methods are types of deterministic interpolation and aim to retain the actual value of the sampling point data [18].The three interpolation methods exhibit distinct spatial distribution effects, as evidenced by the generated spatial distribution map.As shown in Figure 3, the OK method provides the best prediction for spatial distribution of Hg, with a more detailed depiction of varying concentration ranges.The spatial distribution for As, Mn, and Pb predicted by IDW is the best.RBF has the best spatial distribution for the prediction of Cd and Mo.Studies have shown that OK can show a significant smoothing effect in areas with great changes in element content and poor spatial autocorrelation [58].Therefore, the smoothing effect of OK on the spatial distribution of As, Cd, Mn, Pb, and Mo in this study is stronger than that of the other two methods, while OK can predict the spatial distribution of Hg in detail, which may be because of the high Moran's index and high spatial autocorrelation of As, Cd, Mn, Pb, and Mo (Table 14) [59].The variation system (46.545%) of Hg sampling points was higher, which means that the content of Hg in this study changed greatly.At the same time, the nugget coefficient (C o /(C o + C)) was higher (2.580) and Moran's index (0.056) was lower, indicating that the spatial autocorrelation of Hg was poor (Table 14) [39,59].Given that the pollution assessment of most soil toxic elements indicates a light pollution level, this study focuses on Cd and Hg, two elements with a higher pollution degree, to illustrate the impact of OK, IDW, and RBF on the pollution assessment results.The smoothing effect of OK results in the disappearance of strongly polluted areas with the smallest acreage, leading to inaccurate predictions of more severe pollution zones, which negatively impacts local environmental protection and management.Ultimately, RBF has a minimal impact on the pollution assessment of Cd, and similarly, IDW's influence on the assessment of soil Hg pollution is negligible.It is apparent that OK, despite having lower RMSE and IP values, does not predict polluted areas more accurately than the other two methods.This suggests that RMSE and IP, as measures of prediction accuracy, fail to capture the estimation error of local extremes.Furthermore, this study finds that a lower IP value corresponds to a more pronounced smoothing effect.These findings align with the results of a study conducted by Xie et al. (2011) [60].

Analysis of Influencing Factors and Sources
In this study, the PMF model was applied to analyze the sources of toxic element pollution.The analytical findings are presented in Figure 6.The primary sources of pollution comprise four main factors.The average contribution rates of factors 1, 2, 3, and 4 are 30.8%,21.5%, 28.2%, and 19.5%, respectively.
According to Table 12, the main toxic elements in factor 1 are Cd (45.8%) and Hg (71.9%).Studies have shown that industrial activities such as mining, smelting, and urban solid waste are the main sources of Cd in soil [61][62][63][64].Hg pollutants in soil in China are mainly derived from man-made industrial activities such as non-ferrous metal smelting, coal combustion, coal mining, and slag activities [65,66].According to the spatial distribution of Hg hotspots under the influence of factor 1 (Figure 7c), the high hotspot values are mainly concentrated in construction land, which mainly includes rural settlement and industrial and mining land.Simultaneously, considering the significant changes in Hg content and its poor spatial correlation, factor 1 could potentially represent the source from industrial activities.
The main toxic elements in factor 2 are Mn (49.9%) and As (37.6%).However, the Pearson correlation between As and Mn is not high, which suggests As and Mn may not have the same source.The concentration values of Mn in most sampling points are lower than the soil background value of Chengdu.Simultaneously, Mn, being a loading element, is likely to be derived from soil parent material [50,[67][68][69].Considering the spatial distribution of Mn hotspots influenced by factor 2 (Figure 7d), the high hotspots of Mn are sporadic and far away from the construction land.Additionally, the Moran's index of Mn is relatively low (0.065), which means the spatial autocorrelation is poor.The source may be from the natural environment, thus, factor 2 may represent the source of soil parent material.
The main toxic elements in factor 3 are As (50.6%), Pb (40.0%), and Mo (50.5%).In the correlation analysis, the Pearson correlation coefficient of these three toxic elements is also high, which may indicate that As, Mo, and Pb have a similar pollution level or are released from the same pollution source [70,71].Additionally, they are mainly affected by the same factor, which also confirms the possibility for same source of As, Pb, and Mo.According to the spatial distribution of As, Pb, and Mo hotspots under the influence of factor 3 (Figure 7a,e,f), the spatial distribution of the hotspots of the three elements is similar, and the high hotspot values are concentrated in the construction land of the study area.Meanwhile, As and Pb pollution in soil may come from coal combustion and atmospheric deposition of industrial by-products [63,72].Coal combustion in the rural settlement, incorporated into the construction land, may be the primary source of As and Pb in the soil of the study area.Additionally, one of the main sources of soil Pb and Mo pollution is traffic exhaust [73,74].The high hotspot area of Pb and Mo contains several crucial local traffic roads, suggesting that soil Pb and Mo may result from the deposition of traffic exhaust.In summary, factor 3 may represent the source of atmospheric deposition caused by coal burning and automobile exhaust.
The main toxic element in factor 4 is Cd (54.2%).In addition to industrial activities, Cd is a significant element in pesticides, phosphate fertilizers, and animal manure [63,75,76].Simultaneously, according to the spatial distribution of Cd hotspots under the influence of factor 4 (Figure 7b), high hotspot values are concentrated in woodland and cultivated land near the construction land.The main source of pollution may include pesticides and fertilizers.Therefore, factor 4 may represent the source of agricultural activity.

Conclusions
In this study, the toxic element content of soil surface was measured at 228 sampling points in Cangxi County, Guangyuan City.According to the measurement results, three different spatial interpolation methods were used.The results revealed significant errors following the use of the three types of interpolation.Compared with the measured values, the overall interpolation error was significant, but the difference between OK, IDW, and RBF for the same toxic element was relatively minor.These results may be associated with human activities and changes in environmental spatial distribution.Comparing RMSE, ME, and IP of OK, IDW, and RBF, it can be found that OK is more accurate in interpolation, followed by RBF and IDW.Simultaneously, the three interpolation methods have varying degrees of smoothing effects that influence the prediction of the maximum and minimum values, spatial distribution maps and the pollution assessment.In general, the smoothing effect of OK is more convenient than the other two methods, while IDW and RBF are similar.OK has the most significant impact on the prediction of maximum and minimum values, and the impact of IDW is slightly more pronounced than that of RBF.The three interpolation methods indicate similar spatial distribution characteristics.However, due to the differences in predicted values at the same position, the shape and acreage of different content intervals vary.Regarding the impact on pollution assessment, all three methods show a trend of increasing the acreage of the pollution area proportionately, with each method causing a different rate of increase.Based on the accuracy of interpolation, the interpolation results, the spatial distribution maps and the impact on pollution assessment, the most suitable method for each toxic element in soil can be determined.The optimal interpolation method for As, Hg, and Mn is IDW; for Cd and Mo, it is RBF; and for Pb, it is OK.Therefore, it is necessary to select the most suitable interpolation parameters and method according to the environmental factors such as soil characteristics and human factors in the study area.
In analyzing the source of soil toxic elements in the study area, we used correlation analysis, PMF model, and hotspot analysis.The analysis findings show that the main source is human activities, further subdivided into industrial activities (30.8%), potentially a significant source of soil Cd and Hg pollution, atmospheric deposition caused by coal burning and traffic exhaust (21.5%) which may be a major source of soil As, Pb, and Mo pollution and agricultural activities (19.5%) which may be a major source of soil Cd pollution.However, the soil parent material, as a natural source, which may be a major source of soil Mn pollution, only contributes 28.2% of the average contribution rate.In the future, due to the strong impact of human activities, there may be a continuous increase in the concentration of soil toxic elements in the study area.The goal of this study is to provide a scientific foundation for local soil remediation, environmental management, and conservation, aiming to enhance local sustainability.

Figure 1 .
Figure 1.Study area location and surface soil sampling point diagram.

Figure 1 .
Figure 1.Study area location and surface soil sampling point diagram.

Figure 2 .
Figure 2. Cross-validation of three interpolation methods for soil toxic elements.

Figure 2 .
Figure 2. Cross-validation of three interpolation methods for soil toxic elements.

Figure 3 .
Figure 3.Comparison of spatial distribution of soil toxic elements after interpolation by three interpolation methods.

Figure 3 .
Figure 3.Comparison of spatial distribution of soil toxic elements after interpolation by three interpolation methods.

Figure 4 .
Figure 4. Distribution of potential ecological risk pollution areas of toxic elements in soil after interpolation by three interpolation methods.

Figure 5 .
Figure 5. Correlation analysis of soil toxic elements.

Figure 4 .
Figure 4. Distribution of potential ecological risk pollution areas of toxic elements in soil after interpolation by three interpolation methods.For Cd, the three interpolation methods all increased the acreage of moderate pollution area (40 ≤ E i Cd < 80) and reduced the proportion of light pollution area (E i Cd < 40) and strong pollution area (80 ≤ E i Cd < 160).The order of decrement in the proportion of light pollution is IDW (21.103%) > RBF (17.255%) > OK (12.456%).It can be observed in the spatial distribution map that the OK method eliminates the area of strong pollution due to its smoothing effect, which significantly impacts the pollution assessment.For the acreage of light pollution areas, OK is comparable to RBF and the predicted pollution areas interpolated by the two methods is larger than IDW.For Hg, all three methods increased the proportion of light pollution area (E i Hg < 40) and reduced the proportion of moderate pollution area (40 ≤ E i Hg < 80) and strong pollution area (80 ≤ E i Hg < 160).The order of amplification for light pollution is OK (20.936%) > IDW (19.104%) > RBF (18.286%), and the order of decrement for moderate pollution was OK (18.746%) > IDW (16.922%) > RBF (16.089%).OK also does not present the prediction of strong pollution area.For the light and moderate pollution area, the predictions of the three methods are similar.In general, all three methods exhibit the phenomenon of increasing the acreage of the pollution area with larger proportion, and the increases in the three methods are similar.This phenomenon of data concentration and accuracy decline, caused by the decrease in data discreteness, is also a manifestation of the smoothing effect.The greater the increase or decrease, the more pronounced the smoothing effect.

Figure 5 .
Figure 5. Correlation analysis of soil toxic elements.

Figure 5 .
Figure 5. Correlation analysis of soil toxic elements.

Figure 6 .
Figure 6.Contribution rates and average contribution rates of PMF source apportionment fa soil toxic elements in the study area.

Figure 6 .
Figure 6.Contribution rates and average contribution rates of PMF source apportionment factors of soil toxic elements in the study area.

Figure 7 .
Figure 7. (a) The spatial distribution of As hotspots under the influence of factor 3; (b) the spatial distribution of Cd hotspots under the influence of factor 4; (c) the spatial distribution of Hg hotspots under the influence of factor 1; (d) the spatial distribution of Mn hotspots under the influence of factor 2; (e) the spatial distribution of Mo hotspots under the influence of factor 3; and (f) the spatial distribution of Pb hotspots under the influence of factor 3.

Figure 7 .
Figure 7. (a) The spatial distribution of As hotspots under the influence of factor 3; (b) the spatial distribution of Cd hotspots under the influence of factor 4; (c) the spatial distribution of Hg hotspots under the influence of factor 1; (d) the spatial distribution of Mn hotspots under the influence of factor 2; (e) the spatial distribution of Mo hotspots under the influence of factor 3; and (f) the spatial distribution of Pb hotspots under the influence of factor 3.

Table 1 .
Spectroscopic analysis used to determine the total content of each toxic element in soil and their respective method detection limits (MDL).

Table 4 .
Statistical results of soil toxic elements contents and reference soil background values for China and Chengdu area.
[49] data comes from "Soil environmental quality Risk control standard for soil contamination of agricultural land (GB 15618-2018)[49]" published by the Ministry of Ecology and Environment (China) on 22 June 2018.

Table 5 .
Statistical analysis of soil toxic element content.

Table 6 .
Frequency distribution of E i r and RI in soil.

Table 7 .
The theoretical model of variation function of soil toxic element elements and its related parameters.

Table 8 .
The optimal parameters of IDW and RBF for soil toxic element elements.
. For As interpolation, OK yields an ME closest to 0 and has the smallest IP, whereas RBF achieves the smallest RMSE.For Cd interpolation, RBF's ME is closest to 0, with OK and IDW showing similar results, and relatively low RMSE and IP values.When OK is used to interpolate Hg, it results in an ME closest to 0, with its RMSE and IP values being slightly smaller than those achieved by RBF and IDW.For Mn interpolation, IDW's ME is closest to 0, with the smallest RMSE and IP observed when using OK method.Pb interpolation conducted by OK is remarkably accurate, with ME near 0 and the lowest RMSE and IP values.OK's ME for Mo interpolation approaches 0, while RBF and IDW present similar and relatively low RMSE and IP.On average, the ME generated by IDW is relatively close to 0, and the RMSE and IP produced by OK are comparatively small.

Table 9 .
RMSE, ME, and IP of three interpolation methods for soil toxic elements.

Table 9 .
RMSE, ME, and IP of three interpolation methods for soil toxic elements.

Table 10 .
Statistics of interpolation results of different interpolation methods for toxic elements in soil.

Table 12 .
Contribution rates of soil toxic element pollution factors in the study area.

Table 12 .
Contribution rates of soil toxic element pollution factors in the study area.

Table 13 .
Comparison of toxic element concentration and pollution status in the study area of Cangxi and other research reports.

Table 14 .
Spatial autocorrelation parameters of soil toxic elements.