Method and system for soil pollution simulation and migration path evaluation coupled with three-dimensional stratum

Abstract

The application discloses a kind of three-dimensional stratum coupling soil pollution simulation and migration path evaluation method, belong to computer data processing technical field, comprising: step S1, obtain geological exploration and heavy metal sampling data;Step S2, generate each stratum top and bottom plate three-dimensional grid surface, establish the physical space boundary of each heterogeneous stratum, form three-dimensional stratum framework;Step S3, construct three-dimensional geological entity model;Step S4, generate layered pollution grid in geometric shape and space attribute and align with three-dimensional stratum framework;Step S5, pollution grid and entity model depth coupling, analyze the absolute pollution distribution state of multiple heavy metals in different specific stratum;Step S6, quantitatively statistics pollution volume and assess the differential migration law of heavy metal along specific channel in combination with stratum property.The application effectively blocks the interpolation smoothing effect across heterogeneous stratum, so that the pollution simulation result objectively reflects the geological control effect.

Description

Technical Field

[0001] This invention belongs to the field of computer data processing technology, specifically a method and system for three-dimensional stratigraphic coupling of soil pollution simulation and migration path assessment. Background Technology

[0002] In the interdisciplinary field of soil environmental science and geological engineering, the spatial distribution of heavy metal pollution in soil typically exhibits non-uniform characteristics. Current soil pollution status assessment techniques often rely on field survey data, combined with two-dimensional spatial interpolation algorithms to simulate the planar distribution of a single heavy metal element at different sampling depths, or employ conventional three-dimensional spatial interpolation algorithms to perform isotropic concentration estimation across the entire computational domain.

[0003] The vertical distribution and migration of heavy metals in soil and underground media are influenced by physical stratigraphic structure. Different geological bodies, such as artificial fill layers, low-permeability clay zones, and weathered rock layers with developed fissures, exert different physical control effects on the migration, diffusion, and distribution of pollutants.

[0004] In existing technical solutions, the spatial interpolation calculation of heavy metal concentration is usually relatively independent of the three-dimensional physical stratigraphic model. The distance or covariance calculation of the interpolation algorithm does not take the physical boundaries of the strata as constraints, which may lead to the generation of pollution distribution models exhibiting calculation smoothing phenomena across non-permeable strata, making it difficult to accurately reflect the control of the strata on the spatial distribution and migration patterns of heavy metals.

[0005] Therefore, to address the above issues, a method and system for three-dimensional stratigraphic coupling of soil pollution simulation and migration path assessment is provided. Summary of the Invention

[0006] To address the aforementioned problems in the existing technology, this invention provides a method and system for three-dimensional stratigraphic coupling of soil pollution simulation and migration path assessment. Through a controlled physical spatial boundary mechanism, it objectively analyzes the three-dimensional spatial distribution of pollutants within different geological layers and assesses their corresponding migration patterns. This effectively blocks the interpolation smoothing effect across heterogeneous strata, enabling the pollution simulation results to objectively reflect the geological control effect.

[0007] The technical solution to achieve the above objectives is: One aspect of this invention is a three-dimensional stratigraphic coupled method for soil pollution simulation and migration path assessment, comprising: Step S1: Obtain geological survey data and heavy metal sampling data for the study area. After preprocessing, the strata in the study area are generalized into six strata distributed from the surface downwards: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. Among them, the geological survey data includes borehole spatial coordinates, borehole columnar section data, borehole stratification data, strata lithology classification information, groundwater level depth data, and physical, mechanical, and hydrogeological parameters of each rock and soil layer. The heavy metal sampling data includes measured concentration data of arsenic, lead, cadmium, and mercury at different depths. Step S2: Based on the borehole spatial coordinates, an adaptive planar mesh is constructed using the Delaunay triangulation algorithm or the Voronoi diagram mesh generation algorithm. The planar mesh is stretched along the Z-axis by combining digital elevation model data and borehole layer data. Topological anti-intersection and node merging logic is configured to handle the formation pinch-out phenomenon, generate three-dimensional mesh surfaces of the top and bottom plates of each stratum, establish the physical spatial boundaries of each heterogeneous stratum, and form a three-dimensional stratigraphic framework. Step S3: Based on the interface mesh of the three-dimensional stratigraphic framework, the three-dimensional boundary shell of the study area is generated by the convex hull mesh algorithm, the undulating area of ​​the stratigraphic interface is adaptively densified, and the solid filling between the curved surfaces is realized by the three-dimensional mesh solidification algorithm to construct a three-dimensional geological solid model. Step S4: When performing spatial interpolation calculations on heavy metal sampling data to establish a three-dimensional distribution model of soil heavy metal pollution, the layered grid of the above-mentioned three-dimensional stratigraphic framework is forcibly used as the constraint condition for spatial interpolation, so that the interpolation calculation of heavy metal concentration is controlled by the physical spatial boundary, thereby generating a layered pollution grid that is aligned with the three-dimensional stratigraphic framework in terms of geometric shape and spatial properties. Step S5: Spatial coordinate matching is performed between the layered pollution grid and the three-dimensional geological entity model. The heavy metal concentration attribute and the stratum attribute are jointly mapped by the octree index and the axial bounding box collision detection algorithm to form a three-dimensional composite entity model. Based on the stratum attribute and concentration threshold filtering conditions, the absolute distribution of heavy metal pollution in different strata is refined and analyzed. Step S6: Based on the absolute pollution distribution status, quantitatively calculate the pollution volume of each heavy metal in different strata, and combine the strata physical properties to assess the differentiated migration patterns of heavy metals along specific channels and the potential risk of deep pollution.

[0008] Preferably, in step S1, the preprocessing process includes data cleaning, outlier removal, and coordinate unification; spatial coordinate data from different sources and formats are uniformly converted to a preset standard three-dimensional spatial coordinate system to establish a unified spatial calculation benchmark. After data preprocessing, based on borehole columnar section data and stratigraphic lithology classification information, the strata in the study area were generalized into six strata distributed from the surface downwards: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. The generalized strata described above represent the main media units with different permeability and physical barrier properties within the study area, among which, The fill layer is configured as the main infiltration and storage layer for pollutants; Clay and silty clay layers, due to their low permeability coefficient, are configured as barrier layers or slow-seepage layers for pollutants. Weathered silty mudstone layers are configured as potential deep migration channels for pollutants, depending on their degree of weathering and the development of fissures.

[0009] Preferably, in step S2, during the adaptive planar mesh generation process, the algorithm dynamically adjusts the mesh resolution based on the spatial density distribution of borehole data points: In areas with dense borehole distribution or large gradients in formation elevation, the grid side length is automatically reduced to increase node density; in areas with sparse borehole distribution or gentle changes in formation elevation, the grid size is increased accordingly. The side length adjustment logic for adaptive planar meshes is implemented through a preset mesh size control function: Let the side length of the basic grid in the study area be... For any evaluation node in the grid generation algorithm, calculate its preset search radius. Total number of boreholes And the absolute value of the rate of change of the elevation gradient at the location of the node. Then the grid side length The specific mathematical relation configuration is as follows: ; ; In the formula, Represented by natural constant An exponential function with base 0. This is the borehole density weighting coefficient, and its value range is configured as follows: , The elevation gradient weighting coefficient has a value range configured as follows: , The gradient vector of elevation in a horizontal two-dimensional plane. For partial derivative mathematical operations, and These are the x and y components of the horizontal two-dimensional plane, respectively. This refers to the elevation value of the corresponding grid node; When the calculated grid side length Less than the preset minimum grid limit threshold Forced configuration This is to control the upper limit of memory consumption in the computing domain; Further, digital elevation model data of the study area is obtained, and the digital elevation model data is mapped to the nodes of the adaptive planar grid, and the corresponding initial Z-axis coordinates of the ground surface elevation are assigned to the two-dimensional grid nodes. Extract the borehole stratification data of each generalized stratum recorded in the geological exploration borehole data, namely the burial depth data of the top plate and bottom plate; Based on the adaptive planar grid of the ground elevation, a three-dimensional spatial stretching operation is performed along the vertical downward direction of the Z-axis; During the stretching algorithm, the nodes of the adaptive planar mesh are anchored when they reach the elevation of the top and bottom plate interface of each generalized stratum determined by the borehole layering data, thereby generating a three-dimensional mesh surface representing each stratum interface in three-dimensional space. To address the stratigraphic pinch-out phenomenon present in actual geological exploration data, the 3D stratigraphic framework incorporates topological anti-intersection and node merging logic during the stretching construction process, namely: Let the second be an adaptive planar mesh. Node coordinates At this location, the first in the generalized stratigraphic sequence The elevation scale of the top surface of the layer is The bottom elevation scale is ,in, These correspond to six generalized strata, ranging from miscellaneous fill to completely weathered argillaceous siltstone. is the globally unique positive integer index of the node. , This represents the total number of nodes in the two-dimensional planar mesh. , The stretching calculation unit calculates the X-axis and Y-axis coordinate components of the node in a preset standard three-dimensional coordinate system. Vertical thickness at the location The formula for its calculation is: ; The algorithm pre-configures a threshold for determining the minimum formation thickness. ; When logical judgment conditions When established, it is determined that the first occurrence occurred at that spatial coordinate. The strata pinch out; The algorithm forces the node merging instruction to be executed, making The upcoming The top and bottom grid surfaces of the layer are in The topological nodes at the coordinates are forced to coincide in the Z-axis direction; By using this underlying constraint logic, abnormal geometric units with zero or negative thickness are eliminated, thereby generating a three-dimensional stratigraphic framework that satisfies the geometric conditions of a closed manifold without self-intersection.

[0010] Preferably, in step S3, constructing the three-dimensional geological entity model includes: First, the convex hull mesh algorithm is used to determine the spatial envelope of the outermost data points in order to generate the three-dimensional boundary shell of the study area; Secondly, adaptive mesh refinement technology is applied to areas with large undulations or intersections at the stratigraphic interfaces. By increasing the number of polygons on the local mesh surfaces, the geometric shape of the geological interfaces is smoothly expressed, thereby generating the generalized three-dimensional surfaces of the interfaces of the six stratigraphic layers. Finally, using the surface surface generated from the digital elevation model data as the top interface and the three-dimensional surfaces of each internal stratum as the middle and bottom interfaces, the closure and solid filling between adjacent surfaces are achieved through the three-dimensional mesh solidification algorithm, thus constructing a three-dimensional geological solid model. Among them, the three-dimensional mesh solidification algorithm is either a voxelization algorithm or a solid filling algorithm based on the boundary representation method B-Rep; In the three-dimensional geological entity model, the space is divided into discrete voxel units or unstructured tetrahedral / hexahedral mesh units, and each three-dimensional geometric unit is configured and assigned a corresponding stratigraphic attribute code.

[0011] Preferably, in step S4, the heavy metal sampling data includes the spatial coordinates of heavy metals such as arsenic, lead, cadmium, and mercury. and its measured concentration scalar value; A three-dimensional kriging interpolation algorithm is used to estimate concentration using three-dimensional mesh nodes as the computational carrier. The algorithm internally incorporates a spherical semivariogram function model to quantify spatial autocorrelation. The mathematical definition rules are as follows: when hour: ; when hour: ; In the formula, Let Euclidean distance be the three-dimensional straight-line distance between the grid node to be estimated and the known heavy metal sampling point. To quantify the semivariance function value of spatial autocorrelation, The nugget constant represents the spatial discontinuity caused by measurement errors and microscale variations. Its value is obtained by extrapolating the semivariance of the full sample data within the minimum sampling interval using a linear fit. Obtaining the intercept at that point, For the partial sill value, The sill value constitutes the total sample space variance of the entire interpolation calculation. The range threshold represents the maximum effective distance boundary of spatial correlation. It is determined as follows: A scatter plot of the experimental variogram of the measured samples is plotted; a spherical curve is fitted using the nonlinear least squares method; and the x-axis distance value corresponding to when the fitted curve reaches the sill level is extracted as the range threshold. ; After completing the constrained space interpolation calculation, the generated output is a hierarchical pollution grid, which contains a grid structure composed of three-dimensional nodes. Each node stores the concentration values ​​of various heavy metals obtained from the interpolation calculation. Its interpolation process is strictly controlled by the layered grid interface of the three-dimensional stratigraphic framework. The layered contamination grid is consistent and aligned with the three-dimensional stratigraphic framework in terms of geometric boundaries and spatial topology.

[0012] Preferably, step S5 includes: The process involves traversing all nodes or voxels in the layered pollution grid, extracting the heavy metal concentration scalar information they carry, and mapping this concentration information to the corresponding coordinate units in the three-dimensional geological entity model according to the spatial coordinate matching principle. Each smallest computational unit in the three-dimensional composite entity model after mapping contains three-dimensional spatial coordinate data, stratigraphic lithology data, and concentration attribute data of multiple heavy metals. Based on the three-dimensional composite entity model after mapping, the system is configured with spatial attribute query and Boolean operation logic. By setting the stratum attribute filtering conditions and concentration threshold filtering conditions, the spatial distribution of pollutants inside a specific stratum can be analyzed and extracted respectively. Among them, the efficiency of mapping and matching operations in the three-dimensional coordinate system is improved by using octree indexing and axis bounding box collision detection algorithms. The specific mapping extraction steps are as follows: First, using the global spatial cube boundary of the 3D geological entity model as the root node, recursively construct a maximum depth of... The octree index assigns all voxel units carrying formation attributes to the corresponding octree leaf nodes. Secondly, for any target node in the layered pollution grid, let its three-dimensional spatial coordinates be... The system generates this node. With geometric center and side length as The bounding box of the cube ,in, The edge length is the globally unique positive integer index of the target node in the contaminated mesh. Forced configuration to the minimum side length of the reference voxel in the 3D geological solid model times; Finally, by traversing the octree index, the bounding box can be quickly retrieved and queried. The set of leaf nodes that exhibit spatial geometric overlap is identified, and a 3D point containment test algorithm is invoked to determine the unique containing node in 3D space. The target voxel unit.

[0013] Preferably, step S6 includes: Quantitative statistics of contamination volume: First, for each heavy metal element and each generalized stratum, the system receives preset environmental risk control concentration limits as input parameters; Secondly, the statistical algorithm traverses all voxel units within the corresponding stratum, compares the heavy metal concentration attribute values ​​of the voxel units with the set concentration limits, and filters out all target voxel units whose concentration attribute values ​​are greater than or equal to the set limits. Finally, by summing and integrating the geometric volumes of all selected target voxel units, the total volume of contaminated soil for each heavy metal in different generalized strata, including fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone, is quantitatively calculated and output. Let the nominal geometric volume of a single reference voxel be... The scalar concentration of a certain heavy metal contained within this voxel is [value missing]. The preset environmental risk control concentration limit for this type of heavy metal is Define a dimensionless volume correction factor. ; When the concentration inside the target voxel Furthermore, the concentration values ​​of the voxel in the three-dimensional space of all 26 adjacent neighboring voxels are greater than or equal to the concentration of the voxel. When configuring ; When the internal concentration of voxels However, there exists at least one neighboring voxel with a concentration value less than [a certain value]. When the voxel is located on the three-dimensional boundary of the contamination isosurface, the isosurface function is calculated by calling a three-dimensional linear interpolation algorithm. The effective oversized polyhedron volume formed by truncating the target voxel, and the effective oversized polyhedron volume and the nominal geometric volume are then compared. The ratio is assigned to ; in, For arbitrary spatial coordinates within a 3D computational domain generated based on spatial interpolation The continuous heavy metal concentration distribution function at a given location; Total absolute volume of a single type of heavy metal in contaminated soil within a specific generalized stratum The mathematical calculation relation is configured as follows: ; In the formula, This represents the total absolute volume of contaminated soil containing a single type of heavy metal within a specific generalized stratum. The nominal geometric volume of a single reference voxel unit. For the logical conditions to be met within this specific stratum The total number of voxel units, The iterative index of the effective voxels participating in the accumulation of volume integrals; Quantitative assessment of migration patterns: Extract the continuous set of voxels with excessive concentrations within a specific formation, and calculate the concentration gradient vector using the three-dimensional central difference method. Combined with the permeability coefficient scalar of the corresponding strata Effective porosity parameters Concentration limits for heavy metal risk control Construct a normalized dimensionless concentration gradient Define the apparent migration potential vector controlled by formation permeability. The mathematical relationship is: ; ; In the formula, For heavy metal concentration variables in continuous space. , , These represent the horizontal, vertical, and axial coordinate components of the spatial coordinates within the three-dimensional computational domain. These are mathematical operators for partial derivatives; Principal component analysis (PCA) was used to perform dimensionality reduction on the extracted continuous voxel set's 3D coordinate point cloud data. The eigenvector direction corresponding to the first principal component was then extracted and used as the spatial dominance dip vector for high-pollution areas. ; Read the joint and fracture surface attitude information recorded in the geological exploration data, and generate a unit reference vector representing the dip and dip angle of the joint surface. ; Calculate the cosine similarity between two vectors The calculation formula is as follows: ; In the formula, The spatial dominance dip vector for high-pollution areas. As a unit reference vector representing the dip and dip angle of the joint surface, For the dot product of vectors, and Representing vectors respectively and The modulus length; System configuration similarity threshold When logical judgment conditions When it was established, the distribution direction of the pollutants was determined to be highly coincident with the direction of natural joints and fissures, and the assessment and identification results of heavy metal pollutants deep intrusion along specific joint and fissure channels were output. The specific analytical logic includes: In a layer of miscellaneous fill with high surface permeability, the longitudinal infiltration trend of pollutants driven by hydrodynamics was analyzed. At the interface between miscellaneous fill and low-permeability clay layer, the lateral distribution characteristics of the concentration gradient were extracted to assess the lateral migration and local enrichment mechanisms of pollutants after they were blocked. In fully weathered, strongly weathered, and moderately weathered argillaceous siltstone layers, the three-dimensional geometric morphology of high-pollution areas is analyzed. If the high-pollution areas in the weathered argillaceous siltstone layers are distributed in a linear, vein-like, or finger-like pattern extending along a specific spatial dip angle, and this spatial distribution pattern matches the spatial occurrence of the joint and fissure development zone recorded in the geological exploration data, then the phenomenon of heavy metal pollutants locally intruding into the deep geological environment along the dominant flow migration channel of the weathered layer fissures is identified and confirmed. Generation of 3D targeted repair model: Based on the identification of lateral migration and local enrichment mechanisms, as well as the phenomenon of heavy metal pollutants locally intruding into the deep geological environment along the dominant flow migration channels of weathering layer fissures, the three-dimensional spatial coordinates corresponding to the relevant phenomena are further extracted to generate a three-dimensional targeted remediation spatial control model. This model includes a set of coordinates for shallow excavation barrier boundaries and a set of coordinates for deep in-situ injection remediation nodes. The generation steps include: The first step is to generate a set of shallow excavation barrier boundary coordinates based on lateral migration and local enrichment mechanisms, namely: For the lateral enrichment region at the interface between the fill and the low-permeability clay layer identified in the assessment, traverse the region to satisfy the logical conditions. For all target voxel units, extract the horizontal and vertical coordinate point cloud of the target voxel units in the horizontal two-dimensional plane, and apply the convex hull algorithm to generate the outermost two-dimensional closed polygon surrounding the point cloud. And extend a preset safety buffer distance outward along the normal direction of the polygon boundary. The two-dimensional planar boundary for the repair excavation is generated, where the value range is forcibly configured based on the empirical porosity of the site fill soil. ; Simultaneously, the system iterates through the Z-axis coordinate components in the three-dimensional spatial coordinates of the aforementioned target voxel units, extracting the minimum value as the elevation scalar of the excavation barrier bottom interface. The x and y coordinates of the extended two-dimensional plane boundary and Together, we establish and output the set of three-dimensional coordinates of the shallow excavation barrier boundary; The second step involves generating a set of coordinates for deep in-situ injection remediation nodes based on the phenomenon of heavy metal pollutants locally intruding into the deep geological environment through dominant flow migration channels along weathering layer fissures. For satisfying logical judgment conditions Furthermore, the target voxel set within the weathered rock layer where deep intrusion occurred has been identified, assuming that this set contains... Individual unit; The system calculates and extracts the coordinates of the three-dimensional geometric center of the voxel set. ;in, , , This is the arithmetic mean of the X, Y, and Z coordinate components of all voxel units within the set in a preset standard three-dimensional spatial coordinate system; Simultaneously, the system extracts the maximum value of the Z-axis coordinate component from all voxel elements within the set as the elevation of the highest point. Extract the minimum value of the Z-axis coordinate components as the elevation of the lowest point. ; Subsequently, with Starting from the spatial origin, along the spatial dominance dip vector extracted by the aforementioned dimensionality reduction calculation. The determined spatial straight path, defined in and Within the vertical depth range between them, according to the preset injection node spacing Perform equidistant spatial node sampling; where, The spatial linear Euclidean distance between two adjacent injection points is defined as the physical parameter of connectivity of fractures in weathered rock mass, and its value range is forcibly configured as follows: ; The series of discrete three-dimensional coordinate points generated by the above equidistant sampling constitute the set of coordinates of deep in-situ injection repair nodes, which are used to limit the injection depth and spatial target placement range of the agent in the actual site repair process. Output: Based on quantitative pollution volume statistics, identified lateral migration mechanisms, and spatial location and depth information of deep intrusion phenomena along fissures, a pollution risk distribution data array and assessment report for the study area are output.

[0014] A second aspect of the present invention provides a three-dimensional stratigraphic coupled soil pollution simulation and migration path assessment system, comprising: The multi-source data acquisition module is used to read and preprocess geological exploration data and heavy metal sampling data, and to complete stratigraphic generalization and coordinate unification transformation; The physical boundary establishment module is used to construct an adaptive planar mesh, stretch it to generate a three-dimensional stratigraphic framework, and establish the physical spatial boundaries of heterogeneous strata. The solid model building module is used to generate a three-dimensional geological solid model based on a three-dimensional stratigraphic framework and assign stratigraphic attribute codes to each voxel unit; The controlled interpolation and mesh generation module is used to perform three-dimensional kriging interpolation calculations constrained by the stratigraphic boundaries and output a layered contaminated mesh aligned with the three-dimensional stratigraphic framework. The deep coupling analysis module is used to realize the attribute mapping between the layered pollution grid and the three-dimensional geological entity model, and to analyze the pollution distribution status in different strata; The dynamic risk assessment module is used to quantitatively calculate the volume of contaminated material, determine migration paths, and output a three-dimensional targeted remediation spatial control model and risk assessment report.

[0015] Preferably, the multi-source data acquisition module includes: The data parsing unit is used to read geological survey data files and laboratory test result spreadsheets in different formats, and extract the three-dimensional coordinates, layer thickness, lithological parameters, and concentration values ​​of various heavy metals from the borehole. The coordinate transformation unit is used to perform projection transformation operations on the extracted spatial coordinate data and unify it to a preset standard three-dimensional spatial coordinate system.

[0016] Preferably, the dynamic risk assessment module includes: The voxel integral calculation unit is used to perform volumetric summation calculations on the input strata and voxel units above a specific concentration threshold, and outputs statistical values ​​of contaminated soil volume. The spatial vector analysis unit is used to read the spatial distribution geometry of high pollutant concentration areas and the corresponding physical and mechanical parameters of the strata, and to determine and output vector data representing the migration path of lateral retention or infiltration along fracture channels. The control model generation unit is used to generate a three-dimensional targeted repair space control model based on the migration path vector data, which consists of a set of coordinates for shallow excavation barrier boundaries and a set of coordinates for deep in-situ injection repair nodes. The output unit is used to output a pollution risk distribution data array and assessment report for the study area based on quantitative pollution volume statistics, identified lateral migration mechanisms, and spatial location and depth information of deep intrusion phenomena along fissures.

[0017] Compared with the prior art, the beneficial effects of the present invention are: This invention eliminates the numerical smoothing effect across heterogeneous physical strata by using the physical spatial boundary of the formation as a constraint condition for spatial interpolation calculation of heavy metal concentration, thus making the spatial estimation of concentration conform to the actual barrier or permeability characteristics of different formation structures. Deeply superimposing and coupling the controlled-generated layered pollution grid with the three-dimensional geological entity model can objectively analyze the discrete and continuous distribution states of various heavy metals within geological layers with different physical properties. The quantitative statistics of contaminated soil volume based on spatial attributes and the assessment of differentiated migration patterns based on geological physical and mechanical parameters provide technical parameter basis based on three-dimensional spatial coordinates for the design depth of barrier engineering and the delineation of targeted remediation range for precise risk control measures of contaminated sites. The interpolation smoothing effect across heterogeneous strata was effectively blocked, allowing the pollution simulation results to objectively reflect the geological control effect. Attached Figure Description

[0018] 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 stratigraphic coupled method for soil pollution simulation and migration path assessment according to the present invention; Figure 2 This is a block diagram of a three-dimensional stratigraphic coupled soil pollution simulation and migration path assessment system according to the present invention. Figure 3 This is a rendering of a three-dimensional geological entity model constructed based on physical spatial boundaries in an embodiment of the present invention; Figure 4 This is a schematic diagram of the computational flow (node ​​flow) for controlled spatial interpolation and hierarchical contamination mesh generation in an embodiment of the present invention; Figure 5 This is a schematic diagram showing the absolute pollution distribution of four different heavy metal (Pb, Cd, As, Hg) exceeding the standard areas in each generalized stratum in the embodiments of the present invention. Figure 6 This is a visual schematic diagram of the three-dimensional stratigraphic framework and generalized geological stratification in an embodiment of the present invention; Figure 7 This is a schematic diagram simulating the initial distribution pattern of heavy metal pollution sources in shallow fill soil in an embodiment of the present invention; Figure 8 This is a spatial distribution characteristic diagram of heavy metal pollutants penetrating the clay layer and undergoing lateral enrichment in localized areas, as described in an embodiment of the present invention. Figure 9 This is a three-dimensional evolution diagram of the deep intrusion and migration of heavy metal pollutants along the fracture channels of the underlying weathered rock layer in an embodiment of the present invention; Figure 10 This is a detailed module diagram of the multi-source data acquisition module in the three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment system of this invention; Figure 11This is a specific module diagram of the dynamic risk assessment module in the three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment system of the present invention. Detailed Implementation

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

[0020] In the interdisciplinary field of soil environmental science and geological engineering, the spatial distribution information of heavy metal pollution in soil typically exhibits non-uniform characteristics. Existing soil pollution status assessment techniques usually rely on field survey data, combined with two-dimensional spatial interpolation algorithms to simulate the planar distribution of a single heavy metal element at different sampling depths, or employ conventional three-dimensional spatial interpolation algorithms to perform isotropic concentration estimation across the entire computational domain.

[0021] The vertical distribution and migration of heavy metals in soil and subsurface media are influenced by physical stratigraphic structures. For example, different geological bodies, such as fill layers, low-permeability clay layers, and weathered silty mudstone layers with fractures, exert different physical controls on the migration, diffusion, and distribution of pollutants. In existing technical solutions, the spatial interpolation calculation of heavy metal concentrations is usually relatively independent of the three-dimensional physical stratigraphic model. The distance or covariance calculation of the interpolation algorithm does not take the physical boundaries of the strata as constraints, which may lead to the generated pollution distribution model exhibiting computational smoothing across impermeable strata, making it difficult to accurately reflect the controlling role of the strata on the spatial distribution and migration patterns of heavy metals.

[0022] To address the aforementioned technical problems, this invention discloses a three-dimensional stratum coupled soil pollution simulation and migration path assessment method and system. The physical spatial boundary of the strata is used as a constraint condition for the spatial interpolation calculation of heavy metal concentration, thereby analyzing the three-dimensional spatial distribution and volume of pollutants within different geological layers and assessing the migration patterns of heavy metals under corresponding geological conditions.

[0023] like Figure 1 As shown, a three-dimensional stratigraphic coupled method for soil pollution simulation and migration path assessment includes: Step S1: Obtain geological survey data and heavy metal sampling data for the study area. After preprocessing, the strata in the study area are generalized into six strata distributed from the surface downwards: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. Among them, the geological survey data includes borehole spatial coordinates (including longitude, latitude, and elevation information), borehole columnar section data, borehole layer data, strata lithology classification information, groundwater level depth data, and physical, mechanical, and hydrogeological parameters of each rock and soil layer (including but not limited to permeability coefficient, porosity, dry density, and water content). The heavy metal sampling data includes measured concentration data of arsenic (As), lead (Pb), cadmium (Cd), and mercury (Hg) at different depths.

[0024] In this embodiment, the preprocessing process includes data cleaning, outlier removal, and coordinate unification; spatial coordinate data from different sources and formats are uniformly converted to a preset standard three-dimensional spatial coordinate system to establish a unified spatial calculation benchmark. After data preprocessing, based on borehole columnar section data and stratigraphic lithology classification information, the strata in the study area were generalized into six strata distributed from the surface downwards: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. The generalized strata described above represent the main media units with different permeability and physical barrier properties within the study area, among which, The fill layer is configured as the main infiltration and storage layer for pollutants; Clay and silty clay layers, due to their low permeability coefficient, are configured as barrier layers or slow-seepage layers for pollutants. Weathered silty mudstone layers are configured as potential deep migration channels for pollutants, depending on their degree of weathering and the development of fissures.

[0025] Step S2: Based on the borehole spatial coordinates, an adaptive planar mesh is constructed using the Delaunay triangulation algorithm or the Voronoi diagram mesh generation algorithm. The planar mesh is stretched along the Z-axis by combining digital elevation model data and borehole layer data. Topological anti-intersection and node merging logic is configured to handle the formation pinch-out phenomenon, generating three-dimensional mesh surfaces of the top and bottom plates of each stratum, establishing the physical spatial boundaries of each heterogeneous stratum, and forming a three-dimensional stratigraphic framework.

[0026] In this embodiment, during the adaptive planar mesh generation process, the algorithm dynamically adjusts the mesh resolution based on the spatial density distribution of borehole data points: In areas with dense borehole distribution or large gradients in formation elevation, the grid side length is automatically reduced to increase node density; in areas with sparse borehole distribution or gentle changes in formation elevation, the grid size is increased accordingly. The side length adjustment logic for adaptive planar meshes is implemented through a preset mesh size control function: Let the side length of the basic grid in the study area be... For any evaluation node in the grid generation algorithm, calculate its preset search radius. Total number of boreholes And the absolute value of the rate of change of the elevation gradient at the location of the node. Then the grid side length The specific mathematical relation configuration is as follows: ; ; In the formula, Represented by natural constant An exponential function with base 0. This is the borehole density weighting coefficient, and its value range is configured as follows: , The elevation gradient weighting coefficient has a value range configured as follows: , The gradient vector of elevation in a horizontal two-dimensional plane. For partial derivative mathematical operations, and These are the x and y components of the horizontal two-dimensional plane, respectively. This refers to the elevation value of the corresponding grid node; When the calculated grid side length Less than the preset minimum grid limit threshold Forced configuration This is to control the upper limit of memory consumption in the computing domain; Further, digital elevation model data of the study area is obtained, and the digital elevation model data is mapped to the nodes of the adaptive planar grid, and the corresponding initial Z-axis coordinates of the ground surface elevation are assigned to the two-dimensional grid nodes. Extract the borehole stratification data of each generalized stratum recorded in the geological exploration borehole data, namely the burial depth data of the top plate and bottom plate; Based on the adaptive planar grid of the ground elevation, a three-dimensional spatial stretching operation is performed along the vertical downward direction of the Z-axis; During the stretching algorithm, the nodes of the adaptive planar mesh are anchored when they reach the elevation of the top and bottom plate interface of each generalized stratum determined by the borehole layering data, thereby generating a three-dimensional mesh surface representing each stratum interface in three-dimensional space. To address the stratigraphic pinch-out phenomenon present in actual geological exploration data, the 3D stratigraphic framework incorporates topological anti-intersection and node merging logic during the stretching construction process, namely: Let the second be an adaptive planar mesh. Node coordinates At this location, the first in the generalized stratigraphic sequence The elevation scale of the top surface of the layer is The bottom elevation scale is ,in, These correspond to six generalized strata, ranging from miscellaneous fill to completely weathered argillaceous siltstone. is the globally unique positive integer index of the node. , This represents the total number of nodes in the two-dimensional planar mesh. , The stretching calculation unit calculates the X-axis and Y-axis coordinate components of the node in a preset standard three-dimensional coordinate system. Vertical thickness at the location The formula for its calculation is: ; The algorithm pre-configures a threshold for determining the minimum formation thickness. ; When logical judgment conditions When established, it is determined that the first occurrence occurred at that spatial coordinate. The strata pinch out; The algorithm forces the node merging instruction to be executed, making The upcoming The top and bottom grid surfaces of the layer are in The topological nodes at the coordinates are forced to coincide in the Z-axis direction; By using this underlying constraint logic, abnormal geometric units with zero or negative thickness are eliminated, thereby generating a three-dimensional stratigraphic framework that satisfies the geometric conditions of a closed manifold without self-intersection.

[0027] Step S3: Based on the interface mesh of the three-dimensional stratigraphic framework, the three-dimensional boundary shell of the study area is generated by the convex hull mesh algorithm. The undulating area of ​​the stratigraphic interface is adaptively densified. The solid filling between the curved surfaces is realized by the three-dimensional mesh solidification algorithm to construct a three-dimensional geological solid model.

[0028] In this embodiment, constructing a three-dimensional geological entity model includes: First, the convex hull mesh algorithm is used to determine the spatial envelope of the outermost data points in order to generate the three-dimensional boundary shell of the study area; Secondly, adaptive mesh refinement technology is applied to areas with large undulations or intersections at the stratigraphic interfaces. By increasing the number of polygons on the local mesh surfaces, the geometric shape of the geological interfaces is smoothly expressed, thereby generating the generalized three-dimensional surfaces of the interfaces of the six stratigraphic layers. Finally, using the surface surface generated from the digital elevation model data as the top interface and the three-dimensional surfaces of each internal stratum as the middle and bottom interfaces, the closure and solid filling between adjacent surfaces are achieved through the three-dimensional mesh solidification algorithm, thus constructing a three-dimensional geological solid model. Among them, the three-dimensional mesh solidification algorithm is either a voxelization algorithm or a solid filling algorithm based on the boundary representation method B-Rep; like Figure 3 As shown, this is a visualization of the three-dimensional geological entity model constructed in the embodiment. The model intuitively reflects the surface morphology and the three-dimensional spatial superposition relationship and undulation characteristics of heterogeneous strata such as fill soil, clay, silty clay and weathered silty mudstone layers of various levels. In the three-dimensional geological entity model, the space is divided into discrete voxel units or unstructured tetrahedral / hexahedral mesh units. Each three-dimensional geometric unit is configured and assigned a corresponding stratum attribute code.

[0029] Step S4: When performing spatial interpolation calculations on heavy metal sampling data to establish a three-dimensional distribution model of heavy metal pollution in soil, the layered grid of the aforementioned three-dimensional stratigraphic framework is forcibly used as the constraint condition for spatial interpolation. This ensures that the interpolation calculation of heavy metal concentration is controlled by the physical spatial boundary, thereby generating a layered pollution grid that is aligned with the three-dimensional stratigraphic framework in terms of geometric shape and spatial properties.

[0030] In this embodiment, the heavy metal sampling data includes the spatial coordinates of heavy metals such as arsenic, lead, cadmium, and mercury. and its measured concentration scalar value; A three-dimensional kriging interpolation algorithm is used to estimate concentration using three-dimensional mesh nodes as the computational carrier. The algorithm internally incorporates a spherical semivariogram function model to quantify spatial autocorrelation. The mathematical definition rules are as follows: when hour: ; when hour: ; In the formula, Let Euclidean distance be the three-dimensional straight-line distance between the grid node to be estimated and the known heavy metal sampling point. To quantify the semivariance function value of spatial autocorrelation, The nugget constant represents the spatial discontinuity caused by measurement errors and microscale variations. Its value is obtained by extrapolating the semivariance of the full sample data within the minimum sampling interval using a linear fit. Obtaining the intercept at that point, For the partial sill value, The sill value constitutes the total sample space variance of the entire interpolation calculation. The range threshold represents the maximum effective distance boundary of spatial correlation. It is determined as follows: A scatter plot of the experimental variogram of the measured samples is plotted; a spherical curve is fitted using the nonlinear least squares method; and the x-axis distance value corresponding to when the fitted curve reaches the sill level is extracted as the range threshold. ; In conventional kriging interpolation, the estimation of the concentration at unknown points relies on the weights of surrounding known sampling points calculated using a semi-variogram, the value of which typically depends only on the Euclidean distance between two points in space. In this embodiment, the interpolation algorithm is configured to use the layered mesh of the three-dimensional stratigraphic framework generated in step S2 as the boundary constraints for spatial interpolation. Understandably, when calculating the spatial correlation and distance weights between the grid node to be estimated and the known heavy metal sampling points, the distance calculation module or covariance calculation module inside the interpolation algorithm will execute ray casting and cross-detection logic; the algorithm will generate spatial line segments connecting the node to be estimated and the known sampling points, and detect whether the line segment geometrically intersects with any layered grid surface in the three-dimensional stratigraphic framework; if an intersection occurs, it indicates that the node to be estimated and the sampling point belong to different physical strata (for example, one is located in a miscellaneous fill layer and the other is located in a clay layer). When cross-stratum intersection is detected, the algorithm adjusts the input distance parameter in the semivariance function or directly modifies the weight matrix according to the preset constraint rules. One specific configuration method is: when crossing an interface with barrier properties (such as entering a clay layer), the effective calculation distance between the two points is configured to a preset penalty value, or the covariance contribution rate of the known sampling point to the node to be estimated is forced to be zero or a minimum value close to zero. The specific extraction steps and logical judgment rules for this constraint are defined as follows: In the weight matrix solution module, the three-dimensional spatial coordinates of the node to be estimated are set as follows. The three-dimensional spatial coordinates of the known sampling points currently involved in the calculation are: Calculate the initial three-dimensional straight-line distance between the two points as follows: .in, The nodes to be estimated are respectively X, Y, and Z axis coordinate components in a three-dimensional spatial coordinate system; These are known sampling points. The algorithm constructs a spatial line segment vector connecting two points, along with the X, Y, and Z axis coordinate components. And retrieve the set of interface mesh patches in the three-dimensional stratigraphic framework. ; Execute the ray-polygon intersection detection algorithm to determine the spatial line segment vector. Is it related to the set? Any grid patch with physical barrier properties (such as clay layers) intersects geometrically; Define Boolean variables for intersection determination If geometric intersection occurs, then configure Otherwise, configure The algorithm constructs a corrected spatial distance calculation. Its mathematical mapping relationship is configured as follows: ; In the formula, To calculate the distance in the corrected space, Let be the initial three-dimensional straight-line distance between the node to be estimated and the known sampling points. Boolean variables for intersection determination, The preset physical barrier penalty parameters, The value is forcibly configured to be greater than the aforementioned range threshold. At least two orders of magnitude more (for example, if Meters, then configuration (meters); triggered when crossing a physical barrier interface At that time, the corrected distance Much greater than the range threshold After substituting into the aforementioned spherical semivariance function, the known sampling point... Corresponding node to be estimated The spatial covariance value is forced to be calculated as zero; thus, through this judgment condition and the penalty mapping relationship, the numerical smoothing across heterogeneous physical strata is cut off. Through this conditional judgment and weight intervention mechanism based on physical interface cross-detection, the smoothing effect of Kriging interpolation is blocked by the physical boundary; the calculation gradient of the heavy metal concentration scalar field will produce discontinuities or jumps that conform to physical reality when crossing the boundary of heterogeneous strata, that is, the interpolation calculation of heavy metal concentration is controlled by the physical spatial boundary. After completing the constrained space interpolation calculation, the generated output is a hierarchical pollution grid, which contains a grid structure composed of three-dimensional nodes. Each node stores the concentration values ​​of various heavy metals obtained from the interpolation calculation. Its interpolation process is strictly controlled by the layered grid interface of the three-dimensional stratigraphic framework. The layered contamination grid is consistent and aligned with the three-dimensional stratigraphic framework in terms of geometric boundaries and spatial topology.

[0031] A schematic diagram of the computational flow (node ​​flow) for controlled spatial interpolation and hierarchical contamination mesh generation is shown below. Figure 4 As shown, this process is implemented through a visualized node flow. The specific node structure and connection relationships are described below: First, the system inputs multi-source data, which flows through the gridding and stratum processing nodes (reading the ground plane and generating a planar grid) to extract and generate a controlled basic planar grid. At the same time, the system imports the site boundary line through the CAD file import node (importing site boundary CAD) and uses the polygon triangulation node (converting CAD format) to complete the boundary data format conversion and triangulation. Subsequently, the read stratigraphic data passes through the stratigraphic-to-3D node (reads the mesh and generates 3D stratigraphic data) and flows into the decomposition and scaling series nodes (including decomposition and scaling #1 and decomposition and scaling #2), so that the 3D data is "exploded" and "scaled" according to the stratigraphic layers to finely adapt to the 3D geometric features of heterogeneous stratigraphic boundaries; on the pollution concentration estimation side, the heavy metal sampling data enters the 3D estimation #1 node (generates a pollution information model) to perform Kriging 3D estimation calculations; Subsequently, the stratigraphic spatial grid and the estimated pollution model are input together into a series of tangent nodes for calculating the two-dimensional regional distance (cutting the pollution / stratigraphic model according to the CAD range), and spatial trimming is completed according to the predetermined study area boundary. The trimmed model enters the core boundary constraint calculation stage, which includes intersection shell extraction nodes and multiple parallel intersection operation (screening out the pollutant range according to pollution type and threshold) node groups (such as intersection operation #2, intersection operation #3, intersection operation #4). The above nodes use the stratigraphic surface as the geometric barrier judgment surface, calculate the physical intersection conditions, and screen out the range of pollutant voxels above a specific threshold. Finally, the selected heavy metal controlled pollution grids are connected to the corresponding result rendering nodes (i.e. lead element viewer, cadmium element viewer, arsenic element viewer, mercury element viewer), and supplemented with legends, 3D annotations, direction indicators, coordinate axes and other marker nodes. All tributary data are finally aggregated to the ultimate master viewer node (3D visualization window / map rendering) to complete the output of the composite model and 3D rendering.

[0032] Step S5 involves matching the spatial coordinates of the layered pollution grid with the three-dimensional geological entity model. The heavy metal concentration attribute and the stratum attribute are jointly mapped using an octree index and an axial bounding box collision detection algorithm to form a three-dimensional composite entity model. Based on the stratum attribute and concentration threshold filtering conditions, the absolute distribution of heavy metal pollution within different strata is analyzed in detail.

[0033] In this embodiment, step S5 includes: The process involves traversing all nodes or voxels in the layered pollution grid, extracting the heavy metal concentration scalar information they carry, and mapping this concentration information to the corresponding coordinate units in the three-dimensional geological entity model according to the spatial coordinate matching principle. Each smallest computational unit in the three-dimensional composite entity model after mapping contains three-dimensional spatial coordinate data, stratigraphic lithology data, and concentration attribute data of multiple heavy metals. Based on the three-dimensional composite entity model after mapping, the system is configured with spatial attribute query and Boolean operation logic. By setting the stratum attribute filtering conditions and concentration threshold filtering conditions, the spatial distribution of pollutants inside a specific stratum can be analyzed and extracted respectively. Among them, the efficiency of mapping and matching operations in the three-dimensional coordinate system is improved by using octree indexing and axis bounding box collision detection algorithms. The specific mapping extraction steps are as follows: First, using the global spatial cube boundary of the 3D geological entity model as the root node, recursively construct a maximum depth of... The octree index assigns all voxel units carrying formation attributes to the corresponding octree leaf nodes. Secondly, for any target node in the layered pollution grid, let its three-dimensional spatial coordinates be... The system generates this node. With geometric center and side length as The bounding box of the cube ,in, The edge length is the globally unique positive integer index of the target node in the contaminated mesh. Forced configuration to the minimum side length of the reference voxel in the 3D geological solid model times; Finally, by traversing the octree index, the bounding box can be quickly retrieved and queried. The set of leaf nodes that exhibit spatial geometric overlap is identified, and a 3D point containment test algorithm is invoked to determine the unique containing node in 3D space. The target voxel unit.

[0034] Based on the aforementioned attribute-coupled 3D model data structure, the system is configured with spatial attribute query and Boolean operation logic. By setting stratum attribute filtering conditions and concentration threshold filtering conditions, the system can separately analyze and extract the spatial distribution of pollutants within specific strata. For example, by inputting a query command, it can extract all voxel sets where the stratum attribute is fill soil and the Pb concentration is greater than a preset limit, reconstructing the occurrence morphology of heavy metals in the fill soil layer in a 3D visualization form; or it can extract voxel sets where the stratum attribute is clay to analyze the concentration gradient decay state of pollutants within the barrier layer; or it can extract voxel sets where the stratum attribute is moderately weathered silty mudstone to analyze the spatial distribution characteristics of heavy metals in deep rock masses; such as... Figure 5 As shown, this is a detailed analysis of various heavy metals ( Figure 5 Pb in (A) Figure 5 Cd in (B) Figure 5 As in (C) Figure 5 This is a three-dimensional spatial decomposition diagram showing the absolute pollution distribution of Hg (D) within various generalized strata (from top to bottom: fill, clay, silty clay, and various weathered argillaceous siltstone layers). The darker areas visually represent the three-dimensional spatial extension of high-contamination zones where the corresponding heavy metal concentration exceeds its preset risk limit. Through this superposition and coupling step, the discrete and continuous distribution states of various heavy metals in geological strata with different physical properties are clearly analyzed. This step analyzes the discrete and continuous distribution states of various heavy metals in geological strata with different physical properties.

[0035] Step S6: Based on the absolute pollution distribution status, quantitatively calculate the pollution volume of each heavy metal in different strata, and combine the strata physical properties to assess the differentiated migration patterns of heavy metals along specific channels and the potential risk of deep pollution.

[0036] Step S6 includes: Quantitative statistics of contamination volume: First, for each heavy metal element and each generalized stratum, the system receives preset environmental risk control concentration limits as input parameters; Secondly, the statistical algorithm traverses all voxel units within the corresponding stratum, compares the heavy metal concentration attribute values ​​of the voxel units with the set concentration limits, and filters out all target voxel units whose concentration attribute values ​​are greater than or equal to the set limits. Finally, by summing and integrating the geometric volumes of all selected target voxel units, the total volume of contaminated soil for each heavy metal in different generalized strata, including fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone, is quantitatively calculated and output. Let the nominal geometric volume of a single reference voxel be... The scalar concentration of a certain heavy metal contained within this voxel is [value missing]. The preset environmental risk control concentration limit for this type of heavy metal is Define a dimensionless volume correction factor. ; When the concentration inside the target voxel Furthermore, the concentration values ​​of the voxel in the three-dimensional space of all 26 adjacent neighboring voxels are greater than or equal to the concentration of the voxel. When configuring ; When the internal concentration of voxels However, there exists at least one neighboring voxel with a concentration value less than [a certain value]. When the voxel is located on the three-dimensional boundary of the contamination isosurface, the isosurface function is calculated by calling a three-dimensional linear interpolation algorithm. The effective oversized polyhedron volume formed by truncating the target voxel, and the effective oversized polyhedron volume and the nominal geometric volume are then compared. The ratio is assigned to ; in, For arbitrary spatial coordinates within a 3D computational domain generated based on spatial interpolation The continuous heavy metal concentration distribution function at a given location; Total absolute volume of a single type of heavy metal in contaminated soil within a specific generalized stratum The mathematical calculation relation is configured as follows: ; In the formula, This represents the total absolute volume of contaminated soil containing a single type of heavy metal within a specific generalized stratum. The nominal geometric volume of a single reference voxel unit. For the logical conditions to be met within this specific stratum The total number of voxel units, The iterative index of the effective voxels participating in the accumulation of volume integrals; Quantitative assessment of migration patterns: Extract the continuous set of voxels with excessive concentrations within a specific formation, and calculate the concentration gradient vector using the three-dimensional central difference method. Combined with the permeability coefficient scalar of the corresponding strata Effective porosity parameters Concentration limits for heavy metal risk control Construct a normalized dimensionless concentration gradient Define the apparent migration potential vector controlled by formation permeability. The mathematical relationship is: ; ; In the formula, For heavy metal concentration variables in continuous space. , , These represent the horizontal, vertical, and axial coordinate components of the spatial coordinates within the three-dimensional computational domain. These are mathematical operators for partial derivatives; Principal component analysis (PCA) was used to perform dimensionality reduction on the extracted continuous voxel set's 3D coordinate point cloud data. The eigenvector direction corresponding to the first principal component was then extracted and used as the spatial dominance dip vector for high-pollution areas. ; Read the joint and fracture surface attitude information recorded in the geological exploration data, and generate a unit reference vector representing the dip and dip angle of the joint surface. ; Calculate the cosine similarity between two vectors The calculation formula is as follows: ; In the formula, The spatial dominance dip vector for high-pollution areas. As a unit reference vector representing the dip and dip angle of the joint surface, For the dot product of vectors, and Representing vectors respectively and The modulus length; System configuration similarity threshold When logical judgment conditions When it was established, the distribution direction of the pollutants was determined to be highly coincident with the direction of natural joints and fissures, and the assessment and identification results of heavy metal pollutants deep intrusion along specific joint and fissure channels were output. The specific analytical logic includes: In a layer of miscellaneous fill with high surface permeability, the longitudinal infiltration trend of pollutants driven by hydrodynamics was analyzed. At the interface between miscellaneous fill and low-permeability clay layer, the lateral distribution characteristics of the concentration gradient were extracted to assess the lateral migration and local enrichment mechanisms of pollutants after they were blocked. In fully weathered, strongly weathered, and moderately weathered argillaceous siltstone layers, the three-dimensional geometric morphology of high-pollution areas is analyzed. If the high-pollution areas in the weathered argillaceous siltstone layers are distributed in a linear, vein-like, or finger-like pattern extending along a specific spatial dip angle, and this spatial distribution pattern matches the spatial occurrence of the joint and fissure development zone recorded in the geological exploration data, then the phenomenon of heavy metal pollutants locally intruding into the deep geological environment along the dominant flow migration channel of the weathered layer fissures is identified and confirmed. like Figure 6 The figure shows a visualization of the three-dimensional stratigraphic framework and generalized geological strata generated in the embodiment. The figure clearly and intuitively presents the six main generalized stratigraphic units distributed from the surface downwards in the study area: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. The framework accurately restores the three-dimensional spatial superposition relationship of strata with different permeability and barrier properties through the three-dimensional undulation and solid filling of their respective interfaces, and constitutes the physical spatial base boundary for subsequent heavy metal spatial interpolation. When combining a three-dimensional coupling model to perform a refined evaluation of the above-mentioned transfer mechanism, a variety of typical transfer features can be visualized and mapped: like Figure 7As shown, this is a schematic diagram simulating the initial distribution of heavy metal pollution sources in shallow fill soil in the embodiment. The dark area in the figure shows the vertical infiltration trend and initial spatial occurrence of heavy metal pollutants in the fill soil layer with high surface porosity and permeability under the hydrodynamic drive of early natural rainfall and other factors. like Figure 8 The diagram illustrates the spatial distribution characteristics of heavy metal pollutants penetrating the clay layer and undergoing localized lateral enrichment in this embodiment. It is clearly observed that when the pollutant plume migrates to the interface between the fill and the underlying low-permeability clay layer, the vertical infiltration of the pollutants is restricted due to the physical obstruction and low-permeability slowing effect of the clay layer. The concentration gradient exhibits a significant lateral diffusion characteristic at this point, forming a large-scale typical lateral enrichment zone along the undulating interface of the strata. like Figure 9 The diagram illustrates the three-dimensional evolution of heavy metal pollutants' deep intrusion and migration along the fissures of the underlying weathered rock strata in an embodiment of the present invention. Combined with the dark continuum element morphology in the diagram, it can be seen that the high-pollution area not only penetrates the upper silty clay layer but also exhibits a linear, vein-like, or finger-like distribution extending along a specific spatially dominant dip angle in the underlying strongly weathered silty mudstone, moderately weathered silty mudstone, and completely weathered silty mudstone strata. This spatial morphology reconstruction objectively and accurately identifies the concealed phenomenon of heavy metal pollutants locally intruding into the deep geological environment along natural dominant channels such as structural fissures in the weathered rock mass. Generation of 3D targeted repair model: Based on the identification of lateral migration and local enrichment mechanisms, as well as the phenomenon of heavy metal pollutants locally intruding into the deep geological environment along the dominant flow migration channels of weathering layer fissures, the three-dimensional spatial coordinates corresponding to the relevant phenomena are further extracted to generate a three-dimensional targeted remediation spatial control model. This model includes a set of coordinates for shallow excavation barrier boundaries and a set of coordinates for deep in-situ injection remediation nodes. The generation steps include: The first step is to generate a set of shallow excavation barrier boundary coordinates based on lateral migration and local enrichment mechanisms, namely: For the lateral enrichment region at the interface between the fill and the low-permeability clay layer identified in the assessment, traverse the region to satisfy the logical conditions. For all target voxel units, extract the horizontal and vertical coordinate point cloud of the target voxel units in the horizontal two-dimensional plane, and apply the convex hull algorithm to generate the outermost two-dimensional closed polygon surrounding the point cloud. And extend a preset safety buffer distance outward along the normal direction of the polygon boundary. The two-dimensional planar boundary for the repair excavation is generated, where the value range is forcibly configured based on the empirical porosity of the site fill soil. ; Simultaneously, the system iterates through the Z-axis coordinate components in the three-dimensional spatial coordinates of the aforementioned target voxel units, extracting the minimum value as the elevation scalar of the excavation barrier bottom interface. The x and y coordinates of the extended two-dimensional plane boundary and Together, we establish and output the set of three-dimensional coordinates of the shallow excavation barrier boundary; The second step involves generating a set of coordinates for deep in-situ injection remediation nodes based on the phenomenon of heavy metal pollutants locally intruding into the deep geological environment through dominant flow migration channels along weathering layer fissures. For satisfying logical judgment conditions Furthermore, the target voxel set within the weathered rock layer where deep intrusion occurred has been identified, assuming that this set contains... Individual unit; The system calculates and extracts the coordinates of the three-dimensional geometric center of the voxel set. ;in, , , This is the arithmetic mean of the X, Y, and Z coordinate components of all voxel units within the set in a preset standard three-dimensional spatial coordinate system; Simultaneously, the system extracts the maximum value of the Z-axis coordinate component from all voxel elements within the set as the elevation of the highest point. Extract the minimum value of the Z-axis coordinate components as the elevation of the lowest point. ; Subsequently, with Starting from the spatial origin, along the spatial dominance dip vector extracted by the aforementioned dimensionality reduction calculation. The determined spatial straight path, defined in and Within the vertical depth range between them, according to the preset injection node spacing Perform equidistant spatial node sampling; where, The spatial linear Euclidean distance between two adjacent injection points is defined as the physical parameter of connectivity of fractures in weathered rock mass, and its value range is forcibly configured as follows: ; The series of discrete three-dimensional coordinate points generated by the above equidistant sampling constitute the set of coordinates of deep in-situ injection repair nodes, which are used to limit the injection depth and spatial target placement range of the agent in the actual site repair process. Output: Based on quantitative pollution volume statistics, identified lateral migration mechanisms, and spatial location and depth information of deep intrusion phenomena along fissures, a pollution risk distribution data array and assessment report for the study area are output.

[0037] like Figure 2As shown, a three-dimensional stratigraphic coupled soil pollution simulation and migration path assessment system includes: a multi-source data acquisition module 1, a physical boundary establishment module 2, an entity model construction module 3, a controlled interpolation and mesh generation module 4, a deep coupling analysis module 5, and a dynamic risk assessment module 6.

[0038] Multi-source data acquisition module 1 is used to read and preprocess geological exploration data and heavy metal sampling data, and to complete stratigraphic generalization and coordinate unification transformation.

[0039] like Figure 10 As shown, the multi-source data acquisition module 1 includes: a data parsing unit 11 and a coordinate transformation unit 12; The data parsing unit 11 is used to read geological survey data files and laboratory test result spreadsheets in different formats, and extract the three-dimensional coordinates, layer thickness, lithological parameters and concentration values ​​of various heavy metals of the borehole. The coordinate transformation unit 12 is used to perform projection transformation operations on the extracted spatial coordinate data and unify it to the preset standard three-dimensional spatial coordinates.

[0040] Physical boundary establishment module 2 is used to construct an adaptive planar mesh, stretch it to generate a three-dimensional stratigraphic framework, and establish the physical spatial boundary of heterogeneous strata.

[0041] Module 3, the entity model building module, is used to generate a three-dimensional geological entity model based on a three-dimensional stratigraphic framework and assign stratigraphic attribute codes to each voxel unit.

[0042] Controlled interpolation and mesh generation module 4 is used to perform three-dimensional kriging interpolation calculations constrained by the stratigraphic boundary and output a layered contaminated mesh aligned with the three-dimensional stratigraphic framework.

[0043] The deep coupling analysis module 5 is used to realize the attribute mapping between the layered pollution grid and the three-dimensional geological entity model, and to analyze the pollution distribution status in different strata.

[0044] The dynamic risk assessment module 6 is used to quantitatively calculate the volume of contaminated material, determine migration paths, and output a three-dimensional targeted remediation spatial control model and risk assessment report.

[0045] like Figure 11 As shown, the dynamic risk assessment module 6 includes a voxel integral calculation unit 61 and a spatial vector analysis unit 62; Voxel integration calculation unit 61 is used to perform volume accumulation calculation on input strata and voxel units above a specific concentration threshold, and output statistical values ​​of contaminated soil volume. The spatial vector analysis unit 62 is used to read the spatial distribution geometry of high pollutant concentration areas and the physical and mechanical parameters of the corresponding strata, and to determine and output the migration path vector data representing lateral retention or infiltration along fracture channels. The control model generation unit 63 is used to generate a three-dimensional targeted repair space control model based on the migration path vector data, which consists of a set of coordinates of shallow excavation barrier boundaries and a set of coordinates of deep in-situ injection repair nodes. Output unit 64 is used to output a pollution risk distribution data array and assessment report for the study area based on quantitative pollution volume statistics, identified lateral migration mechanisms, and spatial location and depth information of deep intrusion phenomena along fissures.

[0046] 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 stratigraphic coupled method for soil pollution simulation and migration path assessment, characterized in that, include: Step S1: Obtain geological survey data and heavy metal sampling data for the study area. After preprocessing, the strata in the study area are generalized into six strata distributed from the surface downwards: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. Among them, the geological survey data includes borehole spatial coordinates, borehole columnar section data, borehole stratification data, strata lithology classification information, groundwater level depth data, and physical, mechanical, and hydrogeological parameters of each rock and soil layer. The heavy metal sampling data includes measured concentration data of arsenic, lead, cadmium, and mercury at different depths. Step S2: Based on the borehole spatial coordinates, an adaptive planar mesh is constructed using the Delaunay triangulation algorithm or the Voronoi diagram mesh generation algorithm. The planar mesh is stretched along the Z-axis by combining digital elevation model data and borehole layer data. Topological anti-intersection and node merging logic is configured to handle the formation pinch-out phenomenon, generate three-dimensional mesh surfaces of the top and bottom plates of each stratum, establish the physical spatial boundaries of each heterogeneous stratum, and form a three-dimensional stratigraphic framework. Step S3: Based on the interface mesh of the three-dimensional stratigraphic framework, the three-dimensional boundary shell of the study area is generated by the convex hull mesh algorithm, the undulating area of ​​the stratigraphic interface is adaptively densified, and the solid filling between the curved surfaces is realized by the three-dimensional mesh solidification algorithm to construct a three-dimensional geological solid model. Step S4: When performing spatial interpolation calculations on heavy metal sampling data to establish a three-dimensional distribution model of soil heavy metal pollution, the layered grid of the above-mentioned three-dimensional stratigraphic framework is forcibly used as the constraint condition for spatial interpolation, so that the interpolation calculation of heavy metal concentration is controlled by the physical spatial boundary, thereby generating a layered pollution grid that is aligned with the three-dimensional stratigraphic framework in terms of geometric shape and spatial properties. Step S5: Spatial coordinate matching is performed between the layered pollution grid and the three-dimensional geological entity model. The heavy metal concentration attribute and the stratum attribute are jointly mapped by the octree index and the axial bounding box collision detection algorithm to form a three-dimensional composite entity model. Based on the stratum attribute and concentration threshold filtering conditions, the absolute distribution of heavy metal pollution in different strata is refined and analyzed. Step S6: Based on the absolute pollution distribution status, quantitatively calculate the pollution volume of each heavy metal in different strata, and combine the strata physical properties to assess the differentiated migration patterns of heavy metals along specific channels and the potential deep pollution risk. Quantitative statistics of contamination volume: For each heavy metal element and each generalized stratum, the system receives preset environmental risk control concentration limits as input parameters; The statistical algorithm traverses all voxel units within the corresponding stratum, compares the heavy metal concentration attribute values ​​of the voxel units with the set concentration limits, and filters out all target voxel units whose concentration attribute values ​​are greater than or equal to the set limits. By summing and integrating the geometric volumes of all selected target voxel units, the total volume of contaminated soil for each heavy metal in different generalized strata such as miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone is quantitatively calculated and output. Let the nominal geometric volume of a single reference voxel be... The scalar concentration of a certain heavy metal contained within this voxel is [value missing]. The preset environmental risk control concentration limit for this type of heavy metal is Define a dimensionless volume correction factor. ; When the concentration inside the target voxel Furthermore, the concentration values ​​of the voxel in the three-dimensional space of all 26 adjacent neighboring voxels are greater than or equal to the concentration of the voxel. When configuring ; When the internal concentration of voxels However, there exists at least one neighboring voxel with a concentration value less than [a certain value]. When the voxel is located on the three-dimensional boundary of the contamination isosurface, the isosurface function is calculated by calling a three-dimensional linear interpolation algorithm. The effective oversized polyhedron volume formed by truncating the target voxel, and the effective oversized polyhedron volume and the nominal geometric volume are then compared. The ratio is assigned to ; in, For arbitrary spatial coordinates within a 3D computational domain generated based on spatial interpolation The continuous heavy metal concentration distribution function at a given location; Total absolute volume of a single type of heavy metal in contaminated soil within a specific generalized stratum The mathematical calculation relation is configured as follows: ； In the formula, This represents the total absolute volume of contaminated soil containing a single type of heavy metal within a specific generalized stratum. The nominal geometric volume of a single reference voxel unit. For the logical conditions to be met within this specific stratum The total number of voxel units, The iterative index of the effective voxels participating in the accumulation of volume integrals; Quantitative assessment of migration patterns: Extract the continuous set of voxels with excessive concentrations within a specific formation, and calculate the concentration gradient vector using the three-dimensional central difference method. Combined with the permeability coefficient scalar of the corresponding strata Effective porosity parameters Concentration limits for heavy metal risk control Construct a normalized dimensionless concentration gradient Define the apparent migration potential vector controlled by formation permeability. The mathematical relationship is: ； ； In the formula, For heavy metal concentration variables in continuous space. , , These represent the horizontal, vertical, and axial coordinate components of the spatial coordinates within the three-dimensional computational domain. These are mathematical operators for partial derivatives; Principal component analysis (PCA) was used to perform dimensionality reduction on the extracted continuous voxel set's 3D coordinate point cloud data. The eigenvector direction corresponding to the first principal component was then extracted and used as the spatial dominance dip vector for high-pollution areas. ; Read the joint and fracture surface attitude information recorded in the geological exploration data, and generate a unit reference vector representing the dip and dip angle of the joint surface. ; Calculate the cosine similarity between two vectors The calculation formula is as follows: ； In the formula, The spatial dominance dip vector for high-pollution areas. As a unit reference vector representing the dip and dip angle of the joint surface, For the dot product of vectors, and Representing vectors respectively and The modulus length; System configuration similarity threshold When logical judgment conditions When established, the direction of pollution distribution and extension is determined to be highly coincident with the direction of natural joints and fissures, and the assessment and identification results of heavy metal pollutants deeply invading along specific joint and fissure channels are output.

2. The method for three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment according to claim 1, characterized in that, In step S1, the preprocessing process includes data cleaning, outlier removal, and coordinate unification; spatial coordinate data from different sources and formats are uniformly converted to a preset standard three-dimensional spatial coordinate system to establish a unified spatial calculation benchmark. After data preprocessing, based on borehole columnar section data and stratigraphic lithology classification information, the strata in the study area were generalized into six strata distributed from the surface downwards: miscellaneous fill, clay, silty clay, strongly weathered argillaceous siltstone, moderately weathered argillaceous siltstone, and completely weathered argillaceous siltstone. The generalized strata described above represent the main media units with different permeability and physical barrier properties within the study area, among which, The fill layer is configured as the main infiltration and storage layer for pollutants; Clay and silty clay layers, due to their low permeability coefficient, are configured as barrier layers or slow-seepage layers for pollutants. Weathered silty mudstone layers are configured as potential deep migration channels for pollutants, depending on their degree of weathering and the development of fissures.

3. The method for three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment according to claim 1, characterized in that, In step S2, during the adaptive planar mesh generation process, the algorithm dynamically adjusts the mesh resolution based on the spatial density distribution of borehole data points. In areas with dense borehole distribution or large gradients in formation elevation, the grid side length is automatically reduced to increase node density; in areas with sparse borehole distribution or gentle changes in formation elevation, the grid size is increased accordingly. The side length adjustment logic for adaptive planar meshes is implemented through a preset mesh size control function: Let the side length of the basic grid in the study area be... For any evaluation node in the grid generation algorithm, calculate its preset search radius. Total number of boreholes And the absolute value of the rate of change of the elevation gradient at the location of the node. Then the grid side length The specific mathematical relation configuration is as follows: ； ； In the formula, Represented by natural constant An exponential function with base 0. This is the borehole density weighting coefficient, and its value range is configured as follows: , The elevation gradient weighting coefficient has a value range configured as follows: , The gradient vector of elevation in a horizontal two-dimensional plane. For partial derivative mathematical operations, and These are the x and y components of the horizontal two-dimensional plane, respectively. This refers to the elevation value of the corresponding grid node; When the calculated grid side length Less than the preset minimum grid limit threshold Forced configuration This is to control the upper limit of memory consumption in the computing domain; Obtain digital elevation model data of the study area, map the digital elevation model data onto the nodes of the adaptive planar grid, and assign the corresponding initial Z-axis coordinates of the ground surface elevation to the two-dimensional grid nodes; Extract the borehole stratification data of each generalized stratum recorded in the geological exploration borehole data, namely the burial depth data of the top plate and bottom plate; Based on the adaptive planar grid of the ground elevation, a three-dimensional spatial stretching operation is performed along the vertical downward direction of the Z-axis; During the stretching algorithm, the nodes of the adaptive planar mesh are anchored when they reach the elevation of the top and bottom plate interface of each generalized stratum determined by the borehole layering data, thereby generating a three-dimensional mesh surface representing each stratum interface in three-dimensional space. To address the stratigraphic pinch-out phenomenon present in actual geological exploration data, the 3D stratigraphic framework incorporates topological anti-intersection and node merging logic during the stretching construction process, namely: Let the second be an adaptive planar mesh. Node coordinates At this location, the first in the generalized stratigraphic sequence The elevation scale of the top surface of the layer is The bottom elevation scale is ,in, These correspond to six generalized strata, ranging from miscellaneous fill to completely weathered argillaceous siltstone. is the globally unique positive integer index of the node. , This represents the total number of nodes in the two-dimensional planar mesh. , The stretching calculation unit calculates the X-axis and Y-axis coordinate components of the node in a preset standard three-dimensional coordinate system. Vertical thickness at the location The formula for its calculation is: ； The algorithm pre-configures a threshold for determining the minimum formation thickness. ; When logical judgment conditions When established, it is determined that the first occurrence occurred at that spatial coordinate. The strata pinch out; The algorithm forces the node merging instruction to be executed, making The upcoming The top and bottom grid surfaces of the layer are in The topological nodes at the coordinates are forced to coincide in the Z-axis direction; By using this underlying constraint logic, abnormal geometric units with zero or negative thickness are eliminated, thereby generating a three-dimensional stratigraphic framework that satisfies the geometric conditions of a closed manifold without self-intersection.

4. The method for three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment according to claim 1, characterized in that, In step S3, constructing the three-dimensional geological entity model includes: The convex hull mesh algorithm is used to determine the spatial envelope of the outermost data points in order to generate the three-dimensional boundary shell of the study area; Adaptive mesh refinement technology is applied to areas with large undulations or intersections at the stratigraphic interfaces. By increasing the number of polygons on the local mesh surfaces, the geometric shape of the geological interfaces is smoothly expressed, thereby generating the generalized three-dimensional surfaces of the interfaces of the six stratigraphic layers. Using the surface surface generated from digital elevation model data as the top interface, and the three-dimensional surfaces of each internal stratum as the middle and bottom interfaces, a three-dimensional geological solid model is constructed by using a three-dimensional mesh solidification algorithm to achieve closure and solid filling between adjacent surfaces. Among them, the three-dimensional mesh solidification algorithm is either a voxelization algorithm or a solid filling algorithm based on the boundary representation method B-Rep; In the three-dimensional geological entity model, the space is divided into discrete voxel units or unstructured tetrahedral / hexahedral mesh units, and each three-dimensional geometric unit is configured and assigned a corresponding stratigraphic attribute code.

5. The method for three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment according to claim 1, characterized in that, In step S4, the heavy metal sampling data includes the spatial coordinates of heavy metals such as arsenic, lead, cadmium, and mercury. and its measured concentration scalar value; A three-dimensional kriging interpolation algorithm is used to estimate concentration using three-dimensional mesh nodes as the computational carrier. The algorithm internally incorporates a spherical semivariogram function model to quantify spatial autocorrelation. The mathematical definition rules are as follows: when hour: ； when hour: ； In the formula, Let Euclidean distance be the three-dimensional straight-line distance between the grid node to be estimated and the known heavy metal sampling point. To quantify the semivariance function value of spatial autocorrelation, The nugget constant represents the spatial discontinuity caused by measurement errors and microscale variations. Its value is obtained by extrapolating the semivariance of the full sample data within the minimum sampling interval using a linear fit. Obtaining the intercept at that point, For the partial sill value, The sill value constitutes the total sample space variance of the entire interpolation calculation. The range threshold represents the maximum effective distance boundary of spatial correlation. It is determined as follows: A scatter plot of the experimental variogram of the measured samples is plotted; a spherical curve is fitted using the nonlinear least squares method; and the x-axis distance value corresponding to when the fitted curve reaches the sill level is extracted as the range threshold. ; After completing the constrained space interpolation calculation, the generated output is a hierarchical pollution grid, which contains a grid structure composed of three-dimensional nodes. Each node stores the concentration values ​​of various heavy metals obtained from the interpolation calculation. Its interpolation process is strictly controlled by the layered grid interface of the three-dimensional stratigraphic framework. The layered contamination grid is consistent and aligned with the three-dimensional stratigraphic framework in terms of geometric boundaries and spatial topology.

6. The method for three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment according to claim 1, characterized in that, Step S5 includes: The process involves traversing all nodes or voxels in the layered pollution grid, extracting the heavy metal concentration scalar information they carry, and mapping this concentration information to the corresponding coordinate units in the three-dimensional geological entity model according to the spatial coordinate matching principle. Each smallest computational unit in the three-dimensional composite entity model after mapping contains three-dimensional spatial coordinate data, stratigraphic lithology data, and concentration attribute data of multiple heavy metals. Based on the three-dimensional composite entity model after mapping, the system is configured with spatial attribute query and Boolean operation logic. By setting the stratum attribute filtering conditions and concentration threshold filtering conditions, the spatial distribution of pollutants inside a specific stratum can be analyzed and extracted respectively. Among them, the efficiency of mapping and matching operations in the three-dimensional coordinate system is improved by using octree indexing and axis bounding box collision detection algorithms. The specific mapping extraction steps are as follows: First, using the global spatial cube boundary of the 3D geological entity model as the root node, recursively construct a maximum depth of... The octree index assigns all voxel units carrying formation attributes to the corresponding octree leaf nodes. Secondly, for any target node in the layered pollution grid, let its three-dimensional spatial coordinates be... The system generates this node. With geometric center and side length as The bounding box of the cube ,in, The edge length is the globally unique positive integer index of the target node in the contaminated mesh. Forced configuration to the minimum side length of the reference voxel in the 3D geological solid model times; Finally, by traversing the octree index, the bounding box can be quickly retrieved and queried. The set of leaf nodes that exhibit spatial geometric overlap is identified, and a 3D point containment test algorithm is invoked to determine the unique containing node in 3D space. The target voxel unit.

7. The method for three-dimensional stratigraphic coupling soil pollution simulation and migration path assessment according to claim 1, characterized in that, Step S6 include: The specific analytical logic includes: In a layer of miscellaneous fill with high surface permeability, the longitudinal infiltration trend of pollutants driven by hydrodynamics was analyzed. At the interface between miscellaneous fill and low-permeability clay layer, the lateral distribution characteristics of the concentration gradient were extracted to assess the lateral migration and local enrichment mechanisms of pollutants after they were blocked. In fully weathered, strongly weathered, and moderately weathered argillaceous siltstone layers, the three-dimensional geometric morphology of high-pollution areas is analyzed. If the high-pollution areas in the weathered argillaceous siltstone layers are distributed in a linear, vein-like, or finger-like pattern extending along a specific spatial dip angle, and this spatial distribution pattern matches the spatial occurrence of the joint and fissure development zone recorded in the geological exploration data, then the phenomenon of heavy metal pollutants locally intruding into the deep geological environment along the dominant flow migration channel of the weathered layer fissures is identified and confirmed. Generation of 3D targeted repair model: Based on the identification of lateral migration and local enrichment mechanisms, as well as the phenomenon of heavy metal pollutants locally intruding into the deep geological environment along the dominant flow migration channels of weathering layer fissures, the three-dimensional spatial coordinates corresponding to the relevant phenomena are further extracted to generate a three-dimensional targeted remediation spatial control model. This model includes a set of coordinates for shallow excavation barrier boundaries and a set of coordinates for deep in-situ injection remediation nodes. The generation steps include: The first step is to generate a set of shallow excavation barrier boundary coordinates based on lateral migration and local enrichment mechanisms, namely: For the lateral enrichment region at the interface between the fill and the low-permeability clay layer identified in the assessment, traverse the region to satisfy the logical conditions. For all target voxel units, extract the horizontal and vertical coordinate point cloud of the target voxel units in the horizontal two-dimensional plane, and apply the convex hull algorithm to generate the outermost two-dimensional closed polygon surrounding the point cloud. And extend a preset safety buffer distance outward along the normal direction of the polygon boundary. The two-dimensional planar boundary for the repair excavation is generated, where the value range is forcibly configured based on the empirical porosity of the site fill soil. ; Simultaneously, the system iterates through the Z-axis coordinate components in the three-dimensional spatial coordinates of the aforementioned target voxel units, extracting the minimum value as the elevation scalar of the excavation barrier bottom interface. The x and y coordinates of the extended two-dimensional plane boundary and Together, we establish and output the set of three-dimensional coordinates of the shallow excavation barrier boundary; The second step involves generating a set of coordinates for deep in-situ injection remediation nodes based on the phenomenon of heavy metal pollutants locally intruding into the deep geological environment through dominant flow migration channels along weathering layer fissures. For satisfying logical judgment conditions Furthermore, the target voxel set within the weathered rock layer where deep intrusion occurred has been identified, assuming that this set contains... Individual unit; The system calculates and extracts the coordinates of the three-dimensional geometric center of the voxel set. ;in, , , This is the arithmetic mean of the X, Y, and Z coordinate components of all voxel units within the set in a preset standard three-dimensional spatial coordinate system; Simultaneously, the system extracts the maximum value of the Z-axis coordinate component from all voxel elements within the set as the elevation of the highest point. Extract the minimum value of the Z-axis coordinate components as the elevation of the lowest point. ; Subsequently, with Starting from the spatial origin, along the spatial dominance dip vector extracted by the aforementioned dimensionality reduction calculation. The determined spatial straight path, defined in and Within the vertical depth range between them, according to the preset injection node spacing Perform equidistant spatial node sampling; where, The spatial linear Euclidean distance between two adjacent injection points is defined as the physical parameter of connectivity of fractures in weathered rock mass, and its value range is forcibly configured as follows: ; The series of discrete three-dimensional coordinate points generated by the above equidistant sampling constitute the set of coordinates of deep in-situ injection repair nodes, which are used to limit the injection depth and spatial target placement range of the agent in the actual site repair process. Output: Based on quantitative pollution volume statistics, identified lateral migration mechanisms, and spatial location and depth information of deep intrusion phenomena along fissures, a pollution risk distribution data array and assessment report for the study area are output.

8. A three-dimensional stratigraphically coupled soil pollution simulation and migration path assessment system according to any one of claims 1-7, characterized in that, include: The multi-source data acquisition module is used to read and preprocess geological exploration data and heavy metal sampling data, and to complete stratigraphic generalization and coordinate unification transformation; The physical boundary establishment module is used to construct an adaptive planar mesh, stretch it to generate a three-dimensional stratigraphic framework, and establish the physical spatial boundaries of heterogeneous strata. The solid model building module is used to generate a three-dimensional geological solid model based on a three-dimensional stratigraphic framework and assign stratigraphic attribute codes to each voxel unit; The controlled interpolation and mesh generation module is used to perform three-dimensional kriging interpolation calculations constrained by the stratigraphic boundaries and output a layered contaminated mesh aligned with the three-dimensional stratigraphic framework. The deep coupling analysis module is used to realize the attribute mapping between the layered pollution grid and the three-dimensional geological entity model, and to analyze the pollution distribution status in different strata; The dynamic risk assessment module is used to quantitatively calculate the volume of contaminated material, determine migration paths, and output a three-dimensional targeted remediation spatial control model and risk assessment report.

9. The three-dimensional stratigraphic coupled soil pollution simulation and migration path assessment system according to claim 8, characterized in that, The multi-source data acquisition module includes: The data parsing unit is used to read geological survey data files and laboratory test result spreadsheets in different formats, and extract the three-dimensional coordinates, layer thickness, lithological parameters, and concentration values ​​of various heavy metals from the borehole. The coordinate transformation unit is used to perform projection transformation operations on the extracted spatial coordinate data and unify it to a preset standard three-dimensional spatial coordinate system.

10. The three-dimensional stratigraphic coupled soil pollution simulation and migration path assessment system according to claim 8, characterized in that, The dynamic risk assessment module includes: The voxel integral calculation unit is used to perform volumetric summation calculations on the input strata and voxel units above a specific concentration threshold, and outputs statistical values ​​of contaminated soil volume. The spatial vector analysis unit is used to read the spatial distribution geometry of high pollutant concentration areas and the corresponding physical and mechanical parameters of the strata, and to determine and output vector data representing the migration path of lateral retention or infiltration along fracture channels. The control model generation unit is used to generate a three-dimensional targeted repair space control model based on the migration path vector data, which consists of a set of coordinates for shallow excavation barrier boundaries and a set of coordinates for deep in-situ injection repair nodes. The output unit is used to output a pollution risk distribution data array and assessment report for the study area based on quantitative pollution volume statistics, identified lateral migration mechanisms, and spatial location and depth information of deep intrusion phenomena along fissures.

Images