Soil restoration-oriented soil data visual management method and system

Abstract

The invention relates to the technical field of environmental engineering and geoscience statistics, in particular to a soil remediation-oriented soil data visual management method and system.The method comprises the steps that S1, plot survey data are obtained and discretized to generate a three-dimensional voxel grid, constructing a three-dimensional geologic model based on the drilling lithology data, mapping a geologic lithology label to a voxel grid, and establishing a three-dimensional space data structure with stratum constraint; s2, analyzing underground water flow field data to determine a main flow direction azimuth angle of pollutant migration, and accordingly constructing an anisotropic search ellipsoid and setting variation function parameters; according to the method, the hybrid model combining machine learning trend prediction and anisotropic residual error correction is constructed, the global trend is fitted by utilizing the stratum lithology data, and the local deviation is corrected by combining the hydrodynamic factor, so that the smoothing effect of the traditional single interpolation method under the heterogeneous geological condition is improved; and the goodness of fit between a three-dimensional pollutant distribution simulation result and an actual geological environment is improved.

Description

Technical Field

[0001] This invention relates to the fields of environmental engineering and geoscientific statistics, and in particular to a method and system for visual management of soil data for soil remediation. Background Technology

[0002] In pollution investigation and remediation projects, accurately grasping the spatial distribution of soil pollutants and the volume of contaminated soil is the core basis for formulating remediation plans and controlling remediation costs. It usually requires the construction of a three-dimensional spatial distribution model based on limited borehole sampling and detection data, using geostatistical methods or numerical simulation technology, in order to estimate and visualize the concentration of pollutants in unsampled underground areas.

[0003] However, existing technologies often employ purely geometric or single statistical algorithms such as inverse distance weighting or ordinary kriging interpolation when constructing 3D models. These methods often assume that the spatial medium is homogeneous or isotropic, neglecting the control effect of the heterogeneity of strata lithology on the impediment or infiltration of pollutants and the guiding effect of groundwater flow on the migration of pollutants. This results in the simulated pollution plume morphology being too smooth and lacking physical mechanism support, making it difficult to accurately reflect the real pollution boundary under complex geological conditions, and consequently causing a large deviation in the estimation of remediation earthwork volume. Summary of the Invention

[0004] To overcome the above shortcomings, this invention provides a soil data visualization management method and system for soil remediation, aiming to improve the problem that existing technologies often use a single interpolation algorithm based on the isotropic assumption, which ignores the actual control effect of stratigraphic lithological heterogeneity and groundwater flow dynamics on pollutant migration.

[0005] In a first aspect, the present invention provides the following technical solution: a soil data visualization management method for soil remediation, comprising: S1. Obtain land survey data and discretize it to generate a three-dimensional voxel grid. Based on borehole lithology data, construct a three-dimensional geological model and map geological lithology labels to the voxel grid to establish a three-dimensional spatial data structure with stratigraphic constraints. S2. Analyze the groundwater flow field data to determine the azimuth angle of the mainstream direction of pollutant migration, and construct an anisotropic search ellipsoid based on this and set the variogram parameters to establish an anisotropic spatial weighting system. S3. Based on the three-dimensional spatial data structure, extract multi-dimensional feature vectors containing geological and coordinate information, train the model using machine learning regression algorithm, and generate a global pollutant concentration trend field. S4. Calculate the deviation between the measured value of the sampling point and the corresponding value of the trend field to obtain the residual, and use the anisotropic search ellipsoid to perform spatial interpolation on the residual data to generate the residual correction field. S5. The global pollutant concentration trend field and the residual correction field are numerically superimposed to generate a three-dimensional pollution distribution model, and the pollution plume view and remediation earthwork volume statistics are output based on the model.

[0006] Preferably, the establishment of the three-dimensional spatial data structure with stratigraphic constraints specifically includes the following steps: Set the grid resolution and initialize the generation of a zero-value voxel matrix with three-dimensional coordinate indices; The depths of the top and bottom plates of each layer are extracted by analyzing borehole data, and a continuous three-dimensional stratum interface surface model is constructed using interpolation algorithms. Perform spatial topological determination of the centroid coordinates of the voxel and the three-dimensional stratigraphic interface surface model to determine the stratigraphic level to which each voxel belongs; Based on the stratigraphic position, the corresponding lithological and hydrological parameters are matched and written into the voxel attribute table as geological lithology labels.

[0007] Preferably, the construction of the anisotropic search ellipsoid specifically includes the following steps: The hydraulic gradient vector is calculated based on multi-point water level monitoring data to determine the azimuth angle of the main flow direction of pollutant migration; The main azimuth angle is defined as the main axis direction of the search ellipsoid, and the horizontal vector perpendicular to this direction and the vector perpendicular to the ground direction are defined as the horizontal secondary axis direction and the vertical secondary axis direction, respectively. The ratio of the spatial search radius between the primary axis and the secondary axis is set according to the permeability coefficient parameter of the target stratum. By combining the ratio of the principal axis direction to the spatial search radius, the geometric definition parameters of the anisotropic search ellipsoid are generated.

[0008] Preferably, the establishment of the anisotropic spatial weighting system specifically includes the following steps: The experimental semivariance scatter plots of pollutant concentrations are calculated along the principal axis and the vertical direction of the anisotropic search ellipsoid, respectively. The experimental semivariance scatter points were fitted using a theoretical variogram model to obtain the nugget value, sill value, and range parameters in each direction. An anisotropic covariance function is constructed using the nugget value, sill value, and range parameter to generate a spatial weight coefficient matrix.

[0009] Preferably, the generation of the global pollutant concentration trend field specifically includes the following steps: A subset of features with known sampling point locations is selected from the three-dimensional spatial data structure and correlated with the measured concentration to construct a supervised learning dataset; An ensemble learning algorithm based on decision trees was selected and iteratively trained using the supervised learning dataset to obtain a pollutant concentration prediction model. Input the multidimensional feature vectors of all voxels to be predicted in the three-dimensional spatial data structure into the pollutant concentration prediction model to calculate the predicted concentration value; The predicted concentration values ​​are mapped back to a three-dimensional voxel grid using voxel indices to generate a global pollutant concentration trend field.

[0010] Preferably, the calculation of the residual specifically includes the following steps: Traverse the coordinates of the sampling points and map them to the three-dimensional voxel grid index to extract the trend prediction values ​​of the corresponding positions from the global pollutant concentration trend field; Perform numerical difference calculations to subtract the trend prediction value from the measured pollutant concentration to obtain the original residual data; The original residual data is subjected to a normal distribution test to construct a discrete residual point set containing spatial coordinates.

[0011] Preferably, generating the residual correction field specifically includes the following steps: The search range is defined using the anisotropic search ellipsoid with the voxel to be interpolated as the center; Discrete residual points falling within the search range are selected, and interpolation weight coefficients are obtained by solving the Kriging equations based on the anisotropic space weight system. The selected original residual data are weighted and summed using the interpolation weight coefficients to calculate the residual estimate of the current voxel; The estimation calculation is completed by traversing all voxels in the field, and the residual correction field is generated.

[0012] Preferably, the generation of the three-dimensional pollution distribution model specifically includes the following steps: Based on the grid index, the values ​​of the corresponding voxels in the global pollutant concentration trend field and the residual correction field are read, and a linear superposition operation is performed to obtain the composite concentration matrix; The synthetic concentration matrix is ​​subjected to non-negative constraint processing, and the values ​​less than zero in the calculation results are corrected to a preset minimum threshold. The concentration data after non-negative constraint processing is fused with the spatial geometric information of the three-dimensional voxel mesh, and the geological lithology labels are integrated to generate a three-dimensional pollution distribution model.

[0013] Preferably, the output of the pollution plume view and remediation earthwork volume statistics specifically includes the following steps: The three-dimensional pollution distribution model is screened based on the remediation target value to determine the set of voxels of pollutants exceeding the standard; The boundary model of the excessive pollutant voxel set is constructed by applying the isosurface extraction algorithm, and a three-dimensional pollutant plume visualization view is rendered and output. The volume of individual elements is calculated based on grid resolution, and the total earthwork volume for remediation is obtained by summing the accumulated volumes of the excessive pollutant elements. Based on the geological and lithological labels carried by the voxels, the total volume of earthwork to be repaired is classified and statistically analyzed into sub-items of different geological and lithological categories.

[0014] Secondly, the present invention provides the following technical solution: a soil data visualization management system for soil remediation, the system comprising: The data structure construction module is used to discretize and generate a three-dimensional voxel mesh and map geological lithology labels to establish a three-dimensional spatial data structure with stratigraphic constraints. The hydrodynamic rule definition module is used to analyze the groundwater flow field, construct an anisotropic search ellipsoid, and establish an anisotropic spatial weighting system. The trend field prediction module is used to extract multi-dimensional features based on the three-dimensional spatial data structure and generate a global pollutant concentration trend field through machine learning regression algorithms. The residual correction module is used to calculate the predicted residual and to perform spatial interpolation on the residual using the anisotropic search ellipsoid constructed by the hydrodynamic rule definition module to generate the residual correction field. The model reconstruction and output module is used to overlay the global pollutant concentration trend field and the residual correction field to generate a three-dimensional pollution distribution model, and output a three-dimensional pollution plume visualization view and statistics on the amount of earthwork to be repaired.

[0015] The present invention has the following beneficial effects: 1. In this invention, a hybrid model combining machine learning trend prediction and anisotropic residual correction is constructed. The global trend is fitted using stratigraphic lithology data and local deviations are corrected by combining hydrodynamic factors. This improves the smoothing effect of traditional single interpolation methods under heterogeneous geological conditions and enhances the consistency between the three-dimensional pollutant distribution simulation results and the actual geological environment.

[0016] 2. In this invention, an anisotropic spatial weighting system is constructed using groundwater flow field data. The physical direction of pollutant migration with groundwater is incorporated into the interpolation calculation, so that the morphology of the pollution plume can reflect the hydrogeological characteristics extending along the water flow, avoiding the morphological deviation caused by the traditional geometric interpolation ignoring the physical migration law.

[0017] 3. This invention establishes a three-dimensional model based on voxel grids carrying stratigraphic attributes, realizing the spatial alignment of pollutant concentration fields with geological lithology classification. This facilitates the accurate differentiation and statistical analysis of contaminated soil volumes in different geological media, providing quantitative data references for the classification and treatment of remediation projects and the calculation of engineering quantities. Attached Figure Description

[0018] Figure 1 This is a flowchart of a soil data visualization management method for soil remediation proposed in this invention; Figure 2This is an architecture diagram of a soil data visualization management system for soil remediation proposed in this invention. Detailed Implementation

[0019] The technical solutions in 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] Example 1: In a first embodiment of the present invention, the present invention provides a soil data visualization management method for soil remediation, such as... Figure 1 As shown, the method includes the following steps: S1. Obtain land survey data and discretize it to generate a three-dimensional voxel grid. Based on borehole lithology data, construct a three-dimensional geological model and map geological lithology labels to the voxel grid to establish a three-dimensional spatial data structure with stratigraphic constraints. Furthermore, establishing a three-dimensional spatial data structure with stratigraphic constraints specifically includes the following steps: Set the grid resolution and initialize the generation of a zero-value voxel matrix with three-dimensional coordinate indices; The depths of the top and bottom plates of each layer are extracted by analyzing borehole data, and a continuous three-dimensional stratum interface surface model is constructed using interpolation algorithms. Perform spatial topological determination of the centroid coordinates of the voxel and the three-dimensional stratigraphic interface surface model to determine the stratigraphic position of each voxel; Based on the stratigraphic position, the corresponding lithological and hydrological parameters are matched and written into the voxel attribute table as geological lithology labels.

[0021] Specifically, firstly, a three-dimensional spatial discretization process is performed. The boundary coordinate data of the survey area are read, and the minimum and maximum coordinate values ​​of the study area in the east-west, north-south, and vertical directions are determined. According to the preset calculation accuracy requirements, the grid step size in the east-west, north-south, and vertical directions is set respectively. Based on the above boundary coordinate values ​​and grid step size, a three-dimensional grid index system is constructed, and the three-dimensional spatial coordinates of the geometric centroid of each voxel unit in the grid system are calculated. The calculation formula for the centroid coordinates of the voxel unit is as follows: ; In the formula, , , These represent the coordinates of the voxel unit's centroid in the east-west, north-south, and vertical directions, respectively. , , These represent the initial minimum coordinate values ​​of the study area in the three directions, respectively. , , These represent the grid step sizes set in the three directions, respectively. , , These represent the positive integer index numbers of the voxel units in the three dimensions. At the same time, a three-dimensional array with the same dimensions as the grid index system is initialized in memory. The initial values ​​of all array elements are set to zero or null, which will be used to store geological properties later.

[0022] Next, spatial interpolation reconstruction of stratigraphic interfaces is performed, and borehole exploration data within the block is analyzed. For each borehole point, the boundary locations of different stratigraphic horizons are identified. For each stratigraphic level, the top and bottom burial depth data are extracted at all borehole locations. Combined with the absolute elevation of the borehole openings, the burial depth data are converted into absolute elevation data, forming a discrete stratigraphic interface elevation sample set. Using an inverse distance weighted interpolation algorithm, the continuous distribution surface of the stratigraphic interface across the entire field is calculated based on this sample set. The interpolation formula for the top stratigraphic interface is as follows: ; In the formula, Represents plane coordinates The absolute elevation value of the top plate interface of the stratum to be calculated. Indicates the first The absolute elevation of the top plate of the stratum measured at each borehole point. Represents the coordinates of the point to be calculated With the Euclidean distance between the coordinates of each borehole point This represents the total number of valid borehole sample points included in the calculation. This represents the inverse weighted exponent, typically a positive real number. Using the same algorithm and formula, the absolute elevation of the floor interface of this stratigraphic level is calculated. Through iterative calculations, the interface model of the top and bottom slabs of all geological layers within the site is constructed.

[0023] Finally, the geological properties are mapped and assigned, and each voxel in the 3D mesh indexing system is traversed to obtain its centroid coordinates. The planar coordinates of the centroid Substitute the above stratigraphic interface model into the query results for the top elevation at each corresponding location. Elevation of the base plate The execution space contains decision logic; if the centroid of the voxel is perpendicular to the coordinates... If a voxel is determined to be within a stratigraphic level if its base elevation is greater than or equal to that level and less than its top elevation, then the voxel is considered to be within that stratigraphic level. Based on the determined stratigraphic level name, the corresponding lithology code (e.g., CL for clay, SA for sand) and hydraulic permeability coefficient are retrieved from a pre-set geological parameter database. These attribute data are then written into a three-dimensional array with the index [index missing]. The storage unit is used to mark voxels that do not fall within any stratigraphic range as invalid or air cells. After traversal, the three-dimensional array is transformed into a three-dimensional spatial data structure with stratigraphic constraints, which serves as the geometric carrier for subsequent pollutant concentration prediction.

[0024] S2. Analyze the groundwater flow field data to determine the azimuth angle of the mainstream direction of pollutant migration, and construct an anisotropic search ellipsoid based on this and set the variogram parameters to establish an anisotropic spatial weighting system. Furthermore, constructing the anisotropic search ellipsoid specifically includes the following steps: The hydraulic gradient vector is calculated based on multi-point water level monitoring data to determine the azimuth angle of the main flow direction of pollutant migration; The main azimuth angle is defined as the main axis direction of the search ellipsoid, and the horizontal vector perpendicular to this direction and the vector perpendicular to the ground direction are defined as the horizontal secondary axis direction and the vertical secondary axis direction, respectively. The ratio of the spatial search radius between the primary axis and the secondary axis is set according to the permeability coefficient parameter of the target stratum. By combining the ratio of the principal axis direction to the spatial search radius, the geometric definition parameters of the anisotropic search ellipsoid are generated.

[0025] Furthermore, establishing an anisotropic spatial weighting system specifically includes the following steps: Calculate the experimental semivariance scatter points of pollutant concentrations along the principal axis and vertical direction of the anisotropic search ellipsoid, respectively. The experimental semivariance scatter plots were fitted using a theoretical variogram model to obtain the nugget value, sill value, and range parameters in each direction. An anisotropic covariance function is constructed using nugget value, sill value, and range parameter to generate a spatial weight coefficient matrix.

[0026] Specifically, the coordinate information of multiple groundwater monitoring wells within the site and the corresponding groundwater level elevation data at different times were first read to form a hydrological monitoring dataset. In order to determine the direction of the average hydraulic gradient across the entire site, the least squares method was used to fit the discrete water level data into a linear hydraulic gradient plane. A multiple linear regression equation was constructed to describe the functional relationship between groundwater level elevation and planar spatial coordinates. The equation is expressed as follows: ; Symbols in the formula Represents plane coordinates and Estimated groundwater level elevation at the location, Represents the regression intercept constant. This represents the rate of change of water level elevation in the east-west direction. The coefficient representing the rate of change of water level elevation in the north-south direction is used to obtain the hydraulic gradient vector reflecting the overall trend of the groundwater flow field by solving the coefficient of this equation. Based on this hydraulic gradient vector, the azimuth angle of the main flow direction of pollutant migration is calculated, which is the principal axis direction angle of the anisotropic search ellipsoid. The formula for calculating the azimuth angle of the main flow direction is as follows: ; Symbols in the formula This represents the azimuth angle of the mainstream direction of pollutant migration, with a value ranging from zero to 360 degrees. and These are the slope coefficients obtained from the above regression equation.

[0027] Next, the spatial geometric definition of the anisotropic search ellipsoid is constructed. The permeability coefficients of the target stratum along the principal axis and perpendicular to the principal axis are obtained, and their ratio is calculated as the anisotropy ratio. Using the principal axis azimuth and the anisotropy ratio, a coordinate rotation and scaling matrix is ​​constructed to transform the actual physical spatial coordinates into isotropic transformed spatial coordinates. Thus, the search ellipsoid is mathematically defined. The matrix operation formula for calculating the anisotropic distance between two points is as follows: ; Symbols in the formula This represents the anisotropic hydrological distance between two points. This represents the vector difference between the coordinates of two points in the original coordinate system. Indicates based on mainstream azimuth angle The constructed rotation matrix is ​​used to rotate the coordinate system to align with the direction of the water flow. This represents a diagonal scaling matrix constructed based on anisotropic ratios, used to stretch or compress coordinate differences along the secondary axis. The superscript indicates... This represents the matrix transpose operation.

[0028] Subsequently, the experimental semivariance of pollutant concentration was calculated to quantify spatial variability. A series of lag distance steps were set along the principal axis direction determined by the azimuth angle of the main current and along the secondary axis direction perpendicular to this direction. For each lag distance, all sampling point pairs with distance intervals matching that step size were searched across the entire field. The average of half the squares of the pollutant concentration differences between these point pairs was calculated to obtain the experimental semivariance. The formula for calculating the experimental semivariance is as follows: ; Symbols in the formula Indicates that at a lag distance of The experimental semivariance value at time, Indicates the distance lagged behind The number of paired sampling points that are separated. Indicates the location The measured pollutant concentration at the location, Indicates the location The measured pollutant concentration values ​​at the location.

[0029] Finally, a spatial weighting system is established, and a spherical model or exponential model is selected as the theoretical variogram model. The calculated experimental semivariogram scatter points are curve-fitted using the nonlinear least squares method to analyze the nugget value, sill value, and range parameters. Based on the fitted theoretical model, a covariance function is constructed, and a spatial weighting coefficient matrix is ​​generated accordingly. For any two spatial locations, their anisotropic distance is substituted into the covariance function to calculate the correlation. The general expression for the covariance function is as follows: ; Symbols in the formula This represents the covariance between two points. This represents the sill value, i.e., the total variance. Represents anisotropic distance The semivariance value obtained by substituting it into the fitted theoretical variogram model is used to form the coefficient matrix elements of the Kriging equation system, thus establishing the anisotropic spatial weighting system.

[0030] S3. Based on the three-dimensional spatial data structure, extract multi-dimensional feature vectors containing geological and coordinate information, use machine learning regression algorithms to train the model, and generate a global pollutant concentration trend field. Furthermore, generating a global pollutant concentration trend field specifically includes the following steps: A supervised learning dataset is constructed by selecting a feature subset with known sampling point locations from a three-dimensional spatial data structure and correlating it with measured concentrations. An ensemble learning algorithm based on decision trees was selected, and iterative training was performed using a supervised learning dataset to obtain a pollutant concentration prediction model. Input the multidimensional feature vectors of all voxels to be predicted in the three-dimensional spatial data structure into the pollutant concentration prediction model to calculate the predicted concentration value. The predicted concentration values ​​are mapped back to a three-dimensional voxel grid by voxel index to generate a global pollutant concentration trend field.

[0031] Specifically, feature engineering and training dataset construction are performed first. Records of all soil sampling points are read from the land parcel survey database. Each record contains a unique identifier for the sampling point, east-west coordinates, north-south coordinates, vertical sampling depth, and laboratory-measured pollutant concentration values. To associate continuous sampling point coordinates with discrete voxel grids, a floor function algorithm is used to calculate the grid index of each sampling point in the three-dimensional spatial data structure. The formula for calculating the grid index is as follows: ; In the formula , , These represent the integer indices of the sampling points in the voxel grid along the east-west, north-south, and vertical directions, respectively. , , These represent the measured three-dimensional spatial coordinates of the sampling points, , , These represent the initial minimum coordinate values ​​of the 3D voxel mesh in three directions, respectively. , , These represent the preset grid resolution step size, This indicates the floor function.

[0032] Next, based on the calculated grid index, the corresponding voxel unit is retrieved from the previously established 3D spatial data structure with stratigraphic constraints. The stratigraphic horizon code and geological lithology label stored in this voxel unit are then read. Since geological lithology labels are typically text-based classification data, numerical encoding of the classification features is required to enable them to be processed by machine learning regression algorithms. A geological lithology mapping dictionary is established, mapping various lithological descriptions such as clay, silty clay, silt, and sand to consecutive integers that increase from the beginning. Let the mapping function be... For the retrieved text tags Calculate its corresponding numerical characteristics equal For the first For each sampling point, its spatial coordinates, the retrieved stratigraphic horizon code, and the numerically encoded lithological features are combined to construct a one-dimensional feature row vector. This feature vector... The mathematical expression is as follows: ; In the formula , , For spatial coordinates, For stratigraphic integer encoding, For lithological numerical encoding, the feature row vectors of all valid sampling points are stacked along the row direction to form a dimension of The input feature matrix has five rows and five columns. ,in The total number of sampling points is [number], and the measured pollutant concentration value corresponding to each sampling point is extracted and arranged in the same order to form a dimension [dimension]. Target label vector in one row and one column Input feature matrix With target label vector Together, they constitute a supervised learning dataset for subsequent machine learning model training.

[0033] Subsequently, an ensemble learning algorithm based on decision trees was selected for model training, specifically the random forest regression algorithm. This algorithm establishes a nonlinear mapping relationship between input features and pollutant concentrations by integrating multiple independently trained decision tree regressors. Model parameters were initialized, setting the total number of decision trees in the ensemble and the maximum depth of a single tree. During model training, for each decision tree, a bootstrap resampling technique was used to extract a subset of samples from the training dataset. Nodes were then split in the feature space to grow the tree structure based on minimizing the mean squared error. The prediction output of the random forest model was obtained by taking the arithmetic mean of the prediction results from all decision trees. The mathematical formula is as follows:

[0034] In the formula Indicates that the input feature vector is Predicted pollutant concentration values ​​at that time; This represents the total number of decision trees included in the random forest model. Indicates the first Decision trees in parameters The regression prediction function under control, whose parameters include the split node variables, split threshold, and mean response of the leaf nodes in the tree structure. The input is a multidimensional feature vector containing coordinates and geological attributes. By inputting the training dataset into the algorithm for iterative calculation, the internal parameters of the model are optimized until the error of the model on the validation set converges, thereby obtaining the trained pollutant concentration prediction model.

[0035] Finally, batch prediction of global concentration and generation of trend field are performed. All voxel units in the three-dimensional spatial data structure are traversed, and the spatial coordinates, stratigraphic layer codes and geological lithology labels of each voxel are extracted one by one. The global feature matrix to be predicted is constructed in the same format as in the training stage. This global feature matrix to be predicted is input into the trained pollutant concentration prediction model, and forward inference calculation is performed to obtain the pollutant concentration prediction value of each voxel unit in the global field. Using the indexing system of the three-dimensional grid, the calculated prediction values ​​are mapped back to the corresponding voxel positions and stored as three-dimensional scalar field data. This three-dimensional scalar field is the global pollutant concentration trend field, which is continuously distributed in space and reflects the control effect of stratigraphic lithology on pollutant distribution.

[0036] S4. Calculate the deviation between the measured values ​​of the sampling points and the corresponding values ​​of the trend field to obtain the residuals, and use the anisotropic search ellipsoid to perform spatial interpolation on the residual data to generate the residual correction field. Furthermore, calculating the residuals specifically includes the following steps: Traverse the coordinates of the sampling points and map them to the three-dimensional voxel grid index to extract the trend prediction values ​​of the corresponding locations from the global pollutant concentration trend field; Perform numerical difference calculations to subtract the trend prediction value from the measured pollutant concentration to obtain the original residual data; Perform a normality test on the original residual data and construct a discrete residual point set containing spatial coordinates.

[0037] Furthermore, generating the residual correction field specifically includes the following steps: The search range is defined using an anisotropic search ellipsoid with the voxel to be interpolated as the center. Discrete residual points falling within the search range are selected, and interpolation weight coefficients are obtained by solving the Kriging equations based on the anisotropic space weight system. The original residual data selected by the interpolation weighting coefficients are weighted and summed to calculate the residual estimate of the current voxel; The estimation calculation is completed by traversing all voxels in the field, and the residual correction field is generated.

[0038] Specifically, the calculation of prediction residuals and dataset construction are performed first. All known sampling points within the plot are traversed, and the corresponding trend prediction values ​​are retrieved in the global pollutant concentration trend field using the three-dimensional coordinates of the sampling points. Numerical difference operations are then performed to calculate the difference between the measured pollutant concentration and the trend prediction value at each sampling point, thus obtaining the original residual data. The residual calculation formula is as follows: ; Symbols in the formula Indicates the location of the sampling point. The residual value at that point, This indicates the measured concentration of pollutants at that location. This represents the predicted value at this location provided by the global pollutant concentration trend field. The original residual data set obtained by calculation is tested for normality. If the data deviates from the normal distribution, a logarithmic transformation or Box-Cox transformation is performed on it to make it meet the statistical assumptions. The processed residual values ​​are combined with their corresponding three-dimensional spatial coordinates to construct a discrete residual point set, which serves as the source data for subsequent spatial interpolation.

[0039] Subsequently, residual estimates for the unsampled region are calculated using the anisotropic kriging interpolation algorithm. Each voxel cell to be interpolated in the 3D spatial data structure is traversed, with the centroid of that voxel cell as the search origin. The parameters of the anisotropic search ellipsoid constructed in the previous steps are called to define the effective neighborhood search range. The anisotropic distance between the centroid of the voxel to be interpolated and all discrete residual points is calculated. Valid sample points falling within the search ellipsoid range are selected. The calculation of the anisotropic distance must strictly adhere to the principal axis direction and anisotropic ratio parameters defined in step S2. Based on the selected valid sample points, a system of ordinary kriging linear equations is constructed to solve for the interpolation weights. The matrix expression of this system of linear equations is as follows: ; Symbols in the formula Indicates the first The effective sample point and the first The semivariance value between each valid sample point is calculated by substituting the anisotropic distance between the two points into the theoretical variogram model fitted in step S2. Indicates the first The semivariance between each effective sample point and the mass center of the volume to be interpolated. The term to be solved is the first... Interpolation weight coefficients for each valid sample point Denotes the Lagrange multiplier used to satisfy the unbiased estimation constraint. This represents the total number of valid sample points within the search range.

[0040] Finally, the system of equations is solved to generate the residual correction field. The Kriging equations are then solved using Gaussian elimination or matrix inversion to obtain the weighting coefficients for each valid sample point. These weighting coefficients are then used to perform a linear weighted summation of the residual values ​​within the search range, thus obtaining the residual estimate for the current voxel to be interpolated. The residual estimate calculation formula is as follows: ; Symbols in the formula Indicates the position of the voxel to be interpolated The residual estimate at the location, To obtain the weight coefficients, For the neighboring region Given the known residual values ​​of each sample point, the above search, weight calculation, and weighted calculation process is repeated for each voxel in the 3D mesh. The residual estimates of all voxels are filled into the 3D matrix according to the mesh index, and finally a residual correction field covering the entire domain is generated.

[0041] S5. The global pollutant concentration trend field and the residual correction field are numerically superimposed to generate a three-dimensional pollution distribution model, and the pollution plume view and remediation earthwork volume statistics are output based on the model.

[0042] Furthermore, generating a three-dimensional pollution distribution model specifically includes the following steps: Based on the grid index, the values ​​of corresponding voxels in the global pollutant concentration trend field and residual correction field are read, and a linear superposition operation is performed to obtain the composite concentration matrix; The synthesis concentration matrix is ​​subjected to non-negativity constraint processing, and the values ​​less than zero in the calculation results are corrected to the preset minimum threshold. The concentration data after non-negative constraint processing is fused with the spatial geometric information of the three-dimensional voxel mesh, and geological lithology labels are integrated to generate a three-dimensional pollution distribution model.

[0043] Furthermore, the output of the pollution plume view and remediation earthwork volume statistics specifically includes the following steps: The three-dimensional pollution distribution model was screened based on the remediation target value to determine the voxel set of pollutants exceeding the standard; The boundary model of the excessive pollutant voxel set is constructed by applying the isosurface extraction algorithm, and a three-dimensional visualization view of the pollutant plume is rendered and output. The volume of individual elements is calculated based on grid resolution, and the total earthwork volume for remediation is obtained by summing the accumulated volumes of the excessive pollutant elements. Based on the geological and lithological labels carried by the voxels, the total volume of earthwork to be repaired was classified and statistically analyzed into sub-items of different geological and lithological categories.

[0044] Specifically, firstly, numerical reconstruction and non-negativity constraint processing of the three-dimensional pollution distribution model are performed. Based on the unified indexing system of the three-dimensional voxel mesh, the global pollutant concentration trend field generated in step S3 and the residual correction field generated in step S4 are read simultaneously. For each voxel cell, its trend prediction value in the trend field and its residual estimate value in the residual correction field are extracted. A linear superposition operation is performed to obtain a preliminary estimated concentration. The linear superposition formula is as follows: ; In the formula Indicates that the grid index is Preliminary estimated concentration of pollutants at the location, This represents the corresponding trend forecast value. This represents the corresponding residual estimate. Given that Kriging interpolation may produce negative mathematical artifacts in sparse data regions, a non-negativity constraint correction needs to be applied to the initial estimated concentration. A preset minimum threshold is set, typically a value taken from the site's environmental background or the detection limit of the laboratory testing method. If the initial estimated concentration is less than this threshold, it is forcibly assigned to this threshold; if it is greater than or equal to this threshold, the original value is retained. The correction formula is as follows: ; In the formula This indicates the corrected final pollutant concentration. This indicates that a preset minimum threshold is used to write the corrected concentration data into the attribute table of the three-dimensional spatial data structure, and to deeply bind it with the spatial geometric coordinates of the voxels and the geological lithology labels, thereby generating the final three-dimensional pollution distribution model.

[0045] Subsequently, a visualization view of the pollution plume is generated based on the isosurface extraction algorithm. The soil remediation target value is set as the isosurface threshold. The moving cube algorithm is used to traverse the three-dimensional pollution distribution model. This algorithm treats each voxel unit as a cube, compares the pollutant concentration values ​​at the eight vertices of the cube with the remediation target value, determines the intersection points of the isosurface and the cube edges, and connects these intersection points into triangular patches by looking up a pre-set topological state table, thereby reconstructing a three-dimensional geometric mesh surface that fits the boundary of the pollution plume. This geometric mesh surface is imported into the three-dimensional rendering engine, given a semi-transparent material and color mapping, and outputs a plume view that shows the spatial morphology of the pollutants.

[0046] Finally, a layered and quality-based remediation earthwork volume statistics were performed. The physical volume of each voxel was calculated based on the grid resolution using the formula: east-west step size multiplied by north-south step size multiplied by vertical step size. This process was repeated across all voxels to identify the set of voxels with contaminant concentrations exceeding the remediation target value. The geological lithology tag carried by each voxel was read, and the voxel volumes were accumulated according to lithology category to calculate the contaminated earthwork volume for each lithology type. The formula for calculating the earthwork volume for lithological remediation is as follows: ; In the formula Indicates the first Total earthwork volume for lithological remediation Represents the volume of a single voxel. This represents the total number of voxels in the entire field. Let Kronecker function be the function of the first kinematic line. Lithological labels of individual elements With the The value is one when the lithology is consistent; otherwise, it is zero. Let be a step function, when the first... concentration of individual elements Greater than the target value for repair The value is 1 if the time condition is met, and 0 otherwise. Through this calculation, the final output is a statistical table containing the total earthwork volume and the earthwork volume of each sub-item based on lithology.

[0047] Example 2: In constructing 3D models, purely geometric or single statistical algorithms such as inverse distance weighting or ordinary kriging interpolation are often used. These methods tend to assume that the spatial medium is homogeneous or isotropic, neglecting the controlling effect of stratigraphic lithological heterogeneity on pollutant retention or infiltration, and the guiding effect of groundwater flow on pollutant migration. This results in simulated pollution plume morphologies that are too smooth and lack physical mechanism support. To address these issues, this invention provides a soil data visualization management system for soil remediation, the structure of which is as follows: Figure 2 As shown, the system includes: The data structure construction module is used to discretize and generate a three-dimensional voxel mesh and map geological lithology labels to establish a three-dimensional spatial data structure with stratigraphic constraints. The hydrodynamic rule definition module is used to analyze the groundwater flow field, construct an anisotropic search ellipsoid, and establish an anisotropic spatial weighting system. The trend field prediction module is used to extract multi-dimensional features based on the three-dimensional spatial data structure and generate a global pollutant concentration trend field through machine learning regression algorithms. The residual correction module is used to calculate the predicted residual and to perform spatial interpolation on the residual using the anisotropic search ellipsoid constructed by the hydrodynamic rule definition module to generate the residual correction field. The model reconstruction and output module is used to overlay the global pollutant concentration trend field and the residual correction field to generate a three-dimensional pollution distribution model, and output a three-dimensional pollution plume visualization view and statistics on the amount of earthwork to be repaired.

[0048] Specifically, the data structure construction module first reads the site exploration boundary and borehole columnar section data, uses a spatial discretization algorithm to divide the site into a set of three-dimensional voxels with unique indices, and writes the corresponding lithological attributes into the voxel attribute table based on the stratigraphic interpolation results to form the basic data architecture. The hydrodynamic rule definition module calculates the hydraulic gradient vector and the azimuth of the main current direction by parsing groundwater level monitoring data, and sets the major and minor axis ratios and rotation angle parameters of the anisotropic search ellipsoid and establishes a spatial weight matrix. The trend field prediction module calls the decision tree ensemble learning algorithm library, inputs the spatial coordinates of the voxels and lithological coding features for model training and inference, and outputs the global trend prediction value of pollutants across the entire field. The residual correction module performs the difference operation between the measured value and the trend value, and calls the ordinary Kriging interpolation program based on the anisotropic search parameters to generate a residual estimation field covering the entire region. The model reconstruction and output module merges the trend field and residual field through linear superposition operation and performs non-negative correction, uses the moving cube algorithm to render the three-dimensional pollution plume graphics, and simultaneously traverses the voxel grid to count the excess volume under different lithological labels to output a report of remediation earthwork volume.

[0049] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is 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 soil data visualization management method for soil remediation, characterized by, include: S1. Obtain land survey data and discretize it to generate a three-dimensional voxel grid. Based on borehole lithology data, construct a three-dimensional geological model and map geological lithology labels to the voxel grid to establish a three-dimensional spatial data structure with stratigraphic constraints. S2. Analyze the groundwater flow field data to determine the azimuth angle of the mainstream direction of pollutant migration, and construct an anisotropic search ellipsoid based on this and set the variogram parameters to establish an anisotropic spatial weighting system. S3. Based on the three-dimensional spatial data structure, extract multi-dimensional feature vectors containing geological and coordinate information, train the model using machine learning regression algorithm, and generate a global pollutant concentration trend field. S4. Calculate the deviation between the measured value of the sampling point and the corresponding value of the trend field to obtain the residual, and use the anisotropic search ellipsoid to perform spatial interpolation on the residual data to generate the residual correction field. S5. The global pollutant concentration trend field and the residual correction field are numerically superimposed to generate a three-dimensional pollution distribution model, and the pollution plume view and remediation earthwork volume statistics are output based on the model.

2. The soil data visualization management method for soil remediation according to claim 1, wherein, The establishment of a three-dimensional spatial data structure with stratigraphic constraints specifically includes the following steps: Set the grid resolution and initialize the generation of a zero-value voxel matrix with three-dimensional coordinate indices; The depths of the top and bottom plates of each layer are extracted by analyzing borehole data, and a continuous three-dimensional stratum interface surface model is constructed using interpolation algorithms. Perform spatial topological determination of the centroid coordinates of the voxel and the three-dimensional stratigraphic interface surface model to determine the stratigraphic level to which each voxel belongs; Based on the stratigraphic position, the corresponding lithological and hydrological parameters are matched and written into the voxel attribute table as geological lithology labels.

3. The soil data visualization management method for soil remediation according to claim 1, wherein, The construction of the anisotropic search ellipsoid specifically includes the following steps: The hydraulic gradient vector is calculated based on multi-point water level monitoring data to determine the azimuth angle of the main flow direction of pollutant migration; The main azimuth angle is defined as the main axis direction of the search ellipsoid, and the horizontal vector perpendicular to this direction and the vector perpendicular to the ground direction are defined as the horizontal secondary axis direction and the vertical secondary axis direction, respectively. The ratio of the spatial search radius between the primary axis and the secondary axis is set according to the permeability coefficient parameter of the target stratum. By combining the ratio of the principal axis direction to the spatial search radius, the geometric definition parameters of the anisotropic search ellipsoid are generated.

4. The soil data visualization management method for soil remediation according to claim 1, wherein, The establishment of the anisotropic spatial weighting system specifically includes the following steps: The experimental semivariance scatter plots of pollutant concentrations are calculated along the principal axis and the vertical direction of the anisotropic search ellipsoid, respectively. The experimental semivariance scatter points were fitted using a theoretical variogram model to obtain the nugget value, sill value, and range parameters in each direction. An anisotropic covariance function is constructed using the nugget value, sill value, and range parameter to generate a spatial weight coefficient matrix.

5. The soil data visualization management method for soil remediation according to claim 1, wherein, The generation of the global pollutant concentration trend field specifically includes the following steps: A subset of features with known sampling point locations is selected from the three-dimensional spatial data structure and correlated with the measured concentration to construct a supervised learning dataset; An ensemble learning algorithm based on decision trees was selected and iteratively trained using the supervised learning dataset to obtain a pollutant concentration prediction model. Input the multidimensional feature vectors of all voxels to be predicted in the three-dimensional spatial data structure into the pollutant concentration prediction model to calculate the predicted concentration value; The predicted concentration values ​​are mapped back to a three-dimensional voxel grid using voxel indices to generate a global pollutant concentration trend field.

6. The soil data visualization management method for soil remediation according to claim 1, wherein, The calculation of the residual specifically includes the following steps: Traverse the coordinates of the sampling points and map them to the three-dimensional voxel grid index to extract the trend prediction values ​​of the corresponding positions from the global pollutant concentration trend field; Perform numerical difference calculations to subtract the trend prediction value from the measured pollutant concentration to obtain the original residual data; The original residual data is subjected to a normal distribution test to construct a discrete residual point set containing spatial coordinates.

7. The soil data visualization management method for soil remediation according to claim 1, wherein, The generation of the residual correction field specifically includes the following steps: The search range is defined using the anisotropic search ellipsoid with the voxel to be interpolated as the center; Discrete residual points falling within the search range are selected, and interpolation weight coefficients are obtained by solving the Kriging equations based on the anisotropic space weight system. The selected original residual data are weighted and summed using the interpolation weight coefficients to calculate the residual estimate of the current voxel; The estimation calculation is completed by traversing all voxels in the field, and the residual correction field is generated.

8. The soil data visualization management method for soil remediation according to claim 1, wherein, The generation of the three-dimensional pollution distribution model specifically includes the following steps: Based on the grid index, the values ​​of the corresponding voxels in the global pollutant concentration trend field and the residual correction field are read, and a linear superposition operation is performed to obtain the composite concentration matrix; The synthetic concentration matrix is ​​subjected to non-negative constraint processing, and the values ​​less than zero in the calculation results are corrected to a preset minimum threshold. The concentration data after non-negative constraint processing is fused with the spatial geometric information of the three-dimensional voxel mesh, and the geological lithology labels are integrated to generate a three-dimensional pollution distribution model.

9. The soil data visualization management method for soil remediation according to claim 1, wherein, The output of the pollution plume view and the statistical data on remediation earthwork volume specifically includes the following steps: The three-dimensional pollution distribution model is screened based on the remediation target value to determine the set of voxels of pollutants exceeding the standard; The boundary model of the excessive pollutant voxel set is constructed by applying the isosurface extraction algorithm, and a three-dimensional pollutant plume visualization view is rendered and output. The volume of individual elements is calculated based on grid resolution, and the total earthwork volume for remediation is obtained by summing the accumulated volumes of the excessive pollutant elements. Based on the geological and lithological labels carried by the voxels, the total volume of earthwork to be repaired is classified and statistically analyzed into sub-items of different geological and lithological categories.

10. A soil data visualization management system for soil remediation, characterized by, A soil data visualization management method for soil remediation as described in any one of claims 1-9, the system comprising: The data structure construction module is used to discretize and generate a three-dimensional voxel mesh and map geological lithology labels to establish a three-dimensional spatial data structure with stratigraphic constraints. The hydrodynamic rule definition module is used to analyze the groundwater flow field, construct an anisotropic search ellipsoid, and establish an anisotropic spatial weighting system. The trend field prediction module is used to extract multi-dimensional features based on the three-dimensional spatial data structure and generate a global pollutant concentration trend field through machine learning regression algorithms. The residual correction module is used to calculate the predicted residual and to perform spatial interpolation on the residual using the anisotropic search ellipsoid constructed by the hydrodynamic rule definition module to generate the residual correction field. The model reconstruction and output module is used to overlay the global pollutant concentration trend field and the residual correction field to generate a three-dimensional pollution distribution model, and output a three-dimensional pollution plume visualization view and statistics on the amount of earthwork to be repaired.

Images