Method and system for analyzing distribution characteristics of soil heavy metals based on multivariate statistical analysis

Abstract

The application discloses a soil heavy metal distribution characteristic analysis method and system based on multivariate statistical analysis, and relates to the technical field of environmental data analysis.The soil heavy metal distribution characteristic analysis method and system based on multivariate statistical analysis comprises the following steps: S1, preprocessing soil heavy metal distribution characteristic data and terrain migration constraint data; S2, outputting migration marks and generating spatial corresponding sample relations; S3, performing sample participation gating according to covariance structure robustness analysis results, and outputting principal component scores and cluster partition drafts; S4, performing spatial consistency verification analysis, and generating a directional spatial distribution layer and a hot spot candidate area layer; and S5, generating pollution hot spots and treatment partition boundaries.The application solves the problem that, due to the lack of terrain migration constraints and the dependence on isotropic smooth interpolation, the existing analysis method often generates spatial distribution results that are reasonable on the surface but inconsistent with the slope direction, runoff and erosion transport rules, thereby causing deviation in the positioning and partitioning of pollution high-value areas.

Description

Technical Field

[0001] This invention relates to the field of environmental data analysis technology, specifically to a method and system for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis. Background Technology

[0002] Soil heavy metal pollution is characterized by diverse sources, complex migration and transformation processes, and highly uneven spatial distribution, making it a significant concern in ecological environment monitoring and assessment, farmland quality control, mine and industrial site remediation, and regional planning management. In current technologies and industry practices, the analysis of soil heavy metal distribution characteristics typically relies on geochemical surveys and environmental monitoring systems. This involves establishing sampling points within the study area and conducting laboratory tests to generate spatial sample data on heavy metal content. This data is then used in conjunction with geographic information systems for mapping and statistical analysis. With the development of remote sensing mapping, digital terrain modeling, and spatial data processing technologies, more and more studies are using environmental background information such as topography, land use, and meteorological and hydrological data to explain the spatial patterns and diffusion characteristics of heavy metals. This has also driven the development of distribution characteristic analysis from single-concentration mapping to an integrated approach encompassing "statistical comprehensive characterization - spatial inference - risk application," thereby better serving the practical needs of environmental supervision and remediation.

[0003] For example, the invention patent with announcement number CN114443982B discloses a method and system for detecting and analyzing the spatiotemporal distribution characteristics of heavy metals in soil over a large area. It adopts a distributed computer system with a browser-side and server-side architecture to perform grid-based management of the target large area. It integrates multi-objective geochemical surveys, remote sensing, and various environmental and socioeconomic data and standardizes them for database storage. In the spatial dimension, it constructs a regression analysis model and combines regression kriging interpolation for spatial prediction and mapping. In the temporal dimension, it uses a neural network model to predict the trend of change. It improves data support through supplementary sampling and accuracy verification iteration, thereby achieving the effects of analyzing the spatiotemporal distribution characteristics of heavy metals in soil over a large area, continuous mapping, and dynamic monitoring and early warning.

[0004] For example, the invention patent with announcement number CN111581250B discloses a quantitative research method for the change of heavy metals in soil around a mining area with natural factors. It sets up upwind and downwind sampling lines according to the prevailing wind direction and sets up buffer zones in different levels. It samples and measures the heavy metal content at different buffer zones. It quantifies the influence of wind direction and buffer zone topography by calculating the rate of change of content in the downwind direction relative to the upwind direction. It also establishes relationship curves by combining factors such as altitude, slope, buffer zone distance and river distance, and obtains spatial distribution maps by using Kriging interpolation. In this way, it achieves the effect of quantitative characterization and spatial distribution expression of the change of heavy metal content in the surrounding area of ​​the mining area with natural factors.

[0005] While existing soil heavy metal distribution analysis techniques can integrate multi-source monitoring data at the regional scale and generate continuous spatial distribution results, they still primarily rely on statistical correlation analysis and isotropic spatial interpolation methods for spatial inference. These methods often assume spatial proximity of sample points or minimization of regression residuals, failing to adequately characterize the directional migration patterns of pollutants along slope direction, runoff paths, and erosion transport processes under topographical conditions. Furthermore, the lack of a consistency verification mechanism between multivariate statistical analysis and spatial interpolation results regarding physical migration logic means that while the generated spatial distribution maps are mathematically continuous and stable, they are prone to discrepancies with actual pollution migration directions and accumulation characteristics in complex terrain areas. This limits their ability to locate high-pollution areas, delineate boundaries, and ensure the reliability of zoning results, making it difficult to provide sufficient and robust support for refined governance decisions.

[0006] Therefore, in order to address the above problems, there is an urgent need for a method and system for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis. Summary of the Invention

[0007] Technical problems to be solved To address the shortcomings of existing technologies, this invention provides a method and system for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis. This solves the problem that existing analysis methods, due to the lack of topographic migration constraints and reliance on isotropic smooth interpolation, often generate spatial distribution results that are reasonable on the surface but inconsistent with slope aspect, runoff, and erosion transport patterns, leading to deviations in the location and zoning of high-pollution areas.

[0008] Technical solution To achieve the above objectives, this invention provides the following technical solution: a method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis, comprising: S1, collecting soil heavy metal distribution characteristic data and obtaining topographic migration constraint data; preprocessing the soil heavy metal distribution characteristic data and topographic migration constraint data; S2, performing elemental difference topographic migration discrimination on the soil heavy metal distribution characteristic data and topographic migration constraint data, outputting migration markers based on the elemental difference topographic migration discrimination results, and generating spatial corresponding sample relationships constrained by topography; S3, performing covariance structure robustness analysis on the soil heavy metal distribution characteristic data, performing sample participation gating based on the covariance structure robustness analysis results, determining the sample set participating in principal component analysis and cluster partitioning, and outputting principal component scores and initial drafts of cluster partitioning; S4, performing spatial consistency verification analysis on the principal component scores and topographic migration constraint data, generating directional spatial distribution layers and hotspot candidate area layers based on the spatial consistency verification analysis results; S5, fusing the directional spatial distribution layers, hotspot candidate area layers, and initial drafts of cluster partitioning to generate pollution hotspots and remediation zone boundaries, and archiving the data.

[0009] Further, soil heavy metal distribution characteristic data and topographic migration constraint data are collected. The specific preprocessing process for the soil heavy metal distribution characteristic data and topographic migration constraint data is as follows: Soil heavy metal distribution characteristic data is collected, including: sampling point plane coordinate data, sampling depth data, sampling time data, heavy metal element concentration data, and detection limit data; topographic migration constraint data is obtained, including: digital elevation model data, rainfall time series data, and historical rainfall time series data; using coordinate system one and projection transformation algorithms, the sampling point plane coordinate data and digital elevation model data are preprocessed. The data undergoes coordinate benchmark unification processing; connectivity restoration and resolution consistency processing are performed on the digital elevation model data using digital elevation model hole filling and raster resampling algorithms; anomaly identification and replacement processing is performed on heavy metal element concentration data using median absolute deviation anomaly identification algorithm; batch consistency verification and missing data marking are performed on the detection limit data using detection limit consistency verification rules; missing data in rainfall time series data are repaired using a linear interpolation algorithm for missing segments of rainfall; and distribution standardization and linear normalization algorithms are used to normalize and standardize heavy metal element concentration data and detection limit data.

[0010] Further, the specific process for classifying elemental differences in topographic migration based on soil heavy metal distribution characteristic data and topographic migration constraint data is as follows: using sampling time data as the time alignment benchmark, the cumulative rainfall for the time period is obtained by performing sliding accumulation on the rainfall time series data within the rainfall window length; the number of elements is determined by the element dimension of the heavy metal element concentration data; the median of historical rainfall is obtained by taking the median of the historical rainfall time series data within the historical time period constructed with the rainfall window as the scale; the flow direction grid is obtained by using the eight-neighbor flow direction determination algorithm on the digital elevation model data, and the nearest neighbor matching algorithm with projection constraints is performed along the flow direction using the plane coordinate data of the sampling points to obtain a set of sample point pairs, and the same-layer constraint is performed on the set of sample point pairs based on the sampling depth data; the difference between the downstream and upstream sample points of the sample point pair is calculated for the heavy metal element concentration data to form the concentration difference of the sample point pair, and the absolute value of the concentration difference is used in the calculation; The detection limit data is obtained by adding the upstream and downstream corresponding detection limits of the sample point pair to form the composite detection limit of the sample point pair; the ratio of the cumulative rainfall over the period to the median of historical rainfall plus the numerical stability constant is calculated and then incremented by one to obtain the rainfall amplification adjustment term; for each sample point pair in the sample point pair set, all heavy metal elements are traversed according to the heavy metal element index, and the absolute value of the concentration difference of the corresponding sample point pair is divided by the sum of the composite detection limit of the corresponding sample point pair and the numerical stability constant to obtain the detection limit concentration difference; the summation of the detection limit concentration differences of all heavy metal elements within the same sample point pair is divided by the number of heavy metal elements to obtain the multi-element average difference value of the sample point pair; the median of the multi-element average difference values ​​of all sample point pairs in the sample point pair set is taken to obtain the median of the overall migration difference of the sample point pair; the rainfall amplification adjustment term is multiplied by the median of the overall migration difference of the sample point pair and the negative value is taken before performing the natural index calculation to obtain the terrain migration consistency discrimination value.

[0011] Furthermore, the specific process of generating spatially corresponding sample relationships constrained by terrain, based on the terrain migration discrimination results of element differences, is as follows: Real-time comparison of terrain migration consistency discrimination value and terrain migration consistency discrimination threshold: When the terrain migration consistency discrimination value is less than the terrain migration consistency discrimination threshold, a migration inverse logic flag is output, and the migration inverse logic flag is written into the directional spatial distribution generation and inverse logic verification process. During the spatial interpolation calculation process, the isotropic covariance structure is turned off, and only the interpolation calculation path with directional constraints is used to generate the spatial distribution results; A supplementary sampling point suggestion list is generated along the path corresponding to the sample point pair set, and created and archived in the distribution feature database; When the terrain migration consistency discrimination value is greater than or equal to the terrain migration consistency discrimination threshold, a migration consistency flag and a sample point pair set are output, and simultaneously passed to the directional spatial distribution generation and inverse logic verification process as path constraint input.

[0012] Furthermore, the specific process of performing covariance structure robustness analysis on soil heavy metal distribution characteristic data is as follows: A feature sample matrix is ​​constructed from the heavy metal element concentration data and the detection limit data; the robust covariance matrix is ​​obtained by using the minimum covariance determinant robust covariance estimation algorithm on the feature sample matrix; the shrinkage covariance matrix is ​​obtained by using the Redoyt-Wolf shrinkage covariance estimation algorithm on the same feature sample matrix; the kurtosis dominance is obtained by calculating the Fisher kurtosis coefficient of the feature sample matrix according to its elemental dimensions and taking the maximum absolute value; the robust covariance matrix is ​​then calculated. Calculate the determinant of the contracted covariance matrix; add a numerical protection constant to the determinant of the contracted covariance matrix to obtain the stabilized determinant term; divide the determinant of the robust covariance matrix by the stabilized determinant term and add one, then perform a natural logarithm operation on the result to obtain the covariance structure logarithmic difference term; add one to the kurtosis dominance to obtain the kurtosis adjustment term; divide the covariance structure logarithmic difference term by the kurtosis adjustment term to obtain the kurtosis-constrained structural difference; perform a hyperbolic tangent operation on the structural difference to obtain the statistical structural stability value.

[0013] Furthermore, based on the results of the covariance structure robustness analysis, the specific process of implementing sample participation gating and determining the sample set participating in principal component analysis and cluster partitioning, and outputting the principal component score and initial draft of cluster partitioning is as follows: Real-time comparison of statistical structure stability value and statistical structure stability threshold: When the statistical structure stability value is less than the statistical structure stability threshold, it is marked as structurally unstable. After grouping the sample set according to the batch of test reports, the statistical structure stability value is calculated separately. Only the group that meets the requirement that the statistical structure stability value is greater than or equal to the statistical structure stability threshold enters the principal component analysis and cluster partitioning. The group that does not meet the requirement is suspended from participating in the analysis, and the corresponding sample is included in the retesting and quality verification processing path, generating a retesting suggestion list log; When the statistical structure stability value is greater than or equal to the statistical structure stability threshold, principal component analysis and cluster partitioning are performed, and the principal component score and initial draft of cluster partitioning are output.

[0014] Furthermore, the specific process of performing spatial consistency verification analysis on the principal component scores and terrain migration constraint data is as follows: A directional kriging interpolation algorithm is used to obtain the directional interpolated score field; the effective grid cells of the directional interpolated score field constitute a spatial location set; an eight-neighborhood flow direction determination algorithm is used to obtain the flow direction unit vector from the digital elevation model data; a two-dimensional differential gradient calculation algorithm is used to obtain the score field gradient vector; a median absolute deviation algorithm is used to calculate the median absolute deviation of the gradient magnitude from the gradient vector sequence; the direction deviation term is obtained by subtracting the cosine of the score field gradient vector and the flow direction unit vector; the gradient magnitude is then calculated. The absolute deviation of the median is incremented by one, and then the numerical stability constant is added to obtain the gradient robust scale denominator. The directional deviation term is divided by the gradient robust scale denominator to obtain the scale directional deviation. The ratio of the cumulative rainfall over the period to the historical rainfall median plus the numerical stability constant is calculated and incremented by one to obtain the rainfall amplification adjustment term. The scale directional deviation is multiplied by the rainfall amplification adjustment term to obtain the rainfall-adjusted inverse logic metric. The median of the rainfall-adjusted inverse logic metric for each spatial location is taken in the spatial location set to obtain the spatial median inverse logic metric. The spatial median inverse logic metric is negative and the natural exponent operation is performed. The result of the natural exponent operation is subtracted by one to obtain the inverse logic conflict discrimination value.

[0015] Furthermore, the specific process for generating the directional spatial distribution layer and the hotspot candidate region layer based on the spatial consistency verification analysis results is as follows: Real-time comparison of the anti-logic conflict discrimination value and the anti-logic conflict discrimination threshold, performing anti-logic verification processing: When the anti-logic conflict discrimination value is less than the anti-logic conflict discrimination threshold, a directional spatial distribution layer is generated based on the directional interpolation score field and the corresponding spatial grid structure; on the directional spatial distribution layer, regional connectivity analysis is performed based on the principal component scores and spatial adjacency relationships, extracting continuously distributed response regions as hotspot candidate region layers; the directional spatial distribution layer and the hotspot candidate region layer are output; when the anti-logic conflict discrimination value is greater than or equal to the anti-logic conflict threshold, the directional spatial distribution layer and the hotspot candidate region layer are generated ... When determining the threshold for logical conflict, it is marked as anti-logical conflict. The direction-related scale parameters used to control the spatial range along the topographic flow direction and the lateral direction during the directional interpolation process are re-estimated. The re-estimate aims to reduce the anti-logical conflict discrimination value, so that the spatial distribution results after re-estimate gradually converge under the constraints of directional consistency and topographic migration. Based on the relationship between the watershed boundary and flow direction calculated by the digital elevation model, the spatial propagation path across different confluence units in the directional interpolation results is identified, and the propagation path across the watershed is subject to blocking constraints. The directional interpolation score field and anti-logical conflict discrimination value are regenerated, the checklist is output, and the difference layer before and after re-estimate is archived to the distribution feature database.

[0016] Furthermore, the specific process of integrating the directional spatial distribution layer, hotspot candidate area layer, and initial clustering partitioning draft to generate pollution hotspots and governance partition boundaries, and archiving the data, is as follows: Based on the directional spatial distribution layer and hotspot candidate area layer, combined with the initial clustering partitioning draft, the center location of pollution hotspots, the boundary of the hotspot's influence range, and the boundary of the governance partition are generated, and a mapping relationship between the governance partition boundary and the set of sampling points is established; at the same time, the collected data, preprocessing logs, terrain migration consistency discriminant values, statistical structure stability values, anti-logic conflict discriminant values, directional interpolation parameters, and anti-logic reestimation records are archived to the distribution feature database.

[0017] The second aspect of this invention provides a soil heavy metal distribution characteristic analysis system based on multivariate statistical analysis, comprising: a data acquisition and preprocessing module for acquiring soil heavy metal distribution characteristic data and obtaining topographic migration constraint data; preprocessing the soil heavy metal distribution characteristic data and topographic migration constraint data; a topographic migration consistency discrimination module for performing elemental difference topographic migration discrimination on the soil heavy metal distribution characteristic data and topographic migration constraint data, outputting migration markers based on the elemental difference topographic migration discrimination results, and generating spatial correspondence sample relationships constrained by topography; and a constraint-aware multivariate statistical structure extraction module for performing covariance analysis on the soil heavy metal distribution characteristic data. The structural robustness analysis module performs sample participation gating based on the covariance structural robustness analysis results and determines the sample set participating in principal component analysis and cluster partitioning, outputting principal component scores and initial drafts of cluster partitioning. The directional spatial distribution generation and inverse logic verification module performs spatial consistency verification analysis on principal component scores and terrain migration constraint data, generating directional spatial distribution layers and hotspot candidate area layers based on the spatial consistency verification analysis results. The pollution hotspot location, remediation zone delineation, and traceable archiving module integrates the directional spatial distribution layer, hotspot candidate area layer, and initial drafts of cluster partitioning to generate pollution hotspot and remediation zone boundaries, and archives the data.

[0018] Beneficial effects The present invention has the following beneficial effects: (1) This invention strengthens the adaptability of spatial distribution results to real landform and environmental constraints by introducing a collaborative analysis mechanism of topographic spatial structure and environmental background information in the process of analyzing the distribution characteristics of heavy metals in soil. This improves the consistency between spatial distribution results and actual migration patterns and effectively solves the problem of insufficient correlation between spatial analysis results and topographic environment in the prior art.

[0019] (2) By constructing a multivariate statistical structure constraint and analysis process for spatial distribution analysis, this invention comprehensively evaluates the stability and consistency of sample data at the statistical level, thereby enhancing the data quality and structural reliability of the data involved in spatial analysis and effectively solving the problem that the results of multivariate statistical analysis are greatly affected by abnormal data and structural fluctuations in the prior art.

[0020] (3) By introducing directional constraints and consistency discrimination mechanisms in the spatial distribution generation stage, the present invention performs logical verification and adjustment on the spatial inference process, thereby achieving a reasonable balance between directionality and continuity in the spatial distribution expression results, effectively solving the problem of spatial distribution distortion caused by relying on isotropic smooth interpolation in the prior art.

[0021] (4) By linking the spatial distribution analysis results with the process of identifying pollution hotspots and delineating zones, the present invention enables the analysis results to directly support subsequent application stages, thereby realizing the effective transformation of analysis results into application decisions and effectively solving the problem of poor connection between analysis results and actual governance zones in the prior art.

[0022] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0023] Figure 1 This invention relates to a method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis. Figure 2 This invention relates to a soil heavy metal distribution characteristic analysis system based on multivariate statistical analysis. Figure 3 This is a flowchart of the terrain migration consistency determination process of the present invention; Figure 4 This is a comparison diagram of the evolution of the directional index and the anti-logic conflict discriminant value of the present invention; Figure 5 This is a schematic diagram illustrating the relationship between directional deviation and hydrodynamic modulation in this invention. Detailed Implementation

[0024] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0025] Please see Figures 1-5This invention provides a technical solution: a method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis, comprising the following steps: S1, collecting soil heavy metal distribution characteristic data and obtaining topographic migration constraint data; preprocessing the soil heavy metal distribution characteristic data and topographic migration constraint data; S2, performing elemental difference topographic migration discrimination on the soil heavy metal distribution characteristic data and topographic migration constraint data, outputting migration markers based on the elemental difference topographic migration discrimination results, and generating spatial corresponding sample relationships constrained by topography; S3, performing covariance structure robustness analysis on the soil heavy metal distribution characteristic data, performing sample participation gating based on the covariance structure robustness analysis results, determining the sample set participating in principal component analysis and cluster partitioning, and outputting principal component scores and initial drafts of cluster partitioning; S4, performing spatial consistency verification analysis on the principal component scores and topographic migration constraint data, generating directional spatial distribution layers and hotspot candidate area layers based on the spatial consistency verification analysis results; S5, fusing the directional spatial distribution layers, hotspot candidate area layers, and initial drafts of cluster partitioning to generate pollution hotspots and remediation zone boundaries, and archiving the data.

[0026] Specifically, soil heavy metal distribution characteristic data and topographic migration constraint data are collected. The preprocessing process for soil heavy metal distribution characteristic data and topographic migration constraint data is as follows: Soil heavy metal distribution characteristic data are collected, including: sampling point plane coordinate data, sampling depth data, sampling time data, heavy metal element concentration data, and detection limit data. The sampling point plane coordinate data is used to construct the set of sample point pairs and spatial adjacency relationships along the flow direction. The sampling depth data is used for subsequent constraints on samples in the same layer. The sampling time data is used for subsequent rainfall window alignment. The heavy metal element concentration data and detection limit data are used for subsequent discriminant value calculation and multivariate statistical structure extraction.

[0027] The terrain migration constraint data includes: digital elevation model data, time series rainfall data, and historical time series rainfall data. The digital elevation model data is used to calculate the flow direction and constrain the spatial distribution directionality, while the time series rainfall data is used to construct the rainfall-driven adjustment term and participate in consistency and anti-logic verification.

[0028] The plane coordinate data of the sampling points were directly collected by the global satellite navigation system receiver; the sampling depth data were directly collected by the scale of the on-site sampling equipment; the sampling time data were directly collected by the sampling record; the heavy metal element concentration data and the detection limit data were directly obtained from the test report of the testing agency; the digital elevation model data were directly obtained from the UAV airborne lidar mapping results; and the rainfall time series data were directly collected by the meteorological station.

[0029] The coordinate system and projection transformation algorithm are used to unify the coordinate reference of the sampled point plane coordinate data and the digital elevation model data, ensuring consistent overlay indices. The digital elevation model hole-filling algorithm and raster resampling algorithm are used to repair connectivity and unify resolution of the digital elevation model data, ensuring continuous topology for flow direction calculations. The median absolute deviation anomaly identification algorithm is used to identify and replace outliers in heavy metal element concentration data, reducing non-process disturbances to the covariance structure caused by data entry errors and instrument jumps. The detection limit consistency verification rules are used to perform batch consistency verification and missing data marking on the detection limit data. The detection limit consistency verification is performed on a batch-by-batch basis, based on the intra-batch dispersion of the detection limit for the same element. The dispersion threshold is determined by the allowable error of the detection methodology or by statistical analysis of historically stable batches, ensuring that the detection limit constraint terms are consistent. The system traces and repairs missing rainfall time-series data using a linear interpolation algorithm to ensure the calculable rainfall window size. The interpolation interval is determined by adjacent valid observation times, and whether consecutive missing segments are included in the calculation is determined based on the meteorological station sampling frequency and historical missing distribution characteristics. Furthermore, the system normalizes and standardizes heavy metal element concentration data and detection limit data using distribution standardization and linear normalization algorithms. The standardization parameters are derived from the statistical characteristics of the current analysis data, and the normalization target interval is [0, 1]. This eliminates the influence of differences in the dimensions, numerical scales, and detection capabilities of different heavy metal elements on subsequent multivariate statistical analysis and spatial discrimination results, ensuring comparability and numerical stability for each element in covariance estimation, difference measurement, and consistency discrimination. This avoids the non-physical dominance of high-level elements or low-detection-limit elements in the construction of the statistical structure.

[0030] In this implementation plan, by performing unified preprocessing and quality correction on soil heavy metal distribution characteristic data and topographic migration constraint data, data from different sources, scales, and acquisition conditions are consistently expressed at the spatial, temporal, and numerical scale levels, forming a data foundation that is structurally continuous, topologically complete, statistically robust, and traceable. This effectively reduces the interference of non-process factors such as coordinate deviation, resolution inconsistency, outliers, differences in detection limits, and temporal missing data on subsequent analysis results. It ensures that topographic migration constraints, directional spatial distribution analysis, multivariate statistical structure extraction, and discriminant value calculation are all based on reliable, comparable, and physically consistent input data, providing stable data support for subsequent discrimination and analysis.

[0031] Specifically, the process of classifying elemental differences in soil heavy metal distribution characteristics and topographic migration constraints is as follows: Using sampling time data as a time alignment benchmark to unify the reference starting point of rainfall impact and sampling behavior on the time axis, sliding accumulation is performed on the rainfall time series data within the rainfall window length to obtain the cumulative rainfall for that period. The rainfall window length is determined using a multi-scale sliding window correlation analysis method, with a value range between 1 and 30 days. The correlation strength between the cumulative rainfall in different periods and the element sample point concentration difference statistical index is calculated, and the time scale where the correlation strength reaches its maximum value is selected as the rainfall window length to characterize the effective time range of rainfall's impact on soil heavy metal migration. The number of elements is determined by the element dimension of the heavy metal element concentration data. The median of historical rainfall is obtained by taking the median of historical rainfall time series data within the historical period constructed using the rainfall window as the scale, to reduce the interference of extreme rainfall events on the characterization of long-term background rainfall levels. A flow direction grid is obtained from the digital elevation model data using an eight-neighborhood flow direction determination algorithm to explicitly characterize the dominant migration direction of water flow under topographic control, using the plane coordinate data of the sampling points along the flow direction. The projection-constrained nearest neighbor matching algorithm is used to obtain a set of sample point pairs. This algorithm uses the local flow direction as the spatial projection direction and constructs the correspondence between sample points along this direction based on the principle of minimum planar spatial distance. The sample point pair construction process is made subject to the directional consistency constraint of terrain migration. Based on the sampling depth data, a same-layer constraint is applied to the sample point pair set. This same-layer constraint limits the sample point pairs to only samples at the same sampling depth level or with the same soil layer number, avoiding non-physical migration associations introduced by vertical differences. The correlation between downstream and upstream sample points of the sample point pairs is calculated using heavy metal element concentration data. The difference between the concentration differences is used to form the concentration difference between sample pairs, which is used to characterize the spatial variation of elements along the migration path of the terrain. The absolute value of the concentration difference between sample pairs is used in the calculation. The detection limit data are obtained by adding the upstream and downstream corresponding detection limits of the sample pairs to form the composite detection limit of the sample pairs. This is used to construct a detection uncertainty constraint benchmark that matches the concentration difference between sample pairs. The composite detection limit is used to characterize the comprehensive detection uncertainty level of the sample pairs on the corresponding heavy metal elements. It is used to impose scale constraints on the difference in detection capability in the migration difference calculation, so that the sample pairs under different detection conditions are statistically comparable.

[0032] The ratio of cumulative rainfall over a given period to the historical median rainfall plus the numerical stability constant, plus one, yields a rainfall amplification adjustment term. This term characterizes the enhancement or weakening of current rainfall conditions relative to historical background, and explicitly introduces the driving effect of rainfall on heavy metal migration into the migration consistency criterion. For each sample pair in the sample pair set, all heavy metal elements are traversed according to their index. The absolute value of the concentration difference between the corresponding sample pairs is divided by the sum of the lower limit of detection (LD50) and the numerical stability constant of the corresponding sample pairs to obtain the LD50 concentration difference. This suppresses the interference of differences in detection accuracy on migration intensity assessment in multi-element comparisons, ensuring comparability of sample pairs under different detection conditions. The summation of the LD50 concentration differences of all heavy metal elements within the same sample pair, divided by the number of heavy metal elements, yields the multi-element average difference value for the sample pair. By default, the multi-element average difference value uses equal-weighted aggregation to reflect the overall migration differences of the sample pair across multiple heavy metal dimensions. The method is to comprehensively reflect the overall migration difference level of the sample point pair across multiple heavy metal dimensions, while reserving scalability for element participation, thus avoiding the dominance of a single element anomaly in the migration discrimination result. The median of the multi-element average difference values ​​of all sample point pairs in the sample point pair set is taken to obtain the median of the overall migration difference of the sample point pair. This robustly statistically characterizes the typical migration difference features at the regional scale, reducing the impact of local extreme sample point pairs on the overall discrimination. The rainfall amplification adjustment term is multiplied by the median of the overall migration difference of the sample point pair, and the negative value is then applied to the natural exponent calculation. This nonlinear mapping compresses the migration difference scale, ensuring that the migration consistency discrimination result remains numerically stable and monotonically changing, resulting in the topographic migration consistency discrimination value. After exponential mapping, the topographic migration consistency discrimination value falls within a finite interval, used to characterize the monotonic response relationship of migration consistency degree with spatial differences and hydrodynamic conditions, facilitating unified comparison with the migration consistency discrimination logic. The specific calculation formula is as follows: ; In the formula, The topographic migration consistency discriminant value comprehensively characterizes the degree of consistency between the spatial distribution of soil heavy metals and the topographic migration mechanism at the current analysis time. It indicates the number of elements, unifies the dimensions of participation of multiple heavy metal elements in the statistical calculation process, and avoids scale shift of migration difference caused by differences in the number of elements; This represents the cumulative rainfall over a given period, characterizing the overall intensity of rainfall events prior to the sampling time within the selected timescale. It represents the historical median rainfall, used to characterize the background rainfall level at the same time scale, in order to reduce the interference of extreme rainfall events on rainfall driving discrimination; This represents a set of sample pairs, used to define the spatial correspondence between samples participating in migration consistency analysis; It represents the concentration difference between sample points, used to characterize the concentration variation of the same heavy metal element along the topographic migration path between upstream and downstream locations of the sample point pair. This represents the lower limit of detection synthesis amount of the sample pair, used to characterize the level of detection uncertainty of the sample pair for the corresponding heavy metal element; This represents the sample point pair index, a unique number used to identify different spatial correspondences within a set of sample point pairs; This represents the heavy metal element index, used to distinguish the independent participants of different heavy metal elements in the calculation of multi-element migration differences; This represents the analysis time index, used to identify the time reference location corresponding to the terrain migration consistency discriminant value and related statistics; This represents the median operator, used to extract robust typical migration difference statistics at the sample pair set scale to reduce the impact of extreme sample pairs; The numerical stability constant is obtained by using a robust constraint algorithm based on the denominator minimum value of the combined amount of historical rainfall median and sample point detection limit. The value ranges from 0.0005 to 0.001.

[0033] In this implementation scheme, by using sampling time as a unified time reference and introducing rainfall-driven modulation and topographic migration constraints, the spatial variation characteristics of soil heavy metals along the topographically dominated migration path are consistently quantified. This allows the impact of rainfall, topographic control, and differences in multiple heavy metal elements to be comprehensively expressed within the same discrimination framework. This results in a consistent discrimination result of topographic migration that is representative of typical migration behaviors at the regional scale and insensitive to anomalous samples. This effectively characterizes the degree of matching between the current spatial distribution and the topographic migration mechanism, ensuring that the migration discrimination result has physical rationality, statistical robustness, and temporal comparability, and providing a reliable basis for subsequent spatial distribution verification and discrimination decisions.

[0034] Specifically, the process of generating spatially corresponding sample relationships constrained by terrain by outputting migration markers based on element difference terrain migration discrimination results is as follows: Figure 3The flowchart shown illustrates the terrain migration consistency determination process. It compares the terrain migration consistency determination value with the terrain migration consistency determination threshold in real time. When the terrain migration consistency determination value is less than the threshold, it indicates that the current spatial distribution of heavy metals in the soil fails to meet the migration consistency requirements implied by slope aspect control, runoff drive, and erosion transport patterns. A migration inverse logic flag is then output and written into the directional spatial distribution generation and inverse logic verification process. This explicitly triggers the physical constraint enhancement mechanism in the subsequent spatial modeling stage. During spatial interpolation calculations, isotropic covariance structures are disabled, and only interpolation calculation paths with directional constraints are used to generate spatial data. The results of spatial distribution are analyzed; the generation of a spatial distribution layer that appears reasonable on the surface but is physically illogical, based solely on spatial continuity without the constraints of terrain migration, is blocked from the source; a list of suggested supplementary sampling points is generated along the path corresponding to the sample point pair set. The supplementary sampling points are generated based on the spatial distribution density of the sample point pairs on the migration path. When the distance between adjacent sample point pairs exceeds the existing sampling support capacity, supplementary sampling suggestions are generated at the corresponding path position until the end of the path or the end position of the sample point pair coverage. Supplementary sampling suggestions are proposed for positions with sparse samples or insufficient information on the migration path to correct spatial distribution deviations caused by insufficient data support. The results are created and archived in the distribution feature database.

[0035] When the terrain migration consistency discriminant value is greater than or equal to the terrain migration consistency discriminant threshold, a migration consistency marker and a set of sample point pairs are output. This indicates that the current spatial statistical results generally conform to the migration directionality and intensity characteristics under terrain control. The set of sample point pairs is used to limit the spatial corresponding sample relationships participating in multivariate statistical analysis. By clarifying the upstream and downstream, and slope aspect consistent sample associations, physical meaningless cross-slope and counter-current samples are avoided from participating in the statistics. At the same time, it is passed to the directional spatial distribution generation and reverse logic verification processing as path constraint input. The spatial distribution results generated under this constraint can simultaneously meet the requirements of statistical rationality and physical consistency of terrain migration, thereby improving the reliability of locating high-pollution areas and delineating treatment zones.

[0036] In this implementation plan, the determination of whether the spatial distribution results of heavy metals in soil meet the requirements of topographic migration mechanism is transformed into a discriminable and divertable control result. This allows spatial distributions that do not meet migration consistency to be identified and blocked during the generation stage. At the same time, spatial statistical results that meet migration consistency are limited to the sample association range with clear upstream and downstream relationships and slope aspect consistency. Thus, an access and constraint mechanism based on topographic migration law is established throughout the spatial modeling process, ensuring that the final spatial distribution results are consistent in statistical significance and physical mechanism. This effectively improves the reliability and interpretability of identifying and delineating high-pollution areas and treatment zones, and provides a stable and traceable basis for subsequent analysis and decision-making.

[0037] Specifically, the process of performing covariance structure robustness analysis on soil heavy metal distribution characteristic data is as follows: A feature sample matrix is ​​constructed from heavy metal element concentration data and detection limit data. Before entering covariance estimation, the feature sample matrix satisfies the statistical estimability requirements between sample size and feature dimensions, ensuring stable estimation of the covariance structure under the current sample size. Subsequent shrinkage and robust processing reduce the impact of high-dimensional or small sample conditions on the statistical results. By unifying the detection results and detection capability constraints into the same statistical representation space, distortion of the statistical structure between samples due to differences in detection sensitivity is avoided. The feature sample matrix is ​​constructed with sampling points as the sample dimension and heavy metal elements as the feature dimension. During construction, detection limit constraints are introduced. By constraining the data below the detection capability range, the amplified effect of detection noise on covariance estimation is suppressed, forming a standard for multivariate statistical analysis. Standardized sample representation is achieved by using a robust covariance estimation algorithm with minimum covariance determinant to obtain a robust covariance matrix for the feature sample matrix. By automatically selecting the most stable subset in the sample space for covariance estimation, the interference of outliers and extreme samples on the overall covariance structure is effectively reduced. For the same feature sample matrix, a shrinking covariance estimation algorithm with Redoyt-Wolfe is used to obtain a shrinking covariance matrix. An optimal shrinkage tradeoff is introduced between empirical covariance and structured target matrix to improve the numerical stability and generalization ability of covariance estimation under high-dimensional and small sample conditions. The Fisher kurtosis coefficient is calculated for the feature sample matrix according to the element dimension, and the kurtosis dominance is obtained by taking the maximum absolute value. Before participating in subsequent calculations, the kurtosis dominance is scaled by adding one, so that it can be used as a dimensionless adjustment factor to constrain the influence of higher-order distribution forms on the determination of statistical structure, and to characterize the degree of nonnormality of the sample distribution dominated by a single element.

[0038] The determinant of the robust covariance matrix is ​​calculated to characterize the volumetric scale of multivariate features in the overall statistical space, reflecting the global dispersion of the sample covariance structure under anomaly suppression conditions. The determinant of the contracted covariance matrix is ​​also calculated, using the controlled contraction of the covariance structure as a stability benchmark to characterize the scale of the regularized statistical structure under finite sample conditions. A numerical protection constant is added to the determinant of the contracted covariance matrix to obtain a stable determinant term, preventing numerical instability in highly correlated or near-singular conditions. The determinant of the robust covariance matrix is ​​divided by the stable determinant term and incremented by one; the result is then subjected to a natural logarithm operation to obtain the logarithmic difference term of the covariance structure. Mapping and compressing scale differences makes the covariance structure differences exhibit smooth and comparable variations in numerical space. Adding one to the kurtosis dominance yields a kurtosis adjustment term, which explicitly characterizes the influence of higher-order distribution morphology on covariance structure stability, preventing heavy-tailed or leptokurtic distributions from amplifying the structure determination. Dividing the logarithmic covariance structure difference term by the kurtosis adjustment term yields the kurtosis-constrained structural difference quantity. Introducing distribution morphology constraints into the covariance difference measure achieves adaptive suppression of the risk of extreme samples dominating the statistical structure. Performing a hyperbolic tangent operation on the structural difference quantity yields the statistical structure stability value. Bounded mapping of the result using the hyperbolic tangent function ensures that the statistical structure stability value changes continuously within a finite interval. The specific calculation formula is as follows: ; In the formula, It represents the stability value of the statistical structure, which comprehensively characterizes the overall stability of the multivariate statistical structure under the conditions of anomaly suppression and regularization constraints at a given analysis time. The robust covariance matrix is ​​constructed using a robust estimation method to reflect the true covariance structure after suppressing the effects of outliers and extreme samples, emphasizing the robustness of statistical relationships in expressing them. The contracted covariance matrix is ​​used as a stable reference structure obtained by introducing regularization constraints under finite sample conditions, in order to avoid ill-conditioned problems in high-dimensional covariance estimation. It represents the dominance of kurtosis, used to quantify the degree of dominance of higher-order statistical forms in a feature distribution on the overall covariance structure, reflecting the kurtosis and heavy-tailed characteristics of the distribution; This represents the matrix determinant operator, used to characterize the statistical spatial volume scale corresponding to the covariance matrix; This represents the analysis time index, used to distinguish the statistical structure status under different time periods or different data batches; The numerical protection constant is obtained by using the denominator minimization robustness constraint algorithm on the determinant of the robust covariance matrix. The value ranges from 0.0005 to 0.001.

[0039] In this implementation plan, the detection results of heavy metal elements and the constraints of detection capabilities are uniformly incorporated into the same feature representation system before multivariate statistical analysis. By combining anomaly suppression, regularization constraints, and high-order distribution morphology adjustment, the covariance structure among samples is robustly characterized. This ensures that the determination of statistical structure is no longer dominated by individual abnormal samples, detection noise, or numerical instability under high-dimensional small sample conditions. As a result, a stable statistical structure value that is representative of the overall statistical relationship and has adaptive suppression capability for extreme distributions is formed. This ensures that the sample set entering subsequent analysis and partitioning is consistent and reliable in a statistical sense, providing a solid guarantee for the interpretability, reproducibility, and engineering stability of multivariate statistical modeling results.

[0040] Specifically, the process of performing sample participation gating based on the covariance structure robustness analysis results and determining the sample set participating in principal component analysis and cluster partitioning, and outputting the initial draft of principal component scores and cluster partitioning, involves real-time comparison of statistical structure stability values ​​and statistical structure stability thresholds. When the statistical structure stability value is less than the statistical structure stability threshold, it is marked as structurally unstable. This indicates that the multivariate covariance relationship in the current sample set is significantly disturbed by outliers, detection errors, or distribution dispersion. The statistical structure does not yet have the reliability to reflect the true migration and enrichment relationships. The statistical structure stability value is calculated separately for each batch of test reports after grouping the sample set by batch number, test date, and test method identifier issued by the testing institution. Each group must contain a sample size sufficient for minimum statistical calculation to ensure the estimability of the statistical structure stability value. Batch-level splitting reduces the impact of cross-batch systematic bias on the overall statistical structure. The impact of structural determination is mitigated by only including groups whose statistical structural stability values ​​are greater than or equal to the statistical structural stability threshold in principal component analysis and clustering partitioning. This ensures that subsequent feature extraction and partitioning results are based on a statistically consistent and stable subset of samples. Groups that do not meet the criteria are suspended from analysis to prevent them from being misinterpreted as true contamination features in subsequent spatial interpolation and partitioning. The corresponding samples are then included in the retesting and quality verification process, generating a retesting suggestion log. This log must include at least the corresponding batch identifier, the sample number involved, and structural instability markers to support the traceable execution of subsequent retesting and quality verification operations.

[0041] When the statistical structure stability value is greater than or equal to the statistical structure stability threshold, principal component analysis and clustering partitioning are performed. This indicates that the current feature sample matrix has formed a stable and interpretable statistical structure under the joint constraints of robust covariance and contraction covariance, which can be used to reflect the true multi-element covariance pattern. Based on the feature sample matrix, feature decorrelation and dimensionality compression are performed based on the contraction covariance matrix to reduce the interference of multicollinearity and noise dimensions on the spatial partitioning results. The principal component directions that can explain the main covariance structure of the samples are extracted, and the projection results of each sampling point on the principal component directions are used as principal component scores, providing low-noise, high-information-density input variables for subsequent directional spatial analysis. Subsequently, a distance metric between samples is constructed based on the principal component scores, and unsupervised clustering is performed to form partitioning categories with statistical similarity, providing a statistical grouping basis for contaminated spatial partitioning that does not depend on the isotropic assumption. The initial draft of principal component scores and clustering partitioning is output.

[0042] In this implementation plan, a clear structural stability admission mechanism is established before multivariate statistical analysis proceeds to principal component analysis and clustering partitioning by comparing the statistical structural stability value with the statistical structural stability threshold in real time. This ensures that the statistical structure has a discriminative and divisible control result regarding its ability to reflect the true multi-element covariant relationship, thereby effectively isolating abnormal samples, detection errors, and distribution dispersion from interfering with the overall analysis results at the sample level. By implementing batch-level splitting and screening of unstable structures, feature extraction and partitioning analysis are ensured only on a subset of samples with self-consistent and stable statistical structures, avoiding unreliable samples being misinterpreted as true contamination features. At the same time, when the structural stability condition is met, principal component scores and clustering partitioning results with low noise and high information density are output, providing a reliable statistical basis for subsequent spatial distribution analysis and improving the overall reliability, interpretability, and engineering robustness of the contamination partitioning results.

[0043] Specifically, the process of performing spatial consistency verification analysis on principal component scores and terrain migration constraint data is as follows: directional kriging interpolation algorithm is used to obtain the directional interpolated score field; the effective grid cells of the directional interpolated score field constitute a spatial location set; the eight-neighborhood flow direction determination algorithm is used to obtain the flow direction unit vector for the digital elevation model data; the two-dimensional differential gradient calculation algorithm is used to obtain the score field gradient vector; and the median absolute deviation algorithm is used to calculate the median absolute deviation of the gradient magnitude for the score field gradient vector sequence.

[0044] The directional deviation term is obtained by subtracting the cosine of the gradient vector of the score field and the unit vector of the flow direction from one; the absolute deviation of the median gradient magnitude is calculated and incremented by one, then the numerical stability constant is added to obtain the gradient robust scale denominator term; the directional deviation term is divided by the gradient robust scale denominator term to obtain the scale directional deviation; the ratio of the cumulative rainfall over the period to the median of historical rainfall plus the numerical stability constant is calculated and incremented by one to obtain the rainfall amplification adjustment term; the scale directional deviation is multiplied by the rainfall amplification adjustment term to obtain the rainfall-adjusted inverse logic metric; the median of the rainfall-adjusted inverse logic metric corresponding to each spatial location in the spatial location set is taken to obtain the spatial median inverse logic metric; the spatial median inverse logic metric is negative and the natural exponent operation is performed, and the result of the natural exponent operation is subtracted from one to obtain the inverse logic conflict discrimination value.

[0045] Directional kriging interpolation is used to obtain the directional interpolated score field from the principal component scores. By introducing anisotropic structure in the spatial covariance modeling process, the interpolation results have higher continuity and physical rationality along the main migration direction, thus avoiding the weakening of the spatial process directionality by traditional isotropic interpolation. Directional kriging interpolation adopts a preset family of semi-variogram functions and combines directional correlation parameters to characterize spatial correlation. The anisotropic parameters are used to distinguish the spatial influence range along the main migration direction and the lateral direction. The set of spatial locations is composed of the effective grid cells of the directional interpolated score field. Effective grid cells refer to spatial cells with valid values ​​in the interpolation results that are not null values ​​or have no data labels. Null values ​​and boundary cells are excluded. Elements are excluded when constructing the spatial location set to ensure that subsequent spatial consistency and anti-logic verification are only performed within areas with effective statistical significance and spatial coverage. An eight-neighborhood flow direction determination algorithm is used to obtain the flow direction unit vector for the digital elevation model data. By using local elevation differences, the dominant direction of surface runoff is inferred, providing a clear topographic migration reference benchmark for the spatial distribution results. The gradient vector of the score field after directional interpolation is obtained through a two-dimensional differential gradient calculation algorithm to characterize the spatial variation trend and intensity of the principal component scores. The absolute deviation of the median gradient magnitude is obtained by using a median absolute deviation calculation algorithm for the gradient vector sequence of the score field to characterize the overall gradient change scale and suppress scale expansion caused by local extreme gradients or noise.

[0046] The directional deviation term is obtained by subtracting the cosine of the gradient vector and the unit vector of the flow direction from one. This term measures the angular consistency between the gradient change direction and the dominant topographic migration direction, characterizing whether the spatial trend deviates from the natural slope aspect and the main runoff control direction, thus identifying potential physical inconsistency risks. The absolute deviation of the median gradient magnitude is calculated and incremented by one, then the numerical stability constant is added to obtain the gradient robust scale denominator. The median absolute deviation is used to construct a scale constraint term insensitive to local extreme gradients, suppressing excessive amplification of directional deviation caused by interpolation boundary effects or local anomalies. The directional deviation term is divided by the gradient robust scale denominator to obtain the scale directional deviation. The ratio of the cumulative rainfall over a period to the historical median rainfall plus the numerical stability constant is calculated and incremented by one to obtain the rainfall amplification adjustment term. This rainfall amplification adjustment term maintains consistency with the rainfall amplification adjustment term used in the topographic migration consistency judgment in terms of expression, symbol definition, and value domain, used to uniformly characterize the modulating effect of relative rainfall intensity on the migration process. The relative intensity factor of rainfall reflects the amplification effect of hydrodynamic conditions on the migration path and rate of heavy metals, making the directional consistency discrimination process-driven sensitive. Multiplying the scale directional deviation by the rainfall amplification adjustment term yields the rainfall-adjusted inverse logic metric. By coupling spatial directional deviation with the rainfall-driven effect, the higher physical conflict risk corresponding to directional inconsistency under high migration potential conditions is highlighted. The median of the rainfall-adjusted inverse logic metric for each spatial location is taken in the spatial location set to obtain the spatial median inverse logic metric. A median convergence strategy is adopted to reduce the dominant role of local anomalous units on the overall discrimination result, ensuring that the inverse logic assessment reflects the overall regional trend rather than individual noise. The spatial median inverse logic metric is negative and subjected to natural exponentiation. Subtracting the natural exponentiation result from one yields the inverse logic conflict discrimination value. Through exponential mapping, the inverse logic metric is compressed to a finite interval, forming a discrimination value expression with monotonic response characteristics to the degree of directional physical conflict. The specific calculation formula is as follows: ; In the formula, It represents the anti-logic conflict discrimination value, and comprehensively quantifies the degree of deviation between the directional spatial distribution results and the main control mechanism of terrain migration; The resultant field after directional interpolation reflects the spatial response intensity distribution formed under the combined effect of topographic constraints and statistical structure. This represents the flow direction unit vector, used to describe the main direction of local hydrodynamic migration derived from the digital elevation model; This represents the gradient vector of the score field, used to characterize the direction and intensity of the spatial change of the score field after directional interpolation; It represents the absolute deviation of the median gradient magnitude, serving as a robust scale characterization of gradient strength, and is used to mitigate the amplification effect of local anomalous gradients on orientation consistency judgment. It represents the cumulative rainfall over a period of time, used to characterize the overall intensity of hydrodynamic conditions within the analysis period, and has a modulating effect on migration path amplification and spatial redistribution; It represents the historical median rainfall, serving as a robust statistical description of long-term rainfall background levels and used to construct a reference benchmark for relative rainfall intensity; This represents the median operator, which robustly aggregates inverse logical measures at the spatial location set level. This represents a spatial location index, used to identify the specific spatial unit participating in the inverse logic evaluation within the score field after directional interpolation; It represents the set of spatial locations, the range of all valid spatial locations participating in the anti-logic conflict judgment, and is the domain of the spatial statistical aggregation operation; The numerical stability constant is obtained by adding one to the absolute deviation of the median of historical rainfall and the median of gradient magnitude, and applying the minimum value robustness constraint algorithm to the denominator. The value ranges from 0.0005 to 0.001.

[0047] In this embodiment, Table 1 is an example data table for anti-logic conflict discrimination, which records in detail the degree of directional consistency, degree of directional deviation, degree of gradient perturbation, cumulative rainfall, historical rainfall level and the final calculated anti-logic conflict discrimination value under different analysis periods. It is used to quantify the degree of deviation between the directional spatial distribution results and the terrain migration control mechanism. Specifically: for time period T1, the directional deviation is 0.08, the gradient perturbation is 0.15, the cumulative rainfall is 18, the historical rainfall level is 20, and the anti-logic conflict judgment value is 0.07; for time period T2, the directional deviation is 0.25, the gradient perturbation is 0.18, the cumulative rainfall is 22, the historical rainfall level is 20, and the anti-logic conflict judgment value is 0.18; for time period T3, the directional deviation is 0.52, the gradient perturbation is 0.20, the cumulative rainfall is 30, the historical rainfall level is 20, and the anti-logic conflict judgment value is 0.36; for time period T4, the directional deviation is 0.70, the gradient perturbation is 0.25, the cumulative rainfall is 38, the historical rainfall level is 20, and the anti-logic conflict judgment value is 0.58; for time period T5, the directional deviation is 0.85, the gradient perturbation is 0.28, the cumulative rainfall is 45, the historical rainfall level is 20, and the anti-logic conflict judgment value is 0.72.

[0048] Table 1. Example data table for identifying anti-logic conflicts:

[0049] like Figure 4The diagram shows a comparison of the evolution of directional indicators and anti-logic conflict discriminant values. Combined with Table 1, it can be seen that the degree of directional consistency and the degree of directional deviation exhibit opposite trends across different analysis periods, and the anti-logic conflict discriminant value gradually increases with the increase of directional deviation. Specifically, in periods T1 and T2, the degree of directional consistency is high, the degree of directional deviation is low, and the anti-logic conflict discriminant value is at a low level, indicating good consistency between the spatial distribution results of the directional indicators and the main control direction of terrain migration. As we enter periods T3 and T4, the degree of directional consistency decreases significantly, the degree of directional deviation increases significantly, and the anti-logic conflict discriminant value rises synchronously, reflecting a gradual increase in the inconsistency between the spatial distribution results and the terrain migration mechanism. In period T5, the degree of directional deviation reaches a high level, and the anti-logic conflict discriminant value reaches its maximum, indicating that the risk of directional conflict is most prominent at this stage. Overall, this evolutionary comparison diagram intuitively reflects the correspondence between changes in directional indicators and anti-logic conflict discriminant results.

[0050] like Figure 5 The diagram illustrates the relationship between directional deviation and hydrodynamic modulation. By correlating the degree of directional deviation with the anti-logic conflict discrimination value and introducing cumulative rainfall and gradient perturbation as modulation factors, it can be seen that under the same or similar directional deviation conditions, enhanced hydrodynamic conditions further amplify the anti-logic conflict discrimination result. Specifically, during periods with lower cumulative rainfall and smaller gradient perturbation, the anti-logic conflict discrimination value is generally in a lower range. As cumulative rainfall increases and gradient perturbation intensifies, the anti-logic conflict discrimination value exhibits a significant non-linear increasing trend, indicating that hydrodynamic conditions have a significant amplifying and modulating effect on directional conflict. This diagram intuitively reveals the mechanism of change in the anti-logic conflict discrimination value from the perspective of the coupling between directional deviation and hydrodynamics.

[0051] In this implementation scheme, by introducing directional interpolation, topographic migration reference, and rainfall-driven modulation during the spatial modeling stage, the spatial variation trend of principal component scores is physically verified for consistency. This allows the spatial distribution results to be further constrained and judged by the topographic control mechanism in addition to statistical expression, thereby identifying and quantifying potential conflict risks in the spatial distribution that deviate from slope aspect control and runoff migration patterns. Through robust gradient scale constraints and median convergence strategies, the interference of local anomalous gradients, interpolation boundary effects, and noise units on the overall discrimination results is weakened, ensuring that the anti-logic conflict discrimination value stably reflects the degree of directional consistency at the regional scale. By coupling the directional deviation with hydrodynamic conditions, higher sensitivity to the risk of directional inconsistency is ensured under high migration potential conditions. Ultimately, a discrimination result with monotonic response, numerical stability, and physical interpretability is formed for the degree of deviation between the spatial distribution results and the topographic migration control mechanism, providing a reliable basis for subsequent spatial distribution correction, threshold determination, and result publication.

[0052] Specifically, the process of generating a directional spatial distribution layer and a hotspot candidate area layer based on the spatial consistency verification analysis results is as follows: Real-time comparison of the anti-logic conflict discrimination value and the anti-logic conflict discrimination threshold, followed by anti-logic verification processing. When the anti-logic conflict discrimination value is less than the anti-logic conflict discrimination threshold, a directional spatial distribution layer is generated based on the score field after directional interpolation and the corresponding spatial grid structure. During the generation process, the layer maintains the continuity constraint along the terrain flow direction to ensure that the extension direction of the spatial high value area is consistent with the actual slope aspect and confluence direction, avoiding the appearance of regularized and symmetrical pollution spots. On the directional spatial distribution layer, regional connectivity analysis is performed based on the principal component scores and spatial adjacency relationships. The response area is identified by spatial connectivity rather than single-point threshold, suppressing the interference of isolated high values ​​or interpolation noise on the identification of pollution hotspots. The continuously distributed response areas are extracted as hotspot candidate area layers. The directional spatial distribution layer and the hotspot candidate area layer are output.

[0053] When the anti-logic conflict discriminant value is greater than or equal to the anti-logic conflict discriminant threshold, it is marked as an anti-logic conflict, clearly indicating that although the current spatial distribution results are statistically continuous, there is a significant deviation at the directional and topographic migration mechanism level. The directional correlation scale parameters used to control the spatial range along the topographic flow direction and lateral direction in the directional interpolation process are re-evaluated. The re-evaluation aims to reduce the anti-logic conflict discriminant value, so that the re-evaluated spatial distribution results gradually converge under the constraints of directional consistency and topographic migration. The directional correlation scale parameters are used to limit the search radius or influence range along the flow direction and lateral direction in directional kriging interpolation, thereby weakening lateral diffusion and strengthening the constraint strength along the runoff direction, correcting the structural issues caused by isotropic smoothing at the parameter level. Based on the relationship between the watershed boundary and flow direction calculated by the digital elevation model, the structural bias is identified, and spatial propagation paths across different confluence units in the directional interpolation results are identified. The spatial information propagation range is constrained by the topographic structural boundary of the watershed, preventing high pollution values ​​from crossing natural barrier areas during statistical interpolation. Barrier constraints are imposed on cross-watershed propagation paths. Barrier constraints are achieved by treating the watershed boundary as an insurmountable spatial restriction in the interpolation calculation or by introducing penalty weights on cross-boundary propagation paths. The directional interpolation score field and anti-logic conflict discrimination value are regenerated, forming a spatial distribution result after topographic constraint correction and a consistency re-judgment result. A checklist is output and the difference layers before and after reassessment are archived to the distribution feature database.

[0054] This implementation plan establishes a physical consistency verification and self-correction mechanism for directional spatial distribution results, enabling the spatial distribution layer to no longer rely solely on statistical continuity but to be effectively constrained at the directional and topographic migration mechanism levels. When the risk of anti-logical conflict is low, it ensures that the generated directional spatial distribution layer remains continuous along the topographic flow direction and possesses true migration orientation, and stably identifies continuous response areas through regional connectivity analysis, reducing the interference of isolated noise on hotspot determination. When the risk of anti-logical conflict is high, it suppresses lateral non-physical diffusion and cross-confluence unit propagation from both parameter and spatial structure perspectives through adaptive re-estimation of directional scale parameters and watershed boundary constraints, promoting the convergence of spatial distribution results towards the direction that conforms to the main control law of topographic migration, forming traceable and verifiable correction results and difference records, and comprehensively improving the application value of spatial distribution results in terms of physical rationality, result reliability, and engineering controllability.

[0055] Specifically, the process of integrating the directional spatial distribution layer, hotspot candidate area layer, and initial clustering partitioning draft to generate pollution hotspot and remediation zone boundaries, and then archiving the data, is as follows: Based on the directional spatial distribution layer and hotspot candidate area layer, combined with the initial clustering partitioning draft, the center location of the pollution hotspot, the boundary of the hotspot's influence range, and the boundary of the remediation zone are generated. Within the hotspot candidate area, spatial weighted centroid calculation is performed on the principal component scores after directional interpolation to determine the center location of the pollution hotspot. Starting from the center location of the pollution hotspot, the boundary of the hotspot's influence range is determined in the directional spatial distribution layer based on score decay characteristics and spatial connectivity constraints. The influence range boundary is formed by gradually attenuating and truncating the principal component scores along the main expansion direction while satisfying spatial connectivity, thus avoiding the inclusion of isolated units or discontinuous areas. Impact Scope: Within the boundary of the hotspot's impact scope, a preliminary draft of clustering partitioning is introduced to merge and trim spatial fragments, generating spatially continuous governance partition boundaries. The generation of partition boundaries is constrained by maintaining spatial continuity and the principle of minimum segmentation, eliminating partition fragments that are too small in scale or geometrically unstable, and establishing a mapping relationship between the governance partition boundaries and the set of sampling points. The mapping relationship should at least include the spatial identifier of the sampling point, the identifier of the governance partition to which it belongs, and the identifier of the corresponding hotspot. At the same time, the collected data, preprocessing logs, terrain migration consistency discriminant values, statistical structure stability values, anti-logic conflict discriminant values, directional interpolation parameters, and anti-logic reestimation records are archived into the distribution feature database to form a verifiable and reproducible evidence package. The evidence package should at least include the input data index for reproducing the partitioning results, key discriminant value records, and spatial result version identifiers.

[0056] In this implementation plan, by uniformly generating the pollution hotspot center location, impact range boundary, and governance zone boundary based on the directional spatial distribution results and statistical zoning results, the pollution identification results transition from continuous response areas to governance units with clear spatial orientation and engineering boundaries. Furthermore, by combining spatial connectivity and directional attenuation characteristics as constraints, fragmented areas are prevented from being mistakenly classified as independent governance objects, thus forming a structurally complete and clearly defined governance zoning result. By establishing a correspondence between governance zone boundaries and sampling point sets, the zoning results can be directly traced back to the original data source. Simultaneously, key discriminant values, parameter states, and reassessment records throughout the analysis process are systematically archived, constructing a verifiable, traceable, and reproducible evidence package system. This provides stable and reliable data support for subsequent governance decisions, result review, and long-term monitoring.

[0057] like Figure 2 As shown, the second aspect of the present invention provides a soil heavy metal distribution characteristic analysis system based on multivariate statistical analysis, comprising: a data acquisition and preprocessing module for acquiring soil heavy metal distribution characteristic data and obtaining topographic migration constraint data; preprocessing the soil heavy metal distribution characteristic data and topographic migration constraint data; a topographic migration consistency discrimination module for performing elemental difference topographic migration discrimination on the soil heavy metal distribution characteristic data and topographic migration constraint data, outputting migration markers based on the elemental difference topographic migration discrimination results, and generating spatial correspondence sample relationships constrained by topography; and a constraint-aware multivariate statistical structure extraction module for performing covariate analysis on the soil heavy metal distribution characteristic data. The system includes: a robust covariance structure robustness analysis module, which performs sample participation gating based on the covariance structure robustness analysis results and determines the sample set for principal component analysis and cluster partitioning, outputting principal component scores and initial drafts of cluster partitioning; a directional spatial distribution generation and inverse logic verification module, which performs spatial consistency verification analysis on principal component scores and terrain migration constraint data, generating directional spatial distribution layers and hotspot candidate area layers based on the spatial consistency verification analysis results; and a pollution hotspot location, remediation zone delineation, and traceable archiving module, which integrates the directional spatial distribution layer, hotspot candidate area layer, and initial drafts of cluster partitioning to generate pollution hotspot and remediation zone boundaries and archive the data.

[0058] This implementation plan achieves a complete closed-loop control from data acquisition and structural constraint identification to spatial representation and result archiving by collaboratively processing soil heavy metal distribution information and topographic migration constraint information. This ensures that the spatial analysis results possess both statistical consistency and the rationality of topographic migration mechanisms. The overall process effectively suppresses the risk of structural distortion caused by abnormal samples, detection uncertainties, or isotropic assumptions, ensuring that the sample relationships, statistical structures, and spatial distributions involved in the analysis conform to the migration patterns under topographic control. Based on this, the system stably outputs directional spatial distribution results and candidate pollution hotspot areas, and further forms interpretable and verifiable pollution hotspot locations and remediation zoning boundaries. This provides a reliable data foundation and evidence support for subsequent risk assessment, remediation decisions, and result reproduction, effectively improving the scientific rigor and credibility of soil heavy metal spatial analysis and remediation zoning.

[0059] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0060] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis, characterized in that, Includes the following steps: S1. Collect soil heavy metal distribution characteristic data and obtain topographic migration constraint data; preprocess the soil heavy metal distribution characteristic data and topographic migration constraint data. S2 performs elemental difference topographic migration discrimination on soil heavy metal distribution characteristic data and topographic migration constraint data, outputs migration markers based on elemental difference topographic migration discrimination results, and generates spatial correspondence sample relationships constrained by topography. S3 performs covariance structure robustness analysis on soil heavy metal distribution characteristic data, performs sample participation gating based on the covariance structure robustness analysis results, determines the sample set to participate in principal component analysis and cluster partitioning, and outputs principal component scores and initial drafts of cluster partitioning. S4. Perform spatial consistency verification analysis on the principal component scores and terrain migration constraint data, and generate directional spatial distribution layer and hotspot candidate area layer based on the spatial consistency verification analysis results; S5 integrates the directional spatial distribution layer, the hotspot candidate area layer, and the initial draft of the clustering partition to generate pollution hotspots and governance partition boundaries, and then archives the data.

2. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process for collecting soil heavy metal distribution characteristic data and obtaining topographic migration constraint data, and preprocessing the soil heavy metal distribution characteristic data and topographic migration constraint data is as follows: Soil heavy metal distribution characteristic data were collected, including: sampling point plane coordinate data, sampling depth data, sampling time data, heavy metal element concentration data, and detection limit data. Acquire terrain migration constraint data, which includes: digital elevation model data, time-series rainfall data, and historical time-series rainfall data; The coordinate system and projection transformation algorithm are used to unify the coordinate reference of the sampled point plane coordinate data and the digital elevation model data; the digital elevation model hole filling algorithm and raster resampling algorithm are used to repair the connectivity and unify the resolution of the digital elevation model data; the median absolute deviation anomaly identification algorithm is used to identify and replace anomalies in the heavy metal element concentration data; the detection limit consistency verification rule is used to perform batch consistency verification and missing data marking on the detection limit data; the linear interpolation algorithm for missing segments of rainfall is used to repair missing data in the time series of rainfall data; and the distribution standardization and linear normalization algorithms are used to normalize and standardize the heavy metal element concentration data and the detection limit data.

3. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process for determining elemental differences and topographic migration based on soil heavy metal distribution characteristic data and topographic migration constraint data is as follows: Using sampling time data as the time alignment benchmark, the cumulative rainfall for the time period is obtained by performing sliding accumulation on the time series data of rainfall within the rainfall window length; the number of elements is determined by the element dimension of the heavy metal element concentration data; the median of historical rainfall is obtained by taking the median of the historical rainfall time series data within the historical time period constructed with the rainfall window as the scale; the flow direction grid is obtained by using the eight-neighbor flow direction determination algorithm on the digital elevation model data, and the sample point pair set is obtained by performing the projection constraint nearest neighbor matching algorithm along the flow direction using the plane coordinate data of the sampling points, and the same-layer constraint is applied to the sample point pair set based on the sampling depth data; the difference between the downstream and upstream sample points of the sample point pair is calculated for the heavy metal element concentration data to form the sample point pair concentration difference, and the absolute value of the concentration difference is used in the calculation; the detection lower limit data is obtained by adding the corresponding upstream and downstream detection lower limits of the sample point pair to form the sample point pair detection lower limit composite quantity; The rainfall amplification adjustment term is obtained by multiplying the cumulative rainfall over a given period by the median historical rainfall plus the numerical stability constant, and then adding one. For each sample pair in the sample pair set, all heavy metal elements are traversed according to their heavy metal element indices. The absolute value of the concentration difference between the corresponding sample pairs is divided by the sum of the lower limit of detection (LDR) and the numerical stability constant for the corresponding sample pairs to obtain the LDR concentration difference. The LDR concentration differences of all heavy metal elements within the same sample pair are summed and divided by the number of heavy metal elements to obtain the multi-element average difference value of the sample pairs. The median of the multi-element average difference values ​​of all sample pairs in the sample pair set is taken to obtain the median of the overall migration difference of the sample pairs. The rainfall amplification adjustment term is multiplied by the median of the overall migration difference of the sample pairs and the result is negative before performing the natural index calculation to obtain the topographic migration consistency discriminant value.

4. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process of outputting migration markers based on the element difference terrain migration discrimination results and generating spatial correspondence sample relationships constrained by terrain is as follows: Real-time comparison of terrain migration consistency discriminant value and terrain migration consistency threshold: When the terrain migration consistency discriminant value is less than the terrain migration consistency discriminant threshold, a migration inverse logic flag is output. The migration inverse logic flag is written into the directional spatial distribution generation and inverse logic verification process. During the spatial interpolation calculation, the isotropic covariance structure is turned off, and only the interpolation calculation path with directional constraints is used to generate the spatial distribution result. A supplementary sampling point suggestion list is generated along the path corresponding to the sample point set, and it is created and archived into the distribution feature database. When the terrain migration consistency discrimination value is greater than or equal to the terrain migration consistency discrimination threshold, the migration consistency mark and sample point pair set are output and passed to the directional spatial distribution generation and reverse logic verification processing as path constraint input.

5. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process of performing covariance structure robustness analysis on soil heavy metal distribution characteristic data is as follows: A feature sample matrix is ​​constructed from the heavy metal element concentration data and the detection limit data. The robust covariance matrix is ​​obtained by using the minimum covariance determinant robust covariance estimation algorithm on the feature sample matrix. The shrinkage covariance matrix is ​​obtained by using the Redoyt-Wolf shrinkage covariance estimation algorithm on the same feature sample matrix. The kurtosis degree is obtained by calculating the Fisher kurtosis coefficients of the feature sample matrix according to the element dimension and taking the maximum absolute value. Calculate the determinant of the robust covariance matrix; calculate the determinant of the contracted covariance matrix; add a numerical protection constant to the determinant of the contracted covariance matrix to obtain the stabilized determinant term; Divide the determinant of the robust covariance matrix by the stable determinant term and add one, then perform a natural logarithmic operation on the result to obtain the covariance structure logarithmic difference term; add one to the kurtosis dominance to obtain the kurtosis adjustment term; divide the covariance structure logarithmic difference term by the kurtosis adjustment term to obtain the kurtosis-constrained structure difference; perform a hyperbolic tangent operation on the structure difference to obtain the statistical structure stability value.

6. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process of performing sample participation gating based on the covariance structure robustness analysis results, determining the sample set for principal component analysis and cluster partitioning, and outputting the principal component scores and initial draft cluster partitioning is as follows: Real-time comparison of statistical structure stability values ​​and statistical structure stability thresholds: When the statistical structure stability value is less than the statistical structure stability threshold, it is marked as structurally unstable. The sample set is grouped according to the batch of the test report and the statistical structure stability value is calculated separately. Only the group that meets the statistical structure stability value greater than or equal to the statistical structure stability threshold is entered into principal component analysis and cluster partitioning. The group that does not meet the threshold is suspended from participating in the analysis, and the corresponding sample is included in the retest and quality check processing path, and a retest suggestion list log is generated. When the statistical structure stability value is greater than or equal to the statistical structure stability threshold, perform principal component analysis and cluster partitioning, and output the principal component scores and initial drafts of cluster partitioning.

7. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process for performing spatial consistency verification analysis on principal component scores and terrain migration constraint data is as follows: The principal component scores are processed using a directional kriging interpolation algorithm to obtain a directional interpolated score field. A spatial location set is constructed from the effective grid cells of the directional interpolated score field. An eight-neighborhood flow direction determination algorithm is used to determine the flow direction unit vector from the digital elevation model data. The directional interpolated score field is then processed using a two-dimensional differential gradient calculation algorithm to obtain the score field gradient vector. Finally, a median absolute deviation algorithm is used to calculate the median absolute deviation of the gradient magnitude from the gradient vector sequence. The direction deviation term is obtained by subtracting the cosine of the gradient vector of the scoring field and the unit vector of the flow direction from one. Calculate the absolute deviation of the median gradient magnitude and add one, then add the numerical stability constant to obtain the gradient robust scale denominator. Divide the directional deviation term by the gradient robust scale denominator to obtain the scale directional deviation. Calculate the ratio of the cumulative rainfall over the period to the historical rainfall median plus the numerical stability constant and add one to obtain the rainfall amplification adjustment term. Multiply the scale directional deviation by the rainfall amplification adjustment term to obtain the rainfall-adjusted inverse logic metric. Take the median of the rainfall-adjusted inverse logic metric for each spatial location in the spatial location set to obtain the spatial median inverse logic metric. Negate the spatial median inverse logic metric and perform a natural exponent operation. Subtract one from the natural exponent operation result to obtain the inverse logic conflict discrimination value.

8. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process for generating the directional spatial distribution layer and the hotspot candidate region layer based on the spatial consistency verification analysis results is as follows: Real-time comparison of the anti-logic conflict judgment value and the anti-logic conflict judgment threshold, and performance of anti-logic verification processing: When the anti-logic conflict discrimination value is less than the anti-logic conflict discrimination threshold, a directional spatial distribution layer is generated based on the score field after directional interpolation and the corresponding spatial grid structure. On the directional spatial distribution layer, regional connectivity analysis is performed based on the principal component scores and spatial adjacency relationships to extract continuously distributed response regions as hotspot candidate region layers. Output a directional spatial distribution layer and a hotspot candidate region layer; When the anti-logic conflict discrimination value is greater than or equal to the anti-logic conflict discrimination threshold, it is marked as an anti-logic conflict. The direction-related scale parameters used to control the spatial range along the topographic flow direction and the lateral direction in the directional interpolation process are re-estimated. The re-estimated value aims to reduce the anti-logic conflict discrimination value, so that the spatial distribution results after re-estimated gradually converge under the constraints of directional consistency and topographic migration. Based on the watershed boundary and flow direction relationship calculated by the digital elevation model, the spatial propagation path across different confluence units in the directional interpolation results is identified, and the propagation path across the watershed is subject to blocking constraints. The directional interpolation score field and anti-logic conflict discrimination value are regenerated, and the checklist is output and the difference layer before and after re-estimated is archived to the distribution feature database.

9. The method for analyzing the distribution characteristics of heavy metals in soil based on multivariate statistical analysis according to claim 1, characterized in that: The specific process of fusing the directional spatial distribution layer, the hotspot candidate area layer, and the initial clustering partitioning draft to generate pollution hotspot and remediation partition boundaries, and then archiving the data, is as follows: Based on the directional spatial distribution layer and the hotspot candidate area layer, combined with the initial draft of clustering partitioning, the center location of pollution hotspots, the boundary of the hotspot's influence range, and the boundary of the governance partition are generated, and the mapping relationship between the boundary of the governance partition and the set of sampling points is established. At the same time, the collected data, preprocessing logs, terrain migration consistency discriminant values, statistical structure stability values, anti-logic conflict discriminant values, directional interpolation parameters, and anti-logic re-estimation records are archived into the distribution feature database.

10. A soil heavy metal distribution characteristic analysis system based on multivariate statistical analysis, employing the soil heavy metal distribution characteristic analysis method based on multivariate statistical analysis as described in any one of claims 1-9, comprising: The data acquisition and preprocessing module is used to acquire soil heavy metal distribution characteristic data and obtain topographic migration constraint data; and to preprocess the soil heavy metal distribution characteristic data and topographic migration constraint data. The terrain migration consistency discrimination module is used to perform elemental difference terrain migration discrimination on soil heavy metal distribution characteristic data and terrain migration constraint data. Based on the elemental difference terrain migration discrimination results, it outputs migration markers and generates spatial correspondence sample relationships constrained by terrain. The constraint-aware multivariate statistical structure extraction module is used to perform covariance structure robustness analysis on soil heavy metal distribution characteristic data. Based on the covariance structure robustness analysis results, it performs sample participation gating and determines the sample set to participate in principal component analysis and cluster partitioning, and outputs principal component scores and initial drafts of cluster partitioning. The directional spatial distribution generation and anti-logic verification module is used to perform spatial consistency verification analysis on principal component scores and terrain migration constraint data, and generate directional spatial distribution layers and hotspot candidate area layers based on the spatial consistency verification analysis results. The pollution hotspot location, governance zone delineation, and traceable archiving module is used to integrate the directional spatial distribution layer, the hotspot candidate area layer, and the initial draft of the clustering partition to generate the boundaries of pollution hotspots and governance zones, and to archive the data.

Images