Soil heavy metal high-precision prediction method fusing optimal scale screening and spatial layering

Abstract

The invention relates to the technical field of geographic information, and discloses a soil heavy metal high-precision prediction method fusing optimal scale screening and spatial layering. The method aims at solving the problems that a traditional global model is insufficient in prediction precision and unclear in mechanism interpretation in spatial heterogeneity regions such as a plain urban and rural ecotone. According to the technical scheme, the method comprises the steps that firstly, a system collects multi-source environment variables and constructs a multi-scale buffer area candidate variable set, the optimal action scale of each variable is screened out in combination with soil sample data screened through geochemical background values, and an optimal scale variable set is formed; then, driving mode space layering is carried out on the research area through unsupervised clustering based on the variable set, and the heterogeneous area is divided into a plurality of relatively homogeneous sub-areas; and finally, constructing an independent random forest model for each sub-region for prediction and fusion mapping. According to the method, the precision and mechanism interpretation capability of soil heavy metal space prediction are remarkably improved, and a reliable basis can be provided for soil environment risk management and control.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geographic information, in particular to a high-precision prediction method for soil heavy metals fusing optimal scale screening and spatial stratification. BACKGROUND

[0002] The plain urban-rural ecotone is a sensitive area of environmental pollution, and the sources of soil heavy metal pollution in the area are complex and have strong spatial heterogeneity. Accurate understanding of the spatial distribution is the premise of ecological risk assessment and precise management.

[0003] At present, the commonly used spatial prediction methods such as geostatistical interpolation method (such as Kriging method) or global machine learning model (such as random forest) have obvious limitations in such areas. First, these methods often have difficulty in systematically integrating multi-scale artificial pollution source factors (such as the spatial influence of different industry enterprises) and natural factors. Second, the abnormal values caused by high values or sampling errors contained in the model training data caused by geological background will interfere with the model to capture the overall pattern of the region. Most importantly, the traditional method regards the research area controlled by multiple pollution driving modes (such as industrial driving and agricultural driving) and the internal heterogeneity as a homogeneous whole for modeling, ignores the differences in physical and chemical mechanisms of the pollution process, and leads to a bottleneck in prediction accuracy and unclear mechanism of the prediction results.

[0004] Therefore, there is an urgent need for a high-precision prediction method that can finely depict and adapt to spatial heterogeneity and has clear mechanism interpretation ability. SUMMARY

[0005] In view of the shortcomings of the prior art, the technical scheme adopted by the present application to solve the technical problems is: a high-precision prediction method for soil heavy metals fusing optimal scale screening and spatial stratification, comprising the following steps: S1: Multi-source auxiliary variable collection and multi-scale buffer zone construction, collecting human source variables, natural flow variables and natural sink variables related to soil heavy metal enrichment, and calculating the values of at least one distance type or density type human source variable at multiple predetermined buffer zone scales to form a candidate variable set; S2: Representative soil sample collection and chemical analysis, laying sampling points in the target area, collecting soil samples and measuring the content of the target heavy metal; S3: Training data screening based on geochemical background value, calculating the geochemical background value and abnormal threshold value of the soil heavy metal in the study area, removing the sampling points exceeding the abnormal threshold value, and constructing an optimized training data set; S4: Optimal scale screening and driving mode spatial stratification: S4.1: Global optimal scale screening, using the optimized training data set, a global optimal scale is selected for each variable from the candidate variable set through variable importance evaluation or correlation analysis, forming an optimal scale variable set; S4.2: Driving mode space stratification, based on the optimal scale variable set, the entire study area is divided into several sub-regions with relatively homogeneous pollution driving mechanisms; S5: Partitioning to build random forest prediction model, for each sub-region divided in S4.2, an independent random forest regression model is built, and the heavy metal content of the sampling points in the sub-region is used as the dependent variable, and the optimal scale variable set is used as the independent variable for model training; S6: Research area heavy metal spatial distribution mapping, the spatial grid point data of the entire study area is input into the corresponding partition random forest model according to its belonging sub-region, the heavy metal content of each point is predicted, and finally the high-precision spatial distribution map of the whole region is generated.

[0006] Preferably, the S4.2 driving mode space stratification is specifically: directly using the optimal scale variable set to perform unsupervised clustering analysis on all spatial units in the study area, and defining each class obtained by clustering as a pollution driving mode region.

[0007] Preferably, the S4.2 driving mode space stratification is specifically: first, spatial autocorrelation analysis is performed on the heavy metal content data in the optimized training data set to identify spatial aggregation areas; then, the environmental variable values of the optimal scale variable set in these spatial aggregation areas are statistically analyzed to define their dominant pollution driving modes, thereby completing the spatial stratification.

[0008] Preferably, the human source variable in S1 includes the number density and enterprise area of industrial enterprises of different industry types obtained based on GIS and POI data; the natural flow variable includes terrain variables (such as elevation, slope), hydrological variables (such as distance from main rivers and water bodies), and meteorological condition variables (such as wind speed and direction); the natural sink variable includes soil property variables such as clay content, silt content, sand content, pH value, and organic matter content.

[0009] Preferably, the predetermined buffer zone scale in S1 is specifically 400m, 600m, 800m, 1000m, 1200m, 1400m, 1600m, 1800m, 2000m, 2500m, 3000m, 4000m, a total of 12 annular buffer zones with different radii.

[0010] Preferably, in S3, the mean ± 3 standard deviation method or the cumulative frequency method is used to determine the geochemical background value and the anomaly threshold.

[0011] Preferably, in S4.1, the random forest model is used to evaluate the importance of variables to determine the global optimal scale of each variable. The purpose of this step is to perform feature selection, not concentration prediction. The importance score of each candidate variable is calculated by the random forest model, and for different scale variables of the same environmental factor, the scale with the highest importance score is selected as the optimal scale. This ensures that the variables used in subsequent modeling are at their most explanatory spatial scale.

[0012] Further comprising the following step S7: iterative updating mechanism of the model, when new soil sample data or environmental auxiliary variable data is obtained, steps S3 to S5 are repeatedly executed to update the optimal scale variable set and spatial stratification result, and the partition prediction model is retrained.

[0013] The beneficial effects of the present application are as follows: Through the dual driving of "scale optimization" and "mechanism partitioning", the model is more in line with the actual pollution process, which can significantly improve the prediction accuracy. The final partitioning results and variable importance of the entire process can be directly used to interpret the dominant causes of pollution, breaking the "black box" limitation of traditional machine learning models. Through training data screening and partition modeling, the interference of outliers and regional heterogeneity is effectively reduced. The framework is clear and has strong universality. The generated high-precision distribution map and driving mode partition map can provide direct and reliable decision-making basis for environmental management. BRIEF DESCRIPTION OF DRAWINGS

[0014] Fig. 1 is a flowchart of the present application; Fig. 2 is a basic statistical characteristic table of soil properties in the study area of the present application. DETAILED DESCRIPTION

[0015] The present application will be further described in detail below in conjunction with the drawings and specific embodiments. The embodiments of the present application are given for the purpose of illustration and description, and are not exhaustive or limit the present application to the disclosed forms. Many modifications and variations will be apparent to those of ordinary skill in the art. The embodiments are chosen and described in order to better illustrate the principles and practical application of the present application, and to enable those of ordinary skill in the art to understand the present application so as to design various embodiments with various modifications suitable for specific purposes.

[0016] Embodiment: Reference Figs. 1-2 Taking the prediction of lead (Pb) pollution in the urban-rural ecotone soil as an example, the data preparation and preprocessing, and the collection of geographic information data in the multi-source auxiliary variable collection include: Human source variables, obtain industrial enterprise POI data from commercial map platforms, filter out non-ferrous metal smelting, metal products, chemical industry and other heavy industries; obtain traffic road network data. In ArcGIS software, create 400m, 600m, 800m, 1000m, 1200m, 1400m, 1600m, 1800m, 2000m, 2500m, 3000m, 4000m, a total of 12 radius ring buffer zones centered on each sampling point / grid point, and calculate the number and total area of enterprises in each buffer zone as density variables. The reason for choosing this scale is that the typical distance of soil heavy metals (such as Pb) from industrial sources is hundreds of meters (local deposition) to several kilometers (regional migration), and the former small interval ensures fine capture of near-source effects, and the latter large interval avoids redundant calculation and covers the wide-area effect.

[0017] At the same time, calculate the Euclidean distance from each sampling point / grid point to the nearest key industry enterprise and main road as distance type variables. Among them, “key industry enterprises” specifically refer to non-ferrous metal smelting and rolling processing industry, metal products industry, battery manufacturing industry, chemical raw materials and chemical products manufacturing industry, which are considered as key sources due to the emission of heavy metals (such as Pb).

[0018] Natural flow and sink variables, obtain 30m resolution DEM data, derive slope; obtain clay content, organic matter content and other attributes from soil survey data.

[0019] Soil sample collection and testing, based on the spatial distribution of the above variables, 200 surface (0-20cm) soil sampling points are laid out using spatial stratified random sampling method, and coordinates are recorded using GPS. After air drying, grinding and sieving, the content of Pb is measured by inductively coupled plasma mass spectrometry.

[0020] Training data screening, calculate the statistical value of the Pb content of all 200 samples to determine the anomaly threshold. As an implementation, the “mean ± 3 times standard deviation method” can be used, that is, the value of the average plus 3 times the standard deviation is calculated as the upper limit of the anomaly threshold, and the value of the average minus 3 times the standard deviation is calculated as the lower limit of the anomaly threshold, and the sample points with content higher than the upper limit or lower than the lower limit are removed.

[0021] As another optional implementation, the “cumulative frequency method” can also be used to determine the geochemical background value and the anomaly threshold. This method specifically includes the following steps: Firstly, the target heavy metal content values of all samples are sorted in ascending order; then, the cumulative frequency corresponding to each content value is calculated (i.e. the percentage of the number of samples less than or equal to the value in the total number of samples); then, the cumulative frequency curve of the heavy metal content is drawn; usually, the curve will show a clear inflection point, and the low value part before the inflection point mainly represents the geochemical background population. Finally, the content value corresponding to the inflection point, or the content value corresponding to the cumulative frequency of 85% to 95% according to the actual geological background, is set as the upper limit of the abnormal threshold, and the sample points with content higher than the threshold are removed to exclude the interference of possible artificial pollution high value points or abnormal geological high background points, so as to construct an optimized training data set which can better represent the regional natural background pattern. The reason for setting the threshold as the upper limit is that in soil heavy metal analysis, the abnormal value is usually a high value (pollution exceeding the standard), and the low value background represents the natural pattern (reference geochemical statistical method).

[0022] In the optimal scale screening and spatial stratification, the global optimal scale screening includes the Pb content of 197 sample points and all candidate variables (12 scales x (multiple industry types x 2 types of indicators [enterprise number density and enterprise area] ) + other environmental variables); The "multiple industry types" in the application are specifically: non-ferrous metal smelting and rolling processing industry, metal products industry, battery manufacturing industry, chemical raw material and chemical product manufacturing industry. These types are selected based on their typical contribution to soil heavy metal pollution (such as Pb mainly derived from smelting and battery industry).

[0023] After the secondary screening, the 8 core driving variables constituting the final "optimal scale variable set" are specifically: the number of metal products enterprises within the 1000m buffer zone, the land area of battery manufacturing enterprises within the 800m buffer zone, the distance from the nearest trunk road, clay content, organic matter content, elevation, slope, and the distance from the main river. Random forest regression analysis is carried out in the R language environment. The importance is sorted by calculating the average precision reduction percentage of the variables. The results show that for Pb element, the importance of "the number of metal products enterprises within the 1000m buffer zone" is the highest, which is much higher than other scales. Similarly, the optimal scale of variables such as "distance from the nearest trunk road" and "clay content" is determined. In the core driving variable screening, after obtaining the variable set containing the optimal scale of all environmental factors, in order to avoid multiple collinearity between variables and select the core driving factors, secondary screening is carried out: firstly, the variance inflation factor of the variable set is calculated, and the seriously collinear variables are removed; secondly, based on the importance ranking of the random forest model, the top N variables with cumulative contribution degree of importance reaching more than 85% are selected to constitute the final core driving variable set for spatial stratification and prediction, and finally a "optimal scale variable set" containing 8 variables is formed.

[0024] The spatial stratification of driving patterns includes dividing the study area into a grid of 30 m x 30 m, and calculating the values of all variables in the "optimal scale variable set" for each grid. The K-Means clustering analysis of the optimal scale variable data of all grids is performed using the stats package of R language. The elbow rule is used to determine the optimal number of clusters. The clustering results generate three spatially continuous regions. The machine interpretation is combined with the mean characteristics of the core driving variable set in each cluster (partition) after clustering: Partition A (industrial driving dominant area): In this partition, the average values of industrial source variables such as non-ferrous metal smelting enterprises density and battery manufacturing enterprises area are significantly higher than the overall average of the study area (e.g., more than 1 standard deviation above the overall average), while other types of variable values are close to or lower than the average level.

[0025] Partition C (traffic-life source impact area): In this partition, the distance to the main road value is significantly lower than the average level of the study area, while the population density value is relatively high, and its industrial source variable value is at a low level, which effectively distinguishes it from partition A.

[0026] Partition B (natural background-agricultural activity mixed area): In this partition, the values of the above human source variables are at a low level, while the relative importance of natural sink variables (such as clay content) is the highest.

[0027] This definition method based on the relative advantage of variable combination can effectively handle the variable overlap area and scientifically reflect its dominant driving mechanism.

[0028] Accordingly, the spatial stratification of the entire study area into three driving pattern sub-regions is completed.

[0029] As another implementation of the spatial stratification of driving patterns, the spatial autocorrelation analysis of the heavy metal content data in the optimized training data set can be used first. Then, the optimal scale variable set data of all grid points in these spatial aggregation areas is extracted, and the dominant pollution driving pattern is defined through statistical analysis, thereby completing the spatial stratification. The specific method is as follows:

[0030] As an alternative method of spatial stratification of driving patterns, the spatial autocorrelation analysis of the heavy metal content data in the optimized training data set is used first to identify spatial aggregation areas. The specific steps are as follows: First, the heavy metal content data (such as Pb concentration) is logarithmically transformed to achieve normal distribution, and the Kolmogorov-Smirnov test is used to verify it to reduce the impact of data skewness on analysis. The coordinates of the sampling points recorded by GPS are used to construct a spatial weight matrix, and a distance-based neighborhood definition is used, for example, based on the Euclidean distance between sampling points, the points within a threshold distance (e.g., 1000 m based on the scale of the study area) are defined as the neighborhood, otherwise zero.

[0031] Secondly, the global Moran's I index is calculated to assess the overall spatial autocorrelation strength. The global Moran's I is obtained by dividing the sum of the weighted products of the deviations of each sampling point and its neighborhood points' heavy metal content by the sum of all spatial weights (S0), and then dividing by the sum of the squares of all points' deviations.

[0032] Then, the standardized Z-score is calculated, which is the Moran's I minus its expectation value (expectation value is minus one divided by the number of sampling points minus one) divided by the square root of the variance, which is calculated according to the weight matrix. If the absolute value of the Z-score is greater than 1.96 (corresponding to a p-value less than 0.05), there is significant spatial autocorrelation. The spatial correlation map is drawn, which is the curve of Moran's I with distance, to identify the characteristic distance at which Moran's I or Z-score reaches the maximum value, as a scale reference for spatial correlation.

[0033] Then, the local Moran's I is calculated to identify local cluster areas. The local Moran's I is obtained by multiplying the deviation of a single sampling point by the weighted sum of the deviations of its neighborhood points, and then dividing by the overall variance (after standardization). The significance of the local Moran's I (p-value less than 0.05) is tested by Monte Carlo simulation (e.g. 999 permutations), and classified into: high-high cluster (high-value points surrounded by high-value neighborhoods, indicating hotspots), low-low cluster (low-value points surrounded by low-value neighborhoods, indicating coldspots), high-low outlier (high-value points surrounded by low-value neighborhoods), and low-high outlier (low-value points surrounded by high-value neighborhoods). The high-high and low-low areas are defined as spatial clustering areas (e.g. high concentration hotspots correspond to severely polluted areas).

[0034] Subsequently, the environmental variable numerical characteristics of the optimal set of scale variables within these spatial clustering areas are statistically analyzed, such as calculating the mean, median, and standard deviation of the variables within each clustering area, and comparing them with the overall study area. If the industrial source variable (such as non-ferrous metal enterprise density) in a certain clustering area is significantly higher than the mean (more than 1 standard deviation), it is defined as the dominant pollution driving mode as "industrial driving"; similarly, other variables (such as distance from main road, natural sink variables) are analyzed to define other modes (such as "traffic driven" or "natural background driven").

[0035] Accordingly, spatial stratification is completed, and the study area is divided into several sub-areas with relatively homogeneous pollution driving mechanisms. This method can be implemented using GIS software such as the Geostatistical Analyst module of ArcGIS or the spdep package of R language, ensuring the visualization of the analysis results (such as generating a local Moran's I cluster map).

[0036] In the embodiment, the global Moran's I is about 0.45 (Z-score is 3.2, p-value is less than 0.01) for 197 Pb sample points, indicating that there is a positive spatial autocorrelation; the local analysis identifies 3 major high-high hot spots and 2 low-low cold spots, and accordingly, the stratification result is similar to the K-Means method, but more attention is paid to the spatial pattern of heavy metal concentration.

[0037] In the partition modeling, verification and mapping step, the partition construction random forest model includes attributing the 197 sample points to the three partitions according to their spatial positions. A random forest model is constructed for each partition. The sample point Pb content in the partition is used as the dependent variable, and the complete "optimal scale variable set" is used as the independent variable for training. The model parameters are optimized through cross-validation.

[0038] The leave-one-out cross-validation method is used to evaluate the model performance, and the results show that the average determination coefficient R² of the three partition models is 0.76, and the root mean square error RMSE is 4.2 mg / kg. As a comparison, the cross-validation R² of the global single random forest model constructed using all 197 sample points is only 0.59, and the RMSE is 6.8 mg / kg. The precision of the method of the present application is significantly improved.

[0039] Each 30m grid point in the entire study area is input into the corresponding partition random forest model according to its belonging to the driving mode partition, and the Pb content is predicted. Finally, the prediction results of the three partitions are fused to generate a complete and high-precision soil Pb spatial distribution map.

[0040] The soil property indicators such as available phosphorus and clay content in the study area show strong spatial variability, which fully proves that there is significant spatial heterogeneity within the study area. The traditional global model is difficult to depict such complex internal differences, and the present application is precisely to cope with this challenge, thereby providing a theoretical basis for the improvement of the final precision.

[0041] Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art and related fields without creative labor should belong to the scope of protection of the present application. The structures, devices and operation methods not specifically described and explained in the present application, such as without special description and limitation, are implemented according to the conventional means in the art.

Claims

1. A high-precision prediction method for soil heavy metals that integrates optimal scale screening and spatial stratification, characterized in that, Includes the following steps: S1: Multi-source auxiliary variable collection and multi-scale buffer construction. The system collects anthropogenic source variables, natural flow variables and natural sink variables related to soil heavy metal enrichment, and calculates the values ​​of at least one distance-type or density-type anthropogenic source variable at multiple predetermined buffer scales to form a candidate variable set. S2: Collection and analysis of representative soil samples. Sampling points are set up in the target area to collect soil samples and determine the content of target heavy metals. S3: Training data screening based on geochemical background values, calculating the geochemical background values ​​and anomaly thresholds of heavy metals in the soil of the study area, removing sampling points that exceed the anomaly thresholds, and constructing an optimized training dataset; S4: Optimal Scale Selection and Driving Mode Spatial Hierarchy: S4.1: Global optimal scale selection: Using the optimized training dataset, through variable importance assessment or correlation analysis, a globally optimal scale is selected for each variable from the candidate variable set to form an optimal scale variable set; S4.2: Spatial stratification of driving modes: Based on the optimal scale variable set, the entire study area is divided into several sub-regions with relatively homogeneous pollution driving mechanisms. S5: Construct a random forest prediction model by partitioning the region. For each sub-region divided in S4.2, construct an independent random forest regression model and train the model with the heavy metal content of the sampling points in the sub-region as the dependent variable and the optimal scale variable set as the independent variable. S6: Spatial distribution mapping of heavy metals in the study area. The spatial grid point data of the entire study area is input into the corresponding partitioned random forest model according to its sub-region to predict the heavy metal content of each point, and finally fused to generate a high-precision spatial distribution map of the entire area.

2. The method for high-precision prediction of soil heavy metals integrating optimal scale screening and spatial stratification as described in claim 1, characterized in that, The S4.2 driving mode spatial stratification is specifically as follows: using the optimal scale variable set directly, unsupervised clustering analysis is performed on all spatial units in the study area, and each category obtained from the clustering is defined as a pollution driving mode region.

3. The method for high-precision prediction of soil heavy metals integrating optimal scale screening and spatial stratification as described in claim 1, characterized in that, The S4.2 driving mode spatial layering is specifically as follows: First, spatial autocorrelation analysis is performed using the heavy metal content data in the optimized training dataset to identify spatial clustering areas; Subsequently, the numerical characteristics of environmental variables in the optimal scale variable set within these spatial clusters are statistically analyzed to define their dominant pollution-driving modes, thereby completing the spatial stratification.

4. The method for high-precision prediction of soil heavy metals integrating optimal scale screening and spatial stratification as described in claim 1, characterized in that, The human-sourced variables in S1 include the number density of industrial enterprises of different industry types and the area occupied by enterprises, which are obtained based on GIS and POI data.

5. The method for high-precision prediction of soil heavy metals integrating optimal scale screening and spatial stratification as described in claim 1, characterized in that, The predetermined buffer size in S1 is specifically a ring-shaped buffer with 12 radii: 400m, 600m, 800m, 1000m, 1200m, 1400m, 1600m, 1800m, 2000m, 2500m, 3000m, and 4000m.

6. The method for high-precision prediction of soil heavy metals integrating optimal scale screening and spatial stratification as described in claim 1, characterized in that, In step S3, the geochemical background value and anomaly threshold are determined using the mean ± 3 standard deviation method or the cumulative frequency method.

7. The method for high-precision prediction of soil heavy metals integrating optimal scale screening and spatial stratification as described in claim 1, characterized in that, In S4.1, a random forest model is used to evaluate the importance of variables in order to determine the global optimal scale for each variable.

8. The method according to any one of claims 1 to 7, characterized in that, It also includes the following steps: S7: The iterative update mechanism of the model: when new soil sample data or environmental auxiliary variable data are obtained, steps S3 to S5 are repeated to update the optimal scale variable set and spatial stratification results, and the zoning prediction model is retrained.

Images