Three-dimensional space interpolation method and system for soil pollution based on adaptive grid encryption

Abstract

The application discloses a kind of based on adaptive grid encryption soil pollution three-dimensional space interpolation method, belong to soil environmental information processing technical field, including: obtaining multi-source environmental survey data;Adopt convex hull grid division method to construct the initial overall plane grid framework covering target site;Introduce adaptive grid encryption technology, the space dislocation distance between the calculation node and actual sampling point in the region greater than preset threshold is dynamically increased grid subdivision density, obtains optimized calculation grid;With the optimized calculation grid as carrier, a three-dimensional space interpolation algorithm is used for estimation, to obtain concentration distribution data and estimation confidence;According to the estimation confidence, guide the field to carry out dynamic supplementary survey and obtain new data and iterate until the accuracy requirement is met, and construct soil heavy metal pollution three-dimensional space model.The application improves the accuracy of soil pollution space model construction and heavy metal concentration estimation.

Description

Technical Field

[0001] This invention belongs to the field of soil environmental information processing technology, specifically a three-dimensional spatial interpolation method and system for soil pollution based on adaptive grid encryption. Background Technology

[0002] When detecting heavy metal pollution in soil in areas such as legacy industrial sites, two-dimensional spatial interpolation methods are often used to simulate the planar distribution characteristics of a single heavy metal at a specific depth. However, this approach fails to adequately incorporate the vertical distribution parameters of heavy metals and the spatial correlations between multiple complex heavy metal factors at the data processing level.

[0003] When using three-dimensional spatial interpolation methods to obtain a more comprehensive pollution distribution model, the convex hull mesh method is typically used for initial mesh generation. The computational nodes generated by the convex hull mesh method often have spatial coordinate differences between them and the actual sampling points. This positional deviation between the computational nodes and the actual sampling points causes inconsistencies between the basic input coordinates of the interpolation algorithm and the actual physical coordinates. This misalignment of the basic input coordinates directly leads to errors in subsequent concentration interpolation calculations, thus affecting the accuracy of estimating the volume of contaminated soil at the site when calculating the spatial distribution of pollution.

[0004] Therefore, to address the above problems, a method and system for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption is provided. Summary of the Invention

[0005] To address the aforementioned problems in existing technologies, this invention provides a three-dimensional spatial interpolation method and system for soil pollution based on adaptive grid encryption. This method solves the problem of concentration calculation errors caused by spatial misalignment of input coordinates due to grid division in existing three-dimensional spatial interpolation methods, thereby improving the accuracy of soil pollution spatial model construction and heavy metal concentration estimation.

[0006] The technical solution to achieve the above objectives is: One of the present inventions is a three-dimensional spatial interpolation method for soil pollution based on adaptive grid encryption, comprising: Step S1: Grid-based data collection is performed on the target site to obtain multi-source environmental survey data. The multi-source environmental survey data includes at least the three-dimensional spatial coordinates of the actual sampling points, heavy metal concentration information, and stratigraphic information. Step S2: Based on the spatial distribution of actual sampling points, the convex polygon boundary of the target site is determined by Graham scanning method or convex hull mesh generation method, and the initial overall planar mesh framework covering the target site is constructed by Delaunay triangulation algorithm or regular rectangular mesh generation algorithm. Step S3: Calculate the spatial misalignment distance between the computation nodes and actual sampling points of the initial overall planar mesh frame. Introduce adaptive mesh refinement technology to dynamically increase the mesh subdivision density in areas where the spatial misalignment distance is greater than a preset threshold until the generated computation nodes and actual sampling points match in spatial position to obtain the optimized computation mesh. Step S4: Using the optimized computational grid as the computational carrier, a three-dimensional spatial interpolation algorithm is used in conjunction with heavy metal concentration information to perform spatial estimation, thereby obtaining the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimation confidence level. Step S5: Determine whether the current estimation result meets the accuracy requirements based on the estimation confidence level. If not, extract the substandard area as the target area to guide the field to conduct dynamic supplementary exploration to obtain new sampling point data. Then, integrate the new sampling point data into the multi-source environmental survey data and return to iteratively execute step S2. Step S6: If the accuracy requirements are met, a three-dimensional spatial model of soil heavy metal pollution is constructed based on the concentration distribution data. By accurately quantifying the volume of polluted soil, analyzing the complex overlapping relationship between different heavy metal pollutants, and assessing the physical intrusion characteristics of pollutants, data support is provided for subsequent soil classification, removal, and remediation strategies.

[0007] Preferably, in step S1, data collection is conducted using a grid layout strategy with different densities based on the functional attributes of different historical plots within the target site and their corresponding pollution risk levels. The specific data collection process includes: First, the target site is divided into multiple sub-areas according to different types. Grid points are set up for different areas to obtain a sampling point distribution map. For core production areas with high pollution risk, high-density grid point parameters are configured, while for peripheral areas or general storage areas with low pollution risk, low-density grid point parameters are configured. Secondly, after determining the grid layout parameters for each region, the longitude, latitude, and elevation data of each layout point are obtained through the Global Positioning System or real-time dynamic differential positioning technology to establish a three-dimensional spatial coordinate system for the actual sampling points. At each deployment point, undisturbed soil samples were collected vertically using geological drilling equipment. During the drilling process, soil samples were extracted at preset depth intervals, and the geological characteristics of the borehole cores were recorded in detail, including but not limited to soil color, texture, moisture, density, and the boundary depth of different soil and rock layers. By structuring the undisturbed soil sample data obtained from the drilling, three-dimensional geological sequence data containing the deep structure of each sampling point was generated. Next, the collected undisturbed soil samples were sealed on-site and sent to the laboratory for physicochemical analysis. The heavy metal content in the soil samples was determined by inductively coupled plasma mass spectrometry or atomic absorption spectrometry. The obtained heavy metal concentration information included at least the measured mass fraction of one or more heavy metal elements among arsenic, lead, cadmium, and mercury. Finally, the three-dimensional spatial coordinates, heavy metal concentration information, and stratigraphic information containing three-dimensional geological sequence data are formatted and integrated, and stored in a spatial database to form a multi-source environmental survey data set.

[0008] Preferably, in step S2, constructing an initial overall planar grid framework covering the target site specifically includes: First, after acquiring multi-source environmental survey data, the horizontal projection coordinates of all actual sampling points are extracted. ; Secondly, using this point set as input, the Graham scan method or convex hull mesh generation method is used to determine the convex polygon boundary of the target site; Finally, after determining the boundary of the convex polygon, the Delaunay triangulation algorithm or the regular rectangular meshing algorithm is used inside the boundary to generate an initial overall planar mesh framework consisting of multiple mesh edges and multiple mesh nodes.

[0009] Preferably, in step S3, obtaining the optimized computational grid specifically includes: First, extract the set of computation nodes from the initial global planar mesh framework. and the actual set of sampling points ,against Each computing node in ,exist Search for the nearest actual sampling point And the spatial misalignment distance between the two is calculated using the Euclidean distance formula: ; In the formula, , Representing compute nodes The X-axis and Y-axis coordinates on the horizontal projection plane, , These represent the actual sampling points. X-axis and Y-axis coordinates on the horizontal projection plane; Secondly, the maximum spatial misalignment distance dynamically changes with the site conditions and grid density of different projects, and a preset threshold is required. The maximum allowable spatial misalignment error is adaptively set based on the current grid density, with a preset threshold. The functional variables, dynamically calculated based on the grid partitioning reference side length and sampling point density of the target site, are derived and mathematically defined as follows: Set the reference side length for mesh generation as The spatial distribution density parameter of the actual sampling points in the set area is: Define the spatial tolerance coefficient as Then the preset threshold The calculation formula is: ; In the formula, The damping coefficient is fitted using historical survey data, and its value ranges from [value range missing]. Spatial tolerance coefficient Depending on the sensitivity of the interpolation algorithm used to coordinate offset, its empirical range is strictly limited to [specific values]. When the local encryption process reduces the mesh edge length Reduce or local sampling density When it increases, the preset threshold It exhibits exponential adaptive decay, thereby ensuring the uniqueness and determinism of the convergence boundary of the adaptive encryption process at the underlying logic level; Next, traverse all computing nodes, and when a spatial misalignment distance between a computing node and its corresponding actual sampling point is detected... When it is determined that the local area needs adaptive mesh refinement, the specific process of adaptive refinement is as follows: Add new subdivision nodes inside the grid cell containing the misaligned calculation node. That is, use the quadtree recursive subdivision method to divide a rectangular grid cell with misalignment into four smaller sub-cells and recalculate the distance between the sub-cell nodes and the actual sampling points. By locally and recursively splitting the grid cells and densifying the nodes, the coordinates of the newly generated computational nodes gradually approximate the coordinates of the actual sampling points. This adaptive encryption process iterates until the spatial misalignment distance between all computing nodes and actual sampling points is less than or equal to a preset threshold. At this point, the mesh cell subdivision is stopped, the current mesh topology is output, and it is used as the optimized horizontal computational mesh. Finally, a scheme is introduced to divide the computational grid nodes using three-dimensional geological sequence data. That is, after outputting the optimized horizontal computational grid, the three-dimensional geological sequence data is called, and the optimized horizontal computational grid nodes are divided into layers in the vertical direction based on the elevation of each stratum recorded in the geological borehole data. The coordinates of the 2D grid nodes Substitute the spatial surface equations of each fitted stratigraphic interface In the process, the values ​​of each node at different stratigraphic interfaces were calculated. Axis elevation values, among which, A three-dimensional surface relationship function for a specific stratigraphic interface is generated by spatial interpolation fitting based on the stratigraphic elevations recorded in geological borehole data. Based on the above elevation values, the two-dimensional planar grid is stretched and divided into multiple layers of three-dimensional grid units, so that the nodes of the divided three-dimensional computational grid are attached to the actual undulation of the strata, forming a three-dimensional computational grid base containing X, Y, and Z three-dimensional coordinate attributes. The three-dimensional computational grid base containing X, Y, and Z coordinate attributes is the final, complete, and optimized computational grid.

[0010] Preferably, in step S4, obtaining the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimated confidence level specifically includes: After constructing the optimized computational grid, the measured heavy metal concentration data from the actual sampling points are mapped to the corresponding computational grid nodes. A three-dimensional Kriging interpolation method is used, comprehensively considering the spatial distribution of the actual sampling point data in both the vertical and horizontal directions. Based on the covariance function, a mathematical model of the random field is performed to obtain an unbiased estimate of the heavy metal concentration, i.e.: First, based on the known three-dimensional coordinates and concentration values ​​of the sampling points, the experimental semivariogram function is calculated to describe the degree of variation of the variable with increasing spatial distance. The actual sampling point set is paired and grouped according to different spatial distance steps. Statistical calculations are then performed. The specific physical meaning of the experimental semivariogram function is half of the variance of the heavy metal concentration difference between two sampling points spatially separated by a certain distance. Its specific mathematical formula is as follows: ; In the formula, The experimental semivariance function, The step size is the lag in the three-dimensional spatial distance between two sampling points. Spatial distance equal to The number of actual sampled data pairs, This is the index of the actual sampled data pair. , It is a three-dimensional spatial coordinate vector that includes horizontal, vertical, and elevation dimensions. Located at a three-dimensional coordinate point The measured heavy metal concentration at the location, Located at the deviation point Distance is The measured heavy metal concentration at the location; Secondly, a suitable theoretical semivariance model is selected, and the weighted least squares method is used to fit the experimental semivariance data to obtain the nugget value, sill value, and range parameter. Based on these parameters, a covariance function is constructed to mathematically model the three-dimensional random field. The formula for calculating the covariance function is as follows: ; In the formula, Let covariance function be used. To fit the obtained sill values, which are the theoretical maximum estimated variances of all nodes within the target site, This is the theoretical semivariance value; Finally, for the unknown node to be estimated, based on its three-dimensional spatial distance from surrounding known sampling points and the constructed covariance function, a system of Kriging linear equations is established: ; In the formula, The total number of known sampling points involved in the estimation calculation. Label the known sampling points used in the estimation calculation. , Let be the three-dimensional position vector of the unknown node to be estimated. For the first Three-dimensional position vectors of known sampling points The first one to be solved Known sampling points for unknown estimation points The weighting coefficients, The known points are calculated using the semivariance function described above. and The semivariance between them For known points With the point to be estimated The semivariance between them The introduced Lagrange multipliers are used to ensure the unbiased condition that the sum of the weight coefficients is 1; by performing matrix inversion on the above linear equations, all weight coefficients can be analytically obtained. ; By solving this system of equations, the weight coefficients of each known sampling point are obtained, ensuring that the sum of these weight coefficients is 1 and the estimation variance is minimized, thus achieving an unbiased estimation of heavy metal concentration. After solving the Kriging linear equations to obtain the weighting coefficients, the heavy metal concentration of the node to be estimated is calculated as a weighted sum of the measured concentrations at surrounding known sampling points, using the following formula: ; In the formula, For the node to be estimated The estimated concentration at that location, For the first Measured concentration values ​​of heavy metals at known sampling points; And, calculate the Kriging estimate variance for each grid node, using the following formula: ; In the formula, For known points With the point to be estimated The semivariance between them; Define the theoretical maximum estimated variance of all nodes within the target site as: This value is taken from the sill value parameter obtained when fitting the theoretical semivariance model in step S4. ; For any grid node in space ,in These represent the horizontal, vertical, and elevation coordinates of the 3D spatial grid node in the 3D coordinate system, respectively, and their estimated confidence levels. The calculation formula is as follows: ; In the formula, For grid nodes Kriging estimates the variance at the location.

[0011] Preferably, step S5 specifically includes: Set the global spatial interpolation accuracy threshold The range of values ​​is When any node's When this happens, the dynamic supplementary exploration logic is triggered, that is: Iterate through the estimated confidence data of all 3D computational mesh nodes. If a mesh node exists... If so, the current estimation result is determined to not meet the accuracy requirements; The system uses a spatial clustering algorithm to merge adjacent low-confidence grid nodes and extract a set of substandard 3D spatial regions as target areas for further exploration. For the extracted target area to be explored, the system executes the logic of generating target exploration points, calculates the confidence gradient center or spatial geometric center within the target area, and generates guiding coordinate instructions containing the target longitude, latitude and drilling depth. in, The confidence gradient center within the target area is calculated using a spatial weighted calculation with the difference between the confidence level of each grid node and its attainment threshold as the weighting coefficient. The underlying mathematical derivation and parameter definitions are as follows: The extracted substandard 3D spatial target area contains a total of For any given grid node within the target area, 1 node It is the index number, and Its parameters are defined as follows: Indicates the first The three-dimensional coordinates of each grid node in the real physical space correspond to the X-axis coordinate, Y-axis coordinate, and Z-axis depth in the projection plane, respectively. This represents the first variance calculated using Kriging variance. The estimated confidence level of each grid node; This represents the system's preset global spatial interpolation accuracy threshold. Indicates the first The missing confidence values ​​for each node are defined as follows: The larger this value, the higher the uncertainty of the data at that coordinate point, and the greater its weight in the optimization calculation; Based on the above parameters, the derivation formula for the three-dimensional coordinates of the confidence gradient center of the target area is as follows: ; ; ; In the formula, , and These represent the absolute physical coordinates of the target re-exploration point obtained after weighting by confidence error; The coordinates calculated by the above formula will naturally be biased towards the spatial area with the largest estimation error, and these will be used as priority target points for resurvey, generating guiding instructions that include these coordinates. The instruction is sent to the terminal interface of the on-site exploration equipment via the communication network. On-site personnel conduct supplementary drilling and sampling in the target area according to the coordinate instructions, test and analyze the newly added undisturbed soil samples, and obtain data from the newly added sampling points. The newly added sampling point data is added to the spatial database and merged with the original multi-source environmental survey data. The updated data set is then used to re-trigger the execution of steps S2-S4. This dynamic closed-loop process iterates continuously until all 3D computational mesh nodes are reached. All are greater than or equal to Only when the accuracy requirement is met can the iteration be terminated and the process proceed to step S6.

[0012] Preferably, in step S6, after the three-dimensional interpolation accuracy meets the standard, a spatial entity is generated in the three-dimensional geological modeling environment using the spatial interpolation calculation results, specifically including: Import the digital elevation model data of the target site and generate a three-dimensional surface of the ground using an irregular triangular mesh. Extract the various types of layered stratigraphic interface data processed in step S3, and similarly use irregular triangular mesh to construct three-dimensional undulating surfaces of each stratigraphic interface, such as miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. Using the three-dimensional curved surface of the earth's surface as the top plate and the bottom geological interface as the bottom plate, the geological interfaces in the middle of each layer are used as internal dividing surfaces. Through surface stretching and solid filling calculations along the Z-axis, a multi-layered three-dimensional geological solid model is constructed. Based on this, a three-dimensional spatial model of soil heavy metal pollution is constructed by combining the heavy metal concentration interpolation results, namely: The three-dimensional concentration distribution data of each node output in step S4 are converted into a voxel data structure according to the spatial coordinate mapping rules. A preset pollution concentration threshold is set, and the cells in the voxel model whose concentration exceeds the threshold are defined as pollution voxels. The voxel model consists of voxels with a spatial resolution of [missing information]. Composed of a rectangular grid array, , , These represent the physical side lengths of a single voxel unit along the X, Y, and Z axes in a three-dimensional spatial model, respectively. Define three-dimensional spatial logic state functions ; For the specific heavy metal arsenic, a risk control standard value is set at [value missing]. Extract the three-dimensional concentration distribution data of the nodes. The entity mapping logic judgment condition is expressed as follows: when At that time, assign a state value to the voxel unit. , defined as effective voxels that are contaminated; when At that time, assign a state value to the voxel unit. Defined as a pollution-free blank voxel; The underlying operation rule of Boolean intersection is the logical AND operation of corresponding matrix elements, and the exact relation is: ; In the formula, , , , Representing the coordinates in three-dimensional space respectively Logical state values ​​of single heavy metal pollutants (arsenic, lead, cadmium, and mercury) at the location; This is the Boolean logical AND operator; Only when all single-factor state values ​​at a certain spatial voxel coordinate are equal to 1 The output value is 1, thus accurately extracting the three-dimensional boundary coordinate point cloud set of the complex heavy metal pollution area; Similarly, for other heavy metal factors, namely lead, cadmium, and mercury, their corresponding three-dimensional concentration distribution data at the nodes are extracted. , , And set corresponding risk control standard values. , , Based on the entity mapping logic judgment conditions that are completely consistent with those for arsenic, voxel unit state values ​​corresponding to lead, cadmium, and mercury are generated. , , ; These pollutant voxels are rendered into a three-dimensional geological entity model to form a three-dimensional spatial model of soil heavy metal pollution. Subsequently, using the Boolean logic operation module, including union, intersection, and difference operations, spatial entity calculations are performed on the constructed model to define the contaminated boundary volume. Heavy metals Spatial entities whose concentration exceeds their corresponding risk threshold, among which, To characterize the set subscripts for different heavy metal factors, ,in, Intersection operation is used to extract three-dimensional spatial regions that simultaneously meet the standards for multiple heavy metals exceeding limits, obtaining the three-dimensional boundary and coordinate set of the complex heavy metal pollution region: ; In the formula, A collection of three-dimensional spatial entities representing a region contaminated with complex heavy metals; The difference operation is used to extract the boundary of a spatial region where only a single heavy metal exceeds the standard: ; In the formula, This represents the spatial boundary of a region where the single heavy metal arsenic exceeds the standard. , , , Let each represent a set of spatial entities whose concentrations of arsenic, lead, cadmium, and mercury exceed their respective risk thresholds. The set intersection operator. The set union operator. This is the set difference operator; Based on the single heavy metal pollution area and the complex heavy metal pollution area obtained above, the physical volume of the monomeric element unit is set as follows: Total earthwork volume in areas contaminated with complex heavy metals The calculation formula is: ; In the formula, For three-dimensional coordinates The logical state value of the complex heavy metal pollutant at the location. The physical volume of a monomeric unit.

[0013] Preferably, the three-dimensional spatial model of soil heavy metal pollution, which characterizes the range of pollution concentration, and the three-dimensional geological entity model, which characterizes the physical barrier properties, are registered and superimposed in the same three-dimensional coordinate system to extract the spatial gradient vector in the three-dimensional concentration field and analyze the normal direction of the concentration as the coordinate changes. By comparing the distribution location of high-concentration gradient vectors with the stratigraphic interfaces and weathering layer structure in the 3D geological entity model, the spatial correlation data between the two is output, namely: The logic for extracting and determining spatial location correlation data is based on the inner product operation of the concentration field gradient vector and the geological interface normal vector. The specific parameter definitions and derivation relationships are as follows: This represents a three-dimensional continuous heavy metal concentration space function generated based on the optimized computational grid fitting. This represents the spatial gradient vector of the concentration field of a high-concentration voxel unit in this spatial coordinate system. , , Representing the concentration function The partial derivatives along the X, Y, and Z axes are then defined as follows: The spatial gradient vector is defined as... Its physical meaning is the direction vector of the fastest increase in heavy metal concentration at that spatial point; Simultaneously extract the geological entity model features corresponding to the spatial location of this voxel unit: This represents the implicit surface equation of a specific formation physical interface generated by spatial interpolation fitting of borehole formation elevation data; This represents the surface normal vector at that spatial point, representing the physical interface of a specific stratum. , , Describe the formation interface equations respectively The partial derivatives along the X, Y, and Z axes are then defined as follows: The surface normal vector is defined as... Its physical meaning is the direction perpendicular to the geological sequence boundary; Introducing a correlation calculation model: Defined as the correlation alignment coefficient, it characterizes the parallelism between the concentration diffusion direction and the formation interface normal. Its underlying formula is the absolute value of the cosine of the angle between the gradient vector and the normal vector, i.e.: ; Defined as the spatial correlation threshold, the range of values ​​for the preset dimensionless empirical constant is set as follows. to between; Extract all nodes within a high-concentration voxel region and calculate their correlation alignment coefficient. The average value; When the average value satisfies At that time, the computer logic determines that the direction of the heavy metal concentration gradient is highly coincident with the normal direction of the formation interface, thereby determining the objective spatial relationship that the pollutant is physically infiltrating along a specific formation interface. If the high-concentration area extends along a specific lithological interface in terms of spatial morphology, or is concentrated near the joint and fracture zone of weathered silty mudstone, then the corresponding spatial distribution slice view and coordinate trajectory record are output, thereby revealing the spatial evolution and distribution law of different heavy metal pollutants physically invading downward along the stratigraphic interface and fracture channel.

[0014] A second invention provides a three-dimensional spatial interpolation system for soil pollution based on adaptive grid encryption, comprising: The data acquisition module is used to perform gridded data collection and acquire multi-source environmental survey data containing actual sampling point coordinates and stratigraphic information; The initial framework construction module is used to determine the convex polygon boundaries of the target site using the Graham scan method or the convex hull mesh generation method, and to construct an initial global planar mesh framework covering the target site using the Delaunay triangulation algorithm or the regular rectangular mesh generation algorithm. The adaptive misalignment elimination module is used to calculate the spatial misalignment distance between the computation nodes and the actual sampling points of the initial overall planar mesh frame. It introduces adaptive mesh refinement technology to dynamically increase the mesh subdivision density in areas where the spatial misalignment distance is greater than a preset threshold until the generated computation nodes and the actual sampling points match in spatial position, thus obtaining the optimized computation mesh. The three-dimensional spatial interpolation module is used to perform spatial estimation by combining the three-dimensional spatial interpolation algorithm with heavy metal concentration information, using the optimized computational grid as the computational carrier, to obtain the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimation confidence level. The dynamic closed-loop supplementary exploration module is used to determine whether the current estimation result meets the accuracy requirements based on the estimation confidence level. If it does not meet the requirements, the substandard area is extracted as the target area to guide the field to conduct dynamic supplementary exploration to obtain new sampling point data. The new sampling point data is then integrated into the multi-source environmental survey data and closed-loop iteration is performed. The 3D spatial modeling module is used to construct a 3D spatial model of soil heavy metal pollution based on concentration distribution data after the accuracy standard is met. It provides data support for subsequent soil classification, removal and remediation strategies by accurately quantifying the volume of polluted soil, analyzing the complex overlapping relationship between different heavy metal pollutants and assessing the physical intrusion characteristics of pollutants.

[0015] Compared with the prior art, the beneficial effects of the present invention are: This invention obtains an optimized computational mesh by calculating the spatial misalignment distance between the computational nodes and actual sampling points of the initial overall planar mesh frame, and dynamically increasing the mesh subdivision density in areas where the misalignment distance is greater than a preset threshold. This overcomes the problem of inconsistency between input coordinates and physical coordinates caused by the spatial misalignment between computational nodes and sampling points in traditional convex hull meshes. Spatial estimation was performed by combining an optimized computational grid with a three-dimensional spatial interpolation algorithm. When the estimation results did not meet the accuracy requirements, the substandard areas were extracted to guide dynamic supplementary on-site investigation. After the data of the newly added sampling points were merged, closed-loop iterative calculation was performed, which improved the reliability of the concentration interpolation calculation results and the accuracy of the estimation of the volume of contaminated soil in the site. Furthermore, a three-dimensional spatial model of soil heavy metal pollution was constructed using the concentration distribution data that met the accuracy requirements, which realized the objective restoration of the spatial distribution of site pollution. Attached Figure Description

[0016] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is a flowchart of a three-dimensional spatial interpolation method for soil pollution based on adaptive grid encryption according to the present invention; Figure 2 This is a block diagram of a three-dimensional spatial interpolation system for soil pollution based on adaptive grid encryption according to the present invention; Figure 3 This is a schematic image of the target site and the distribution of each historical plot in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the use of different density grid layouts and geological exploration point distributions for different land parcels in an embodiment of the present invention; Figure 5 This is a magnified view of the grid subdivision node distribution and the spatial misalignment sample error of the actual sampling points in an embodiment of the present invention; Figure 6This is a solid model diagram of the three-dimensional spatial distribution of lead (Pb) element concentration in an embodiment of the present invention; Figure 7 This is a three-dimensional geological entity model of the site constructed by combining multiple geological sections in an embodiment of the present invention; Figure 8 This is a three-dimensional spatial distribution diagram of arsenic (As) exceeding the risk control standard value in an embodiment of the present invention; Figure 9 This is a three-dimensional spatial distribution diagram of lead (Pb) exceeding the risk control standard value in an embodiment of the present invention; Figure 10 This is a three-dimensional spatial distribution diagram of cadmium (Cd) exceeding the risk control standard value in an embodiment of the present invention; Figure 11 This is a three-dimensional spatial distribution diagram of mercury (Hg) exceeding the risk control standard value in an embodiment of the present invention; Figure 12 This is a three-dimensional entity distribution map of a single heavy metal contaminated area extracted through Boolean logic operations in an embodiment of the present invention; Figure 13 This is a three-dimensional solid distribution diagram of the complex heavy metal contaminated area and its corresponding volume result in an embodiment of the present invention; Figure 14 This is a circular diagram showing the statistical proportion of soil volume contaminated with single heavy metals and various types of composite heavy metals in embodiments of the present invention. Figure 15 This is a slice view showing the layered intrusion distribution of different heavy metals along different strata interfaces in an embodiment of the present invention. Detailed Implementation

[0017] 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.

[0018] In existing technologies, the spatial distribution detection of heavy metal pollution in soil from legacy industrial sites primarily employs two-dimensional spatial interpolation methods. This method mainly simulates the planar distribution characteristics of a single heavy metal at a specific depth, without incorporating vertical distribution parameters or the spatial correlation between multiple composite heavy metal factors. When using three-dimensional spatial interpolation methods, the convex hull method is typically used for initial mesh generation. The convex hull method generates spatial coordinate differences between the computational nodes and the actual sampling points, with the maximum misalignment typically reaching approximately 1.6 meters. This positional deviation between the computational nodes and the actual sampling points leads to inconsistencies between the basic input coordinates of the interpolation algorithm and the actual physical coordinates, resulting in errors in the concentration interpolation calculations and affecting the accuracy of the estimated volume of contaminated soil at the site.

[0019] To address the aforementioned technical problems, this invention provides a three-dimensional spatial interpolation method and system for soil pollution based on adaptive grid density. This method eliminates spatial discrepancies between computational nodes and actual sampling points by dynamically adjusting the grid density, and combines three-dimensional spatial interpolation with dynamic supplementary survey iterations to output a three-dimensional spatial pollution model that meets a preset accuracy. The specific implementation process is as follows: like Figure 1 As shown, a three-dimensional spatial interpolation method for soil pollution based on adaptive grid encryption includes: Step S1: Grid-based data collection is performed on the target site to obtain multi-source environmental survey data. The multi-source environmental survey data includes at least the three-dimensional spatial coordinates of the actual sampling points, heavy metal concentration information, and stratigraphic information.

[0020] In this embodiment, based on the functional attributes of different historical plots within the target site and their corresponding pollution risk levels, a grid deployment strategy with different densities is adopted for data collection. The specific data collection process includes: First, such as Figure 3 and Figure 4 As shown, the target site is divided into multiple sub-areas (such as plots) according to different types. Grid points are set up for different areas to obtain a sampling point distribution map. For core production areas with high pollution risk, high-density grid point parameters are configured, such as setting the grid spacing to 100m×100m. For peripheral areas or general storage areas with low pollution risk, low-density grid point parameters are configured, such as setting the grid spacing to 400m×400m. Secondly, after determining the grid layout parameters for each region, the longitude, latitude, and elevation data of each layout point are obtained through the Global Positioning System or real-time dynamic differential positioning technology, and a three-dimensional spatial coordinate system for the actual sampling points is established (for example, using the CGCS2000 coordinate system and the 1985 National Elevation Datum). At each deployment point, undisturbed soil samples are collected vertically using geological drilling equipment. During the drilling process, soil samples are extracted at preset depth intervals (e.g., every 0.5m or 1.0m), and the geological characteristics of the borehole core are recorded in detail, including but not limited to soil color, texture, moisture, density, and the boundary depth of different soil and rock layers. By structuring the undisturbed soil sample data obtained from the drilling, three-dimensional geological sequence data containing the deep structure of each sampling point is generated. Next, the collected undisturbed soil samples were sealed on-site and sent to the laboratory for physicochemical analysis. The heavy metal content in the soil samples was determined by inductively coupled plasma mass spectrometry or atomic absorption spectrometry. The heavy metal concentration information obtained included at least the measured mass fraction of one or more heavy metal elements among arsenic (As), lead (Pb), cadmium (Cd), and mercury (Hg). Finally, the three-dimensional spatial coordinates, heavy metal concentration information, and stratigraphic information containing three-dimensional geological sequence data are formatted and integrated, and stored in a spatial database to form a multi-source environmental survey data set.

[0021] Step S2: Based on the spatial distribution of actual sampling points, the convex polygon boundary of the target site is determined by Graham scanning method or convex hull mesh generation method, and the initial overall planar mesh framework covering the target site is constructed by Delaunay triangulation algorithm or regular rectangular mesh generation algorithm.

[0022] In this embodiment, constructing an initial overall planar grid framework covering the target site specifically includes: First, after acquiring multi-source environmental survey data, the horizontal projection coordinates of all actual sampling points are extracted. ; Secondly, using the point set as input, the Graham scan method or convex hull mesh generation method is used to determine the convex polygon boundary of the target site. Taking the Graham scan method as an example, firstly, the point with the smallest ordinate in the point set is found as the reference point, and all other points are sorted according to the polar angle of the line connecting them to the reference point. Then, a node stack is established, and each point is traversed in the sorted order. By calculating the vector cross product, it is determined whether the path formed by the current point and the two points at the top of the stack is a left turn. If it is a right turn, the top element of the stack is popped. This process continues until all points have been traversed. The node sequence retained in the stack constitutes the minimum convex polygon boundary of the target site. Finally, after determining the boundary of the convex polygon, the Delaunay triangulation algorithm or the regular rectangular mesh partitioning algorithm is used inside the boundary to generate an initial overall planar mesh framework composed of multiple mesh edges and multiple mesh nodes. The boundary of this initial mesh framework is limited by the outermost coordinates of the sampling point set, avoiding data extrapolation calculations in the external area without sampling point coverage. At this time, the generated computational mesh nodes are coordinate points calculated based on pure geometric topology rules.

[0023] Step S3: Calculate the spatial misalignment distance between the computation nodes and actual sampling points of the initial overall planar mesh frame. Introduce adaptive mesh refinement technology to dynamically increase the mesh subdivision density in areas where the spatial misalignment distance is greater than a preset threshold until the generated computation nodes and actual sampling points match in spatial position, and obtain the optimized computation mesh.

[0024] In this embodiment, obtaining the optimized computational grid specifically includes: First, extract the set of computation nodes from the initial global planar mesh framework. and the actual set of sampling points ,against Each computing node in ,exist Search for the nearest actual sampling point And the spatial misalignment distance between the two is calculated using the Euclidean distance formula: ; In the formula, , Representing compute nodes The X-axis and Y-axis coordinates on the horizontal projection plane, , These represent the actual sampling points. The X and Y coordinates on the horizontal projection plane; such as Figure 5 As shown, the left side shows the overall grid partitioning and node distribution of the target site, while the right side shows a magnified view of the local area, which intuitively demonstrates the spatial misalignment and sample error between the computational nodes generated by the convex hull grid method and the actual sampling points. It is precisely to eliminate this error that the subsequent adaptive grid densification logic is triggered. Secondly, the maximum spatial misalignment distance dynamically changes with the site conditions and grid density of different projects, and a preset threshold is required. The maximum allowable spatial misalignment error is adaptively set based on the current grid density. For example, in a region with a grid density of 100m × 100m, a preset threshold is used. The value can be set to 0.1m; in areas with a grid density of 400m×400m, a preset threshold is used. It can be set to 0.5m, with a preset threshold. The functional variables, dynamically calculated based on the grid partitioning reference side length and sampling point density of the target site, are derived and mathematically defined as follows: Set the reference side length for mesh generation as The spatial distribution density parameter of the actual sampling points in the set area is: Define the spatial tolerance coefficient as Then the preset threshold The calculation formula is: ; In the formula, The damping coefficient is fitted using historical survey data, and its value ranges from [value range missing]. Spatial tolerance coefficient Depending on the sensitivity of the interpolation algorithm used to coordinate offset, its empirical range is strictly limited to [specific values]. When the local encryption process reduces the mesh edge length Reduce or local sampling density When it increases, the preset threshold It exhibits exponential adaptive decay, thereby ensuring the uniqueness and determinism of the convergence boundary of the adaptive encryption process at the underlying logic level; Next, traverse all computing nodes, and when a spatial misalignment distance between a computing node and its corresponding actual sampling point is detected... When it is determined that the local area needs adaptive mesh refinement, the specific process of adaptive refinement is as follows: Add new subdivision nodes inside the grid cell containing the misaligned calculation node. That is, use the quadtree recursive subdivision method to divide a rectangular grid cell with misalignment into four smaller sub-cells and recalculate the distance between the sub-cell nodes and the actual sampling points. By locally and recursively splitting the grid cells and densifying the nodes, the coordinates of the newly generated computational nodes gradually approximate the coordinates of the actual sampling points. This adaptive encryption process iterates until the spatial misalignment distance between all computing nodes and actual sampling points is less than or equal to a preset threshold. At this point, the mesh cell subdivision is stopped, the current mesh topology is output, and it is used as the optimized horizontal computational mesh. Finally, a scheme for dividing the computational grid nodes using three-dimensional geological sequence data is introduced. Specifically, after outputting the optimized horizontal computational grid, the three-dimensional geological sequence data is called, and based on the elevations of various strata recorded in the geological borehole data (such as the elevation of the bottom plate of miscellaneous fill, the bottom plate of clay layer, and the top plate of weathered silty mudstone), the optimized horizontal computational grid nodes are layered and segmented in the vertical direction (Z-axis). The coordinates of the 2D grid nodes Substitute the spatial surface equations of each fitted stratigraphic interface In the process, the values ​​of each node at different stratigraphic interfaces were calculated. Axis elevation values, among which, A three-dimensional surface relationship function for a specific stratigraphic interface is generated by spatial interpolation fitting based on the stratigraphic elevations recorded in geological borehole data. Based on the above elevation values, the two-dimensional planar grid is stretched and divided into multiple layers of three-dimensional grid units, so that the nodes of the divided three-dimensional computational grid are attached to the actual undulation of the strata, forming a three-dimensional computational grid base containing X, Y, and Z three-dimensional coordinate attributes. The three-dimensional computational grid base containing X, Y, and Z coordinate attributes is the final, complete, and optimized computational grid.

[0025] Step S4: Using the optimized computational grid as the computational carrier, a three-dimensional spatial interpolation algorithm is used in conjunction with heavy metal concentration information to perform spatial estimation, thereby obtaining the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimation confidence level.

[0026] In this embodiment, obtaining the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimated confidence level specifically includes: After constructing the optimized computational grid, the measured heavy metal concentration data from the actual sampling points are mapped to the corresponding computational grid nodes. A three-dimensional Kriging interpolation method is used, comprehensively considering the spatial distribution of the actual sampling point data in both the vertical and horizontal directions. Based on the covariance function, a mathematical model of the random field is performed to obtain an unbiased estimate of the heavy metal concentration, i.e.: First, based on the known three-dimensional coordinates and concentration values ​​of the sampling points, the experimental semivariogram function is calculated to describe the degree of variation of the variable with increasing spatial distance. The actual sampling point set is paired and grouped according to different spatial distance steps. Statistical calculations are then performed. The specific physical meaning of the experimental semivariogram function is half of the variance of the heavy metal concentration difference between two sampling points spatially separated by a certain distance. Its specific mathematical formula is as follows: ; In the formula, The experimental semivariance function, The step size is the lag in the three-dimensional spatial distance between two sampling points. Spatial distance equal to The number of actual sampled data pairs, This is the index of the actual sampled data pair. , It is a three-dimensional spatial coordinate vector that includes horizontal, vertical, and elevation dimensions. Located at a three-dimensional coordinate point The measured heavy metal concentration at the location, Located at the deviation point Distance is The measured heavy metal concentration at the location; Secondly, a suitable theoretical semivariance model is selected, and the weighted least squares method is used to fit the experimental semivariance data to obtain the nugget value, sill value, and range parameter. Based on these parameters, a covariance function is constructed to mathematically model the three-dimensional random field. The formula for calculating the covariance function is as follows: ; In the formula, Let covariance function be used. To fit the obtained sill values, which are the theoretical maximum estimated variances of all nodes within the target site, This is the theoretical semivariance value; Finally, for the unknown node to be estimated, based on its three-dimensional spatial distance from surrounding known sampling points and the constructed covariance function, a system of Kriging linear equations is established: ; In the formula, The total number of known sampling points involved in the estimation calculation. Label the known sampling points used in the estimation calculation. , Let be the three-dimensional position vector of the unknown node to be estimated. For the first Three-dimensional position vectors of known sampling points The first one to be solved Known sampling points for unknown estimation points The weighting coefficients, The known points are calculated using the semivariance function described above. and The semivariance between them For known points With the point to be estimated The semivariance between them The introduced Lagrange multipliers are used to ensure the unbiased condition that the sum of the weight coefficients is 1; by performing matrix inversion on the above linear equations, all weight coefficients can be analytically obtained. ; By solving this system of equations, the weight coefficients of each known sampling point are obtained, ensuring that the sum of these weight coefficients is 1 and the estimation variance is minimized, thus achieving an unbiased estimation of heavy metal concentration. After solving the Kriging linear equations to obtain the weighting coefficients, the heavy metal concentration of the node to be estimated is calculated as a weighted sum of the measured concentrations at surrounding known sampling points, using the following formula: ; In the formula, For the node to be estimated The estimated concentration at that location, For the first Measured concentration values ​​of heavy metals at known sampling points; And, calculate the Kriging estimate variance for each grid node, using the following formula: ; In the formula, For known points With the point to be estimated The semivariance between them; Define the theoretical maximum estimated variance of all nodes within the target site as: This value is taken from the sill value parameter obtained when fitting the theoretical semivariance model in step S4. ; For any grid node in space ,in These represent the horizontal, vertical, and elevation coordinates of the 3D spatial grid node in the 3D coordinate system, respectively, and their estimated confidence levels. The calculation formula is as follows: ; In the formula, For grid nodes Kriging estimates the variance at the location.

[0027] Step S5: Determine whether the current estimation result meets the accuracy requirements based on the estimation confidence level. If not, extract the substandard area as the target area to guide the on-site dynamic supplementary exploration to obtain new sampling point data. Then, integrate the new sampling point data into the multi-source environmental survey data and return to iterative execution step S2.

[0028] In this embodiment, a global spatial interpolation accuracy threshold is set. The range of values ​​is When any node's When this happens, the dynamic supplementary exploration logic is triggered, that is: Iterate through the estimated confidence data of all 3D computational mesh nodes. If a mesh node exists... If so, the current estimation result is determined to not meet the accuracy requirements; The system uses a spatial clustering algorithm to merge adjacent low-confidence grid nodes and extract a set of substandard 3D spatial regions as target areas for further exploration. For the extracted target area to be explored, the system executes the logic of generating target exploration points, calculates the confidence gradient center or spatial geometric center within the target area, and generates guiding coordinate instructions containing the target longitude, latitude and drilling depth. in, The confidence gradient center within the target area is calculated using a spatial weighted calculation with the difference between the confidence level of each grid node and its attainment threshold as the weighting coefficient. The underlying mathematical derivation and parameter definitions are as follows: The extracted substandard 3D spatial target area contains a total of For any given grid node within the target area, 1 node It is the index number, and Its parameters are defined as follows: Indicates the first The three-dimensional coordinates of each grid node in the real physical space correspond to the X-axis coordinate, Y-axis coordinate, and Z-axis depth in the projection plane, respectively. This represents the first variance calculated using Kriging variance. The estimated confidence level of each grid node; This represents the system's preset global spatial interpolation accuracy threshold. Indicates the first The missing confidence values ​​for each node are defined as follows: The larger this value, the higher the uncertainty of the data at that coordinate point, and the greater its weight in the optimization calculation; Based on the above parameters, the derivation formula for the three-dimensional coordinates of the confidence gradient center of the target area is as follows: ; ; ; In the formula, , and These represent the absolute physical coordinates of the target re-exploration point obtained after weighting by confidence error; The coordinates calculated by the above formula will naturally be biased towards the spatial area with the largest estimation error, and these will be used as priority target points for resurvey, generating guiding instructions that include these coordinates. The instruction is sent to the terminal interface of the on-site exploration equipment via the communication network. On-site personnel conduct supplementary drilling and sampling in the target area according to the coordinate instructions, test and analyze the newly added undisturbed soil samples, and obtain data from the newly added sampling points. The newly added sampling point data is added to the spatial database and merged with the original multi-source environmental survey data. The updated data set is then used to re-trigger the process of steps S2-S4. This dynamic closed-loop process iterates continuously until all 3D computational mesh nodes are reached. All are greater than or equal to Only when the accuracy requirement is met can the iteration be terminated and the process proceed to step S6.

[0029] Step S6: If the accuracy requirements are met, a three-dimensional spatial model of soil heavy metal pollution is constructed based on the concentration distribution data. By accurately quantifying the volume of polluted soil, analyzing the complex overlapping relationship between different heavy metal pollutants, and assessing the physical intrusion characteristics of pollutants, data support is provided for subsequent soil classification, removal, and remediation strategies.

[0030] In this embodiment, after the three-dimensional interpolation accuracy meets the standard, spatial entities are generated in the three-dimensional geological modeling environment using the spatial interpolation calculation results, specifically including: Import the digital elevation model data of the target site and generate a three-dimensional surface of the ground using an irregular triangular mesh. Extract the various types of strata interface data processed in step S3, and similarly use irregular triangular mesh to construct three-dimensional undulating surfaces of each strata interface, such as miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. Using the three-dimensional curved surface of the earth's surface as the top plate and the bottom geological interface as the bottom plate, the geological interfaces in the middle of each layer are used as internal dividing surfaces. Through surface stretching and solid filling calculations along the Z-axis, a multi-layered three-dimensional geological solid model is constructed. Based on this, a three-dimensional spatial model of soil heavy metal pollution was constructed by combining the heavy metal concentration interpolation results, such as... Figure 6 and Figure 7 As shown, a comprehensive solid model is presented, combining multiple geological sections of the site (such as fill, clay, weathered silty mudstone, etc.) with the three-dimensional spatial distribution of lead (Pb) element concentration. The 3D concentration distribution data of each node output in step S4 is converted into a voxel data structure according to the spatial coordinate mapping rules. A preset pollution concentration threshold is set, and the cells in the voxel model whose concentration exceeds the threshold are defined as pollution voxels. The voxel model consists of voxels with a spatial resolution of [missing information]. Composed of a rectangular grid array, , , These represent the physical side lengths of a single voxel unit along the X, Y, and Z axes in a three-dimensional spatial model, respectively. Define three-dimensional spatial logic state functions ; For the specific heavy metal arsenic, a risk control standard value is set at [value missing]. Extract the three-dimensional concentration distribution data of the nodes. The entity mapping logic judgment condition is expressed as follows: when At that time, assign a state value to the voxel unit. , defined as effective voxels that are contaminated; when At that time, assign a state value to the voxel unit. Defined as a pollution-free blank voxel; The underlying operation rule of Boolean intersection is the logical AND operation of corresponding matrix elements, and the exact relation is: ; In the formula, , , , Representing the coordinates in three-dimensional space respectively Logical state values ​​of single heavy metal pollutants (arsenic, lead, cadmium, and mercury) at the location; This is the Boolean logical AND operator; Only when all single-factor state values ​​at a certain spatial voxel coordinate are equal to 1 The output value is 1, thus accurately extracting the three-dimensional boundary coordinate point cloud set of the complex heavy metal pollution area; in the visual display of the model, different color scales are assigned according to different heavy metal types and concentration ranges to realize the digital restoration of the pollution distribution of different plots and different depths. Similarly, for other heavy metal factors, namely lead, cadmium, and mercury, their corresponding three-dimensional concentration distribution data at the nodes are extracted. , , And set corresponding risk control standard values. , , Based on the entity mapping logic judgment conditions that are completely consistent with those for arsenic, voxel unit state values ​​corresponding to lead, cadmium, and mercury are generated. , , ; Based on the above logical state function and entity mapping logic, such as Figure 8 , Figure 9 , Figure 10 and Figure 11 The three-dimensional spatial distribution entities of each heavy metal factor exceeding the corresponding risk control standard value can be extracted and rendered separately, namely pollutant voxels. These pollutant voxels are then rendered into a three-dimensional geological entity model to form a three-dimensional spatial model of soil heavy metal pollution. Subsequently, using the Boolean logic operation module, including union, intersection, and difference operations, spatial entity calculations are performed on the constructed model to define the contaminated boundary volume. Heavy metals Spatial entities whose concentration exceeds their corresponding risk threshold, among which, To characterize the set subscripts for different heavy metal factors, ,in, Intersection operation is used to extract three-dimensional spatial regions that simultaneously meet the standards for multiple heavy metals exceeding limits, obtaining the three-dimensional boundary and coordinate set of the complex heavy metal pollution region: ; In the formula, A collection of three-dimensional spatial entities representing a region contaminated with complex heavy metals; The difference operation is used to extract the boundary of a spatial region where only a single heavy metal exceeds the standard: ; In the formula, This represents the spatial boundary of a region where the single heavy metal arsenic exceeds the standard. , , , Let each represent a set of spatial entities whose concentrations of arsenic, lead, cadmium, and mercury exceed their respective risk thresholds. The set intersection operator. The set union operator. This is the set difference operator; By summing and integrating the number of voxel units extracted from each polluted entity, and multiplying by the physical volume of each voxel unit, the soil volume data of single heavy metal pollution areas and complex heavy metal pollution areas, as well as their proportion of the total polluted soil volume, are output. This provides data support for the classification and removal of polluted soil. Figure 12 and Figure 13 As shown, the distribution of entities in single-heavy-metal-contaminated areas extracted through Boolean logic operations, and the distribution of entities in areas with complex heavy-metal pollution (intersection of multiple heavy-metal factors) and their corresponding volume results are presented intuitively. Simultaneously, earthwork volume data generated based on the above Boolean operation results are shown below. Figure 14 As shown, a detailed statistical chart of the volume percentage of soil contaminated with single heavy metals and various types of compound heavy metals was output, accurately quantifying the amount of site contamination. Based on the single heavy metal pollution area and the complex heavy metal pollution area obtained above, the physical volume of the monomeric element unit is set as follows: Total earthwork volume in areas contaminated with complex heavy metals The calculation formula is: ; In the formula, For three-dimensional coordinates The logical state value of the complex heavy metal pollutant at the location. The physical volume of a monomeric unit.

[0031] The three-dimensional spatial model of soil heavy metal pollution, which characterizes the range of pollution concentration, and the three-dimensional geological entity model, which characterizes the physical barrier properties, are registered and superimposed in the same three-dimensional coordinate system. The spatial gradient vector in the three-dimensional concentration field is extracted, and the normal direction of the concentration as the coordinate changes is analyzed. By comparing the distribution location of high-concentration gradient vectors with the stratigraphic interfaces and weathering layer structure in the 3D geological entity model, the spatial correlation data between the two is output, namely: The logic for extracting and determining spatial location correlation data is based on the inner product operation of the concentration field gradient vector and the geological interface normal vector. The specific parameter definitions and derivation relationships are as follows: This represents a three-dimensional continuous heavy metal concentration space function generated based on the optimized computational grid fitting. This represents the spatial gradient vector of the concentration field of a high-concentration voxel unit in this spatial coordinate system. , , Representing the concentration function The partial derivatives along the X, Y, and Z axes are then defined as follows: The spatial gradient vector is defined as... Its physical meaning is the direction vector of the fastest increase in heavy metal concentration at that spatial point; Simultaneously extract the geological entity model features corresponding to the spatial location of this voxel unit: This represents the implicit surface equation of a specific formation physical interface generated by spatial interpolation fitting of borehole formation elevation data; This represents the surface normal vector at that spatial point, representing the physical interface of a specific stratum. , , Describe the formation interface equations respectively The partial derivatives along the X, Y, and Z axes are then defined as follows: The surface normal vector is defined as... Its physical meaning is the direction perpendicular to the geological sequence boundary; Introducing a correlation calculation model: Defined as the correlation alignment coefficient, it characterizes the parallelism between the concentration diffusion direction and the formation interface normal. Its underlying formula is the absolute value of the cosine of the angle between the gradient vector and the normal vector, i.e.: ; Defined as the spatial correlation threshold, the range of values ​​for the preset dimensionless empirical constant is set as follows. to between; Extract all nodes within a high-concentration voxel region and calculate their correlation alignment coefficient. The average value; When the average value satisfies At that time, the computer logic determines that the direction of the heavy metal concentration gradient is highly coincident with the normal direction of the formation interface, thereby determining the objective spatial relationship that the pollutant is physically infiltrating along a specific formation interface. If high-concentration areas extend along specific lithological interfaces or are concentrated near jointed fracture zones of weathered argillaceous siltstone, then corresponding spatial distribution slice views and coordinate trajectory records are output, thereby revealing the spatial evolution and distribution patterns of different heavy metal pollutants physically intruding downwards along stratigraphic interfaces and fracture channels. Figure 15 As shown, the output shows a layered distribution slice view along different strata interfaces such as miscellaneous fill, clay, silty clay, and various weathered argillaceous siltstones.

[0032] like Figure 2 As shown, a three-dimensional spatial interpolation system for soil pollution based on adaptive grid encryption includes: a data acquisition module 1, an initial framework construction module 2, an adaptive misalignment elimination module 3, a three-dimensional spatial interpolation module 4, a dynamic closed-loop supplementary investigation module 5, and a three-dimensional spatial modeling module 6. Data acquisition module 1 is used to perform gridded data collection and acquire multi-source environmental survey data containing actual sampling point coordinates and stratigraphic information; Initial framework construction module 2 is used to determine the convex polygon boundary of the target site using the Graham scan method or convex hull mesh generation method, and to construct an initial overall planar mesh framework covering the target site using the Delaunay triangulation algorithm or regular rectangular mesh generation algorithm. The adaptive misalignment elimination module 3 is used to calculate the spatial misalignment distance between the computation nodes and the actual sampling points of the initial overall planar mesh frame. It introduces adaptive mesh refinement technology to dynamically increase the mesh subdivision density in areas where the spatial misalignment distance is greater than a preset threshold until the generated computation nodes and the actual sampling points match in spatial position, thereby obtaining the optimized computation mesh. The three-dimensional spatial interpolation module 4 is used to perform spatial estimation by using the optimized computational grid as the computational carrier and combining the three-dimensional spatial interpolation algorithm with the heavy metal concentration information, so as to obtain the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimation confidence. The dynamic closed-loop supplementary exploration module 5 is used to determine whether the current estimation result meets the accuracy requirements based on the estimation confidence level. If it does not meet the requirements, the substandard area is extracted as the target area to guide the field to carry out dynamic supplementary exploration to obtain new sampling point data. The new sampling point data is then integrated into the multi-source environmental survey data and closed-loop iteration is performed. The 3D spatial modeling module 6 is used to construct a 3D spatial model of soil heavy metal pollution based on concentration distribution data after the accuracy meets the standard. It provides data support for subsequent soil classification, removal and remediation strategies by accurately quantifying the volume of polluted soil, analyzing the complex overlapping relationship between different heavy metal pollutants and assessing the physical intrusion characteristics of pollutants.

[0033] Finally, it should be noted that the above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A three-dimensional spatial interpolation method for soil pollution based on adaptive grid densification, characterized in that, include: Step S1: Grid-based data collection is performed on the target site to obtain multi-source environmental survey data. The multi-source environmental survey data includes at least the three-dimensional spatial coordinates of the actual sampling points, heavy metal concentration information, and stratigraphic information. Step S2: Based on the spatial distribution of actual sampling points, the convex polygon boundary of the target site is determined by Graham scanning method or convex hull mesh generation method, and the initial overall planar mesh framework covering the target site is constructed by Delaunay triangulation algorithm or regular rectangular mesh generation algorithm. Step S3: Calculate the spatial misalignment distance between the computation nodes and actual sampling points of the initial overall planar mesh frame. Introduce adaptive mesh refinement technology to dynamically increase the mesh subdivision density in areas where the spatial misalignment distance is greater than a preset threshold until the generated computation nodes and actual sampling points match in spatial position to obtain the optimized computation mesh. Step S4: Using the optimized computational grid as the computational carrier, a three-dimensional spatial interpolation algorithm is used in conjunction with heavy metal concentration information to perform spatial estimation, thereby obtaining the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimation confidence level. Step S5: Determine whether the current estimation result meets the accuracy requirements based on the estimation confidence level. If not, extract the substandard area as the target area to guide the field to conduct dynamic supplementary exploration to obtain new sampling point data. Then, integrate the new sampling point data into the multi-source environmental survey data and return to iteratively execute step S2. Step S6: If the accuracy requirements are met, a three-dimensional spatial model of soil heavy metal pollution is constructed based on the concentration distribution data. By accurately quantifying the volume of polluted soil, analyzing the complex overlapping relationship between different heavy metal pollutants, and assessing the physical intrusion characteristics of pollutants, data support is provided for subsequent soil classification, removal, and remediation strategies.

2. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, In step S1, data is collected using a grid layout strategy with different densities based on the functional attributes of different historical plots within the target site and their corresponding pollution risk levels. The specific data collection process includes: The target site is divided into multiple sub-areas according to different types. Grid points are set up for different areas to obtain a sampling point distribution map. For core production areas with high pollution risk, high-density grid point parameters are configured, while for peripheral areas or general storage areas with low pollution risk, low-density grid point parameters are configured. After determining the grid layout parameters for each region, the longitude, latitude, and elevation data of each layout point are obtained through the Global Positioning System or real-time dynamic differential positioning technology to establish a three-dimensional spatial coordinate system for the actual sampling points. At each deployment point, undisturbed soil samples were collected vertically using geological drilling equipment. During the drilling process, soil samples were extracted at preset depth intervals, and the geological characteristics of the borehole cores were recorded in detail, including but not limited to soil color, texture, moisture, density, and the boundary depth of different soil and rock layers. By structuring the undisturbed soil sample data obtained from the drilling, three-dimensional geological sequence data containing the deep structure of each sampling point was generated. The collected undisturbed soil samples were sealed on-site and sent to the laboratory for physicochemical analysis. The heavy metal content in the soil samples was determined by inductively coupled plasma mass spectrometry or atomic absorption spectrometry. The obtained heavy metal concentration information includes at least the measured mass fraction of one or more heavy metal elements among arsenic, lead, cadmium, and mercury. The above-mentioned three-dimensional spatial coordinates, heavy metal concentration information, and stratigraphic information containing three-dimensional geological sequence data are formatted and integrated, and stored in a spatial database to form a multi-source environmental survey data set.

3. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, In step S2, constructing an initial overall planar grid framework covering the target site specifically includes: After acquiring multi-source environmental survey data, the horizontal projection coordinates of all actual sampling points were extracted. ; Using this point set as input, the convex polygon boundary of the target site is determined by Graham scan method or convex hull mesh generation method; After determining the boundary of the convex polygon, the Delaunay triangulation algorithm or the regular rectangular mesh partitioning algorithm is used inside the boundary to generate an initial overall planar mesh framework composed of multiple mesh edges and multiple mesh nodes.

4. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, In step S3, obtaining the optimized computational grid specifically includes: Extracting the set of computational nodes from the initial global planar mesh framework and the actual set of sampling points ,against Each computing node in ,exist Search for the nearest actual sampling point And the spatial misalignment distance between the two is calculated using the Euclidean distance formula: ； In the formula, , Representing compute nodes The X-axis and Y-axis coordinates on the horizontal projection plane, , These represent the actual sampling points. X-axis and Y-axis coordinates on the horizontal projection plane; The maximum spatial misalignment distance dynamically varies with the site conditions and grid density of different projects, with a preset threshold. The maximum allowable spatial misalignment error is adaptively set based on the current grid density, with a preset threshold. The functional variables, dynamically calculated based on the grid partitioning reference side length and sampling point density of the target site, are derived and mathematically defined as follows: Set the reference side length for mesh generation as The spatial distribution density parameter of the actual sampling points in the set area is: Define the spatial tolerance coefficient as Then the preset threshold The calculation formula is: ； In the formula, The damping coefficient is fitted using historical survey data, and its value ranges from [value range missing]. Spatial tolerance coefficient Depending on the sensitivity of the interpolation algorithm used to coordinate offset, its empirical range is strictly limited to [specific values]. When the local encryption process reduces the mesh edge length Reduce or local sampling density When it increases, the preset threshold It exhibits exponential adaptive decay, thereby ensuring the uniqueness and determinism of the convergence boundary of the adaptive encryption process at the underlying logic level; Traverse all computing nodes, and when a spatial misalignment distance between a computing node and its corresponding actual sampling point is detected... When it is determined that the local area needs adaptive mesh refinement, the specific process of adaptive refinement is as follows: Add new subdivision nodes inside the grid cell containing the misaligned calculation node. That is, use the quadtree recursive subdivision method to divide a rectangular grid cell with misalignment into four smaller sub-cells and recalculate the distance between the sub-cell nodes and the actual sampling points. By locally and recursively splitting the grid cells and densifying the nodes, the coordinates of the newly generated computational nodes gradually approximate the coordinates of the actual sampling points. This adaptive encryption process iterates until the spatial misalignment distance between all computing nodes and actual sampling points is less than or equal to a preset threshold. At this point, the mesh cell subdivision is stopped, the current mesh topology is output, and it is used as the optimized horizontal computational mesh. The scheme of introducing three-dimensional geological sequence data to divide the computational grid nodes involves, after outputting the optimized horizontal computational grid, calling the three-dimensional geological sequence data and, based on the elevations of various strata recorded in the geological borehole data, dividing the optimized horizontal computational grid nodes into layers in the vertical direction. The coordinates of the 2D grid nodes Substitute the spatial surface equations of each fitted stratigraphic interface In the process, the values ​​of each node at different stratigraphic interfaces were calculated. Axis elevation values, among which, A three-dimensional surface relationship function for a specific stratigraphic interface is generated by spatial interpolation fitting based on the stratigraphic elevations recorded in geological borehole data. Based on the above elevation values, the two-dimensional planar grid is stretched and divided into multiple layers of three-dimensional grid units, so that the nodes of the divided three-dimensional computational grid are attached to the actual undulation of the strata, forming a three-dimensional computational grid base containing X, Y, and Z three-dimensional coordinate attributes. The three-dimensional computational grid base containing X, Y, and Z coordinate attributes is the final, complete, and optimized computational grid.

5. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, In step S4, obtaining the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimated confidence level specifically includes: After constructing the optimized computational grid, the measured heavy metal concentration data from the actual sampling points are mapped to the corresponding computational grid nodes. A three-dimensional Kriging interpolation method is used, comprehensively considering the spatial distribution of the actual sampling point data in both the vertical and horizontal directions. Based on the covariance function, a mathematical model of the random field is performed to obtain an unbiased estimate of the heavy metal concentration, i.e.: Based on the known three-dimensional coordinates and concentration values ​​of the sampling points, the experimental semivariogram function is calculated to describe the degree of variation of the variable with increasing spatial distance. The actual sampling point set is paired and grouped according to different spatial distance steps, and statistical calculations are performed on this basis. The specific physical meaning of the experimental semivariogram function is half of the variance of the heavy metal concentration difference between two sampling points spatially separated by a certain distance. Its specific mathematical formula is as follows: ； In the formula, The experimental semivariance function, The step size is the lag in the three-dimensional spatial distance between two sampling points. Spatial distance equal to The number of actual sampled data pairs, This is the index of the actual sampled data pair. , It is a three-dimensional spatial coordinate vector that includes horizontal, vertical, and elevation dimensions. Located at a three-dimensional coordinate point The measured heavy metal concentration at the location, Located at the deviation point Distance is The measured heavy metal concentration at the location; A suitable theoretical semivariance model is selected, and the weighted least squares method is used to fit the experimental semivariance data to obtain the nugget value, sill value, and range parameter. Based on these parameters, a covariance function is constructed to mathematically model the three-dimensional random field. The formula for calculating the covariance function is as follows: ； In the formula, Let covariance function be used. To fit the obtained sill values, which are the theoretical maximum estimated variances of all nodes within the target site, This is the theoretical semivariance value; For the unknown node to be estimated, based on its three-dimensional spatial distance from surrounding known sampling points and the constructed covariance function, a system of Kriging linear equations is established: ； In the formula, The total number of known sampling points involved in the estimation calculation. Label the known sampling points used in the estimation calculation. , Let be the three-dimensional position vector of the unknown node to be estimated. For the first Three-dimensional position vectors of known sampling points The first one to be solved Known sampling points for unknown estimation points The weighting coefficients, The known points are calculated using the semivariance function described above. and The semivariance between them For known points With the point to be estimated The semivariance between them The introduced Lagrange multipliers are used to ensure the unbiased condition that the sum of the weight coefficients is 1; by performing matrix inversion on the above linear equations, all weight coefficients can be analytically obtained. ; By solving this system of equations, the weight coefficients of each known sampling point are obtained, ensuring that the sum of these weight coefficients is 1 and the estimation variance is minimized, thus achieving an unbiased estimation of heavy metal concentration. After solving the Kriging linear equations to obtain the weighting coefficients, the heavy metal concentration of the node to be estimated is calculated as a weighted sum of the measured concentrations at surrounding known sampling points, using the following formula: ； In the formula, For the node to be estimated The estimated concentration at that location, For the first Measured concentration values ​​of heavy metals at known sampling points; And, calculate the kriging estimate variance for each grid node, using the following formula: ； In the formula, For known points With the point to be estimated The semivariance between them; Define the theoretical maximum estimated variance of all nodes within the target site as: This value is taken from the sill value parameter obtained when fitting the theoretical semivariance model in step S4. ; For any grid node in space ,in These represent the horizontal, vertical, and elevation coordinates of the 3D spatial grid node in the 3D coordinate system, respectively, and their estimated confidence levels. The calculation formula is as follows: ； In the formula, For grid nodes Kriging estimates the variance at the location.

6. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, Step S5 specifically includes: Set the global spatial interpolation accuracy threshold The range of values ​​is When any node's When this happens, the dynamic supplementary exploration logic is triggered, that is: Iterate through the estimated confidence data of all 3D computational mesh nodes. If a mesh node exists... If so, the current estimation result is determined to not meet the accuracy requirements; The system uses a spatial clustering algorithm to merge adjacent low-confidence grid nodes and extract a set of substandard 3D spatial regions as target areas for further exploration. For the extracted target area to be explored, the system executes the logic of generating target exploration points, calculates the confidence gradient center or spatial geometric center within the target area, and generates guiding coordinate instructions containing the target longitude, latitude and drilling depth. in, The confidence gradient center within the target area is calculated using a spatial weighted calculation with the difference between the confidence level of each grid node and its attainment threshold as the weighting coefficient. The underlying mathematical derivation and parameter definitions are as follows: The extracted substandard 3D spatial target area contains a total of For any given grid node within the target area, 1 node It is the index number, and Its parameters are defined as follows: Indicates the first The three-dimensional coordinates of each grid node in the real physical space correspond to the X-axis coordinate, Y-axis coordinate, and Z-axis depth in the projection plane, respectively. This represents the first variance calculated using Kriging variance. The estimated confidence level of each grid node; This represents the system's preset global spatial interpolation accuracy threshold. Indicates the first The missing confidence values ​​for each node are defined as follows: The larger this value, the higher the uncertainty of the data at that coordinate point, and the greater its weight in the optimization calculation; Based on the above parameters, the derivation formula for the three-dimensional coordinates of the confidence gradient center of the target area is as follows: ； ； ； In the formula, , and These represent the absolute physical coordinates of the target re-exploration point obtained after weighting by confidence error; The coordinates calculated by the above formula will naturally be biased towards the spatial area with the largest estimation error, and these will be used as priority target points for resurvey, generating guiding instructions that include these coordinates. The instruction is sent to the terminal interface of the on-site exploration equipment via the communication network. On-site personnel conduct supplementary drilling and sampling in the target area according to the coordinate instructions, test and analyze the newly added undisturbed soil samples, and obtain data from the newly added sampling points. The newly added sampling point data is added to the spatial database and merged with the original multi-source environmental survey data. The updated data set is then used to re-trigger the execution of steps S2-S4. This dynamic closed-loop process iterates continuously until all 3D computational mesh nodes are reached. All are greater than or equal to Only when the accuracy requirement is met can the iteration be terminated and the process proceed to step S6.

7. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, In step S6, after the three-dimensional interpolation accuracy meets the standard, a spatial entity is generated in the three-dimensional geological modeling environment using the spatial interpolation calculation results. Specifically, this includes: Import the digital elevation model data of the target site and generate a three-dimensional surface of the ground using an irregular triangular mesh. Extract the various types of layered stratigraphic interface data processed in step S3, and similarly use irregular triangular mesh to construct three-dimensional undulating surfaces of each stratigraphic interface, such as miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. Using the three-dimensional curved surface of the earth's surface as the top plate and the bottom geological interface as the bottom plate, the geological interfaces in the middle of each layer are used as internal dividing surfaces. Through surface stretching and solid filling calculations along the Z-axis, a multi-layered three-dimensional geological solid model is constructed. A three-dimensional spatial model of soil heavy metal pollution was constructed by combining the heavy metal concentration interpolation results, namely: The three-dimensional concentration distribution data of each node output in step S4 are converted into a voxel data structure according to the spatial coordinate mapping rules. A preset pollution concentration threshold is set, and the cells in the voxel model whose concentration exceeds the threshold are defined as pollution voxels. The voxel model consists of voxels with a spatial resolution of [missing information]. Composed of a rectangular grid array, , , These represent the physical side lengths of a single voxel unit along the X, Y, and Z axes in a three-dimensional spatial model, respectively. Define three-dimensional spatial logic state functions ; For the specific heavy metal arsenic, a risk control standard value is set at [value missing]. Extract the three-dimensional concentration distribution data of the nodes. The entity mapping logic judgment condition is expressed as follows: when At that time, assign a state value to the voxel unit. , defined as effective voxels that are contaminated; when At that time, assign a state value to the voxel unit. Defined as a pollution-free blank voxel; The underlying operation rule of Boolean intersection is the logical AND operation of corresponding matrix elements, and the exact relation is: ； In the formula, , , , Representing the coordinates in three-dimensional space respectively Logical state values ​​of single heavy metal pollutants (arsenic, lead, cadmium, and mercury) at the location; This is the Boolean logical AND operator; Only when all single-factor state values ​​at a certain spatial voxel coordinate are equal to 1 The output value is 1, thus accurately extracting the three-dimensional boundary coordinate point cloud set of the complex heavy metal pollution area; For other heavy metal factors, namely lead, cadmium, and mercury, extract their corresponding node three-dimensional concentration distribution data. , , And set corresponding risk control standard values. , , Based on the entity mapping logic judgment conditions that are completely consistent with those for arsenic, voxel unit state values ​​corresponding to lead, cadmium, and mercury are generated. , , ; These pollutant voxels are rendered into a three-dimensional geological entity model to form a three-dimensional spatial model of soil heavy metal pollution. Subsequently, using the Boolean logic operation module, including union, intersection, and difference operations, spatial entity calculations are performed on the constructed model to define the contaminated boundary volume. Heavy metals Spatial entities whose concentration exceeds their corresponding risk threshold, among which, To characterize the set subscripts for different heavy metal factors, ,in, Intersection operation is used to extract three-dimensional spatial regions that simultaneously meet the standards for multiple heavy metals exceeding limits, obtaining the three-dimensional boundary and coordinate set of the complex heavy metal pollution region: ； In the formula, A collection of three-dimensional spatial entities representing a region contaminated with complex heavy metals; The difference operation is used to extract the boundary of a spatial region where only a single heavy metal exceeds the standard: ； In the formula, This represents the spatial boundary of a region where the single heavy metal arsenic exceeds the standard. , , , Let each represent a set of spatial entities whose concentrations of arsenic, lead, cadmium, and mercury exceed their respective risk thresholds. The set intersection operator. The set union operator. This is the set difference operator; Based on the single heavy metal pollution area and the complex heavy metal pollution area obtained above, the physical volume of the monomeric element unit is set as follows: Total earthwork volume in areas contaminated with complex heavy metals The calculation formula is: ； In the formula, For three-dimensional coordinates The logical state value of the complex heavy metal pollutant at the location. The physical volume of a monomeric unit.

8. The method for three-dimensional spatial interpolation of soil pollution based on adaptive grid encryption according to claim 1, characterized in that, In step S6, the three-dimensional spatial model of soil heavy metal pollution representing the pollution concentration range and the three-dimensional geological entity model representing the physical barrier characteristics are registered and superimposed in the same three-dimensional coordinate system, and the spatial gradient vector in the three-dimensional concentration field is extracted to analyze the normal direction of concentration as the coordinate changes. By comparing the distribution location of high-concentration gradient vectors with the stratigraphic interfaces and weathering layer structure in the 3D geological entity model, the spatial correlation data between the two is output, namely: The logic for extracting and determining spatial location correlation data is based on the inner product operation of the concentration field gradient vector and the geological interface normal vector. The specific parameter definitions and derivation relationships are as follows: This represents a three-dimensional continuous heavy metal concentration space function generated based on the optimized computational grid fitting. This represents the spatial gradient vector of the concentration field of a high-concentration voxel unit in this spatial coordinate system. , , Representing the concentration function The partial derivatives along the X, Y, and Z axes are then defined as follows: The spatial gradient vector is defined as... Its physical meaning is the direction vector of the fastest increase in heavy metal concentration at that spatial point; Simultaneously extract the geological entity model features corresponding to the spatial location of this voxel unit: This represents the implicit surface equation of a specific formation physical interface generated by spatial interpolation fitting of borehole formation elevation data; This represents the surface normal vector at that spatial point, representing the physical interface of a specific stratum. , , Describe the formation interface equations respectively The partial derivatives along the X, Y, and Z axes are then defined as follows: The surface normal vector is defined as... Its physical meaning is the direction perpendicular to the geological sequence boundary; Introducing a correlation calculation model: Defined as the correlation alignment coefficient, it characterizes the parallelism between the concentration diffusion direction and the formation interface normal. Its underlying formula is the absolute value of the cosine of the angle between the gradient vector and the normal vector, i.e.: ； Defined as the spatial correlation threshold, the range of values ​​for the preset dimensionless empirical constant is set as follows. to between; Extract all nodes within a high-concentration voxel region and calculate their correlation alignment coefficient. The average value; When the average value satisfies At that time, the computer logic determines that the direction of the heavy metal concentration gradient is highly coincident with the normal direction of the formation interface, thereby determining the objective spatial relationship that the pollutant is physically infiltrating along a specific formation interface. If the high-concentration area extends along a specific lithological interface in terms of spatial morphology, or is concentrated near the joint and fracture zone of weathered silty mudstone, then the corresponding spatial distribution slice view and coordinate trajectory record are output, thereby revealing the spatial evolution and distribution law of different heavy metal pollutants physically invading downward along the stratigraphic interface and fracture channel.

9. A three-dimensional spatial interpolation system for soil pollution based on adaptive grid encryption, according to any one of claims 1-8, characterized in that, include: The data acquisition module is used to perform gridded data collection and acquire multi-source environmental survey data containing actual sampling point coordinates and stratigraphic information; The initial framework construction module is used to determine the convex polygon boundaries of the target site using the Graham scan method or the convex hull mesh generation method, and to construct an initial global planar mesh framework covering the target site using the Delaunay triangulation algorithm or the regular rectangular mesh generation algorithm. The adaptive misalignment elimination module is used to calculate the spatial misalignment distance between the computation nodes and the actual sampling points of the initial overall planar mesh frame. It introduces adaptive mesh refinement technology to dynamically increase the mesh subdivision density in areas where the spatial misalignment distance is greater than a preset threshold until the generated computation nodes and the actual sampling points match in spatial position, thus obtaining the optimized computation mesh. The three-dimensional spatial interpolation module is used to perform spatial estimation by combining the three-dimensional spatial interpolation algorithm with heavy metal concentration information, using the optimized computational grid as the computational carrier, to obtain the concentration distribution data of each heavy metal in three-dimensional space and its corresponding estimation confidence level. The dynamic closed-loop supplementary exploration module is used to determine whether the current estimation result meets the accuracy requirements based on the estimation confidence level. If it does not meet the requirements, the substandard area is extracted as the target area to guide the field to conduct dynamic supplementary exploration to obtain new sampling point data. The new sampling point data is then integrated into the multi-source environmental survey data and closed-loop iteration is performed. The 3D spatial modeling module is used to construct a 3D spatial model of soil heavy metal pollution based on concentration distribution data after the accuracy standard is met. It provides data support for subsequent soil classification, removal and remediation strategies by accurately quantifying the volume of polluted soil, analyzing the complex overlapping relationship between different heavy metal pollutants and assessing the physical intrusion characteristics of pollutants.

Images