An Uncertainty-Based Method for 3D Reconstruction and Intelligent Sampling of Underground Space

Abstract

This invention discloses an uncertainty-based method for three-dimensional reconstruction and intelligent sampling of underground space, comprising the following steps: Step 1: Constructing a multi-source heterogeneous data acquisition system and performing data preprocessing and spatiotemporal registration; Step 2: Extracting multimodal heterogeneous features and unifying spatial mapping based on deep learning; Step 3: Constructing a high-fidelity three-dimensional reconstruction volumetric model of underground space based on adaptive octree subdivision; Step 4: Constructing a high-fidelity intelligent sampling database containing a Bayesian dynamic update layer; Step 5: Performing uncertainty assessment-driven sampling-assisted decision-making to update the database. This invention solves the problems of insufficient depth of multi-source data fusion, low fidelity due to rigid three-dimensional modeling mesh, and lack of an uncertainty-based active sampling feedback mechanism in the prior art, achieving accurate characterization of complex strata and pollution distribution and intelligent feedback of sampling decisions.

Description

Technical Field

[0001] This invention relates to the field of environmental analysis technology, and in particular to an uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space. Background Technology

[0002] Currently, with the support of hydraulic direct-push acquisition technology, the acquired samples generally have high fidelity and can be used for subsequent analysis. However, in complex site environments, due to the combined effects of variable groundwater and geological conditions and the heterogeneous distribution of pollutants, the high fidelity of soil samples obtained solely from discrete sampling points is still insufficient. Further technological support is needed to deeply integrate multi-source heterogeneous data, accurately reconstruct the in-situ physical properties of the site that couple macroscopic and microscopic dimensions, dynamically store the data, and implement intelligent sampling decision-making.

[0003] In response to the above industry needs, technical personnel in related fields have made some improvements, but the existing technical solutions still have significant limitations in practical applications.

[0004] For example, patent application CN116737817B proposes a "method, device, and computer-readable storage medium for fusing multi-source heterogeneous data." Although it proposes a general fusion framework that includes spatiotemporal, government, and IoT data, and achieves data integration through coordinate unification and format conversion, this solution is mainly aimed at urban surface facilities and does not fully consider the special characteristics of underground environments. Specifically, this technology lacks a fusion mechanism for "process data" specific to contaminated sites, and does not include sampling operating parameters (such as drilling pressure and rotation speed) and drilling rig behavior characteristics in the modeling scope, resulting in the model being unable to reflect the impact of sampling disturbances on the true value of the data; moreover, its fusion results are only used for visualization and lack an analysis module based on the degree of conflict of heterogeneous data, making it difficult to directly guide the optimization of subsequent sampling point layout.

[0005] Furthermore, patent application CN118378319B proposes a "three-dimensional unit body attribute modeling method and system based on multi-source heterogeneous information fusion." Although this scheme generates a three-dimensional model using drilling-while-sense data in tunnel engineering, it is difficult to directly apply to environmental geotechnical engineering. On the one hand, this scheme focuses on the mechanical stability of the surrounding rock and lacks a description of the coupling mechanism between the chemical field (pollution concentration) and the physical field (stratum structure), making it unable to accurately depict the diffusion morphology of the pollution plume. On the other hand, its data source is limited to local construction data and does not incorporate large-scale remote sensing image textures and macroscopic information such as environmental monitoring (e.g., groundwater level), making it difficult to achieve unified modeling of "macro background - micro structure." In addition, such models are usually statically constructed and lack a posterior update mechanism based on Bayesian theory. When new sampled data is acquired, the uncertainty field of the model cannot be corrected in real time, causing the database to fail to achieve version iteration and accuracy convergence as the project progresses.

[0006] Furthermore, patent application CN118817738B proposes a "three-dimensional reconstruction method for soil." While this patent utilizes CT scanning and porosity correction technology to achieve high-precision reconstruction of the pore structure of soil samples at the millimeter level, the technical solution in this application still has shortcomings in macroscopic site applications. This method is only suitable for microscopic soil sample analysis at the laboratory scale, and the reconstructed information depends entirely on expensive CT scans. When facing macroscopic contaminated sites of tens of thousands of square meters, it cannot be directly applied due to excessive computational load and scale effects. Moreover, this method cannot solve the problem of gridded sampling points in macroscopic sites, making it difficult to guide the optimization of on-site sampling points and the decision-making of operating parameters through uncertainty assessment.

[0007] In terms of site data management, while existing environmental information management systems have achieved digital data storage, their databases often remain static and inflexible when faced with high-fidelity sampling requirements. Existing systems typically serve only as static data repositories, resulting in the loss of crucial drilling process data and a lack of a dynamic update kernel based on Bayesian theory, making it impossible to use new data to retrospectively correct older models. In summary, existing technologies generally suffer from insufficient data fusion depth, poor model scale adaptability, a lack of intelligent sampling decisions, and static, inflexible databases. Summary of the Invention

[0008] The purpose of this invention is to overcome the above-mentioned shortcomings and provide an uncertainty-based method for three-dimensional reconstruction and intelligent sampling of underground space. This method addresses the problems of insufficient depth of multi-source data fusion, low fidelity due to rigid three-dimensional modeling mesh, and lack of an active sampling feedback mechanism based on uncertainty in existing technologies. It enables accurate characterization of complex strata and pollution distribution, as well as intelligent feedback for sampling decisions.

[0009] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: an uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space, comprising the following steps: Step 1: Construct a multi-source heterogeneous data acquisition system and perform data preprocessing and spatiotemporal registration; Step 2: Multimodal heterogeneous feature extraction and unified spatial mapping based on deep learning; Step 3: Construct a high-fidelity 3D reconstruction volumetric model of underground space based on adaptive octree subdivision; Step 4: Construct a high-fidelity intelligent sampling database containing a Bayesian dynamic update layer; Step 5: Perform uncertainty assessment-driven sampling-assisted decision-making to update the database.

[0010] Preferably, in step 1, data preprocessing and spatiotemporal registration specifically include: acquiring borehole sample data, drilling rig behavior data, remote sensing image data, and environmental monitoring data; for drilling rig behavior and environmental monitoring data with timestamps, using basis functions for time interpolation alignment; for spatial data, unifying the coordinate reference through projection transformation and geometric fine correction; for missing values ​​in the data, constructing a three-dimensional data tensor, and using a low-rank tensor completion algorithm based on kernel norm minimization for data recovery; standardizing the data format units and using cross-validation to verify the accuracy and completeness of the data, ensuring the quality of modeling input.

[0011] Preferably, for spatial data, unifying the coordinate datum through projection transformation and geometric calibration specifically includes: (1) The borehole opening coordinates, remote sensing image positioning, and environmental monitoring station coordinates are transformed to a unified reference coordinate system through projection transformation; (2) Perform geometric fine correction on multispectral remote sensing images and DEM elevation data using ground control points to eliminate imaging geometric distortion; (3) Use bilinear interpolation or cubic convolution to resample the raster and unify the standard grid size of raster data from different sources and with different resolutions; (4) For non-geographic reference data, establish a precise mapping relationship with the corrected geographic coordinates using "bore number-sampling depth" or "time stamp-spatial location" as key fields.

[0012] Preferably, in step 2, the multimodal heterogeneous feature extraction and unified spatial mapping specifically includes: establishing parallel feature extraction channels: extracting physicochemical property features from sample data, corresponding to sample attribute channels, drilling rig behavior channels, remote sensing texture channels, and environmental monitoring channels respectively; calculating the spatial gradient of pressure and the rate of change of torque from drilling rig behavior data to quantify abrupt changes in formation mechanics; extracting texture and spectral features from remote sensing data; performing feature aggregation on environmental monitoring data; constructing a deep learning fusion subsystem to perform cascaded fusion of multi-source features using an adaptive weighting mechanism; extracting and normalizing high-order features based on tensor decomposition; and using a spatially constrained kernel spatial fuzzy clustering algorithm, introducing a Gaussian radial basis kernel function and a spatial neighborhood penalty term to generate a multidimensional attribute mapping matrix containing information on formation structure, pollution distribution, and sampling perturbation.

[0013] Preferably, the drilling rig behavior feature extraction process is as follows: (1) Perform gradient analysis on the pressure curve inside the borehole, calculate the pressure change rate and pressure mutation index, and use them to identify changes in formation resistance; (2) Calculate the torque change rate and energy consumption index for the rotational torque time series to reflect the drilling difficulty and formation cohesion; (3) Collect data on drilling speed and location, calculate average drilling speed, drilling acceleration, and number of drilling stops to reflect the smoothness of the drilling process; (4) Calculate the sampling disturbance index and disturbance degree. for: ; In the formula, For the first Soil disturbance index at each sampling depth location; and These are the drilling pressures at the current moment and the previous moment, respectively; and These are the drilling rig torques at the current moment and the previous moment, respectively; and These are the rated maximum drilling pressure and maximum torque for this drilling rig model, respectively. For the current time window Standard deviation of inward advance speed; This represents the average advance rate within the current time window. These are the weighting coefficients for the pressure, torque, and speed terms, respectively, and they satisfy... ; (5) Use clustering or decision tree methods to classify the drilling status and encode the classification results in One-Hot encoding; (6) Output drilling rig behavior feature vector .

[0014] Preferably, the remote sensing texture feature extraction process is as follows: (1) For multispectral remote sensing images, the gray-level co-occurrence matrix (GLCM) method is applied to extract five types of texture features: contrast, homogeneity, energy, correlation and entropy. Under specified directions, i.e. 0° / 45° / 90° / 135° and distance, the gray-level co-occurrence matrix is ​​calculated by pixel pairs, and the corresponding texture features are derived from it. (2) Use the Sobel operator to perform boundary detection, generate an edge gradient map, and calculate the gradient magnitude and gradient direction; (3) Calculate spectral indices based on multispectral images to reflect soil cover and moisture content characteristics; (4) Integrate the features obtained in steps (1) to (3) to output the remote sensing texture feature vector. .

[0015] Preferably, step 3 specifically includes: constructing a three-dimensional spatial reconstruction domain and loading DEM topography, stratigraphic interfaces, and groundwater levels as macroscopic constraints; performing adaptive octree subdivision: calculating the pollution concentration gradient norm, stratigraphic interface curvature, and drilling resistance mutation index within the voxels; if the above indicators exceed a preset threshold, recursively splitting the voxels into sub-voxels until the minimum resolution is reached, thereby generating a non-uniform voxel mesh; constructing an anisotropic three-dimensional reconstruction field: calculating the experimental variogram to identify the main control direction of pollutant diffusion, and constructing an anisotropic search ellipsoid with the major axis along the main control direction; performing voxel assignment: following the hard data locking principle, sequential Gaussian simulation is used for continuous attributes, and sequential indicator simulation is used for discrete attributes to generate a high-fidelity underground space three-dimensional reconstruction voxel model.

[0016] Preferably, constructing a three-dimensional spatial reconstruction domain and loading DEM topography, stratigraphic interfaces, and groundwater level as macroscopic constraints specifically includes: Step 3-1-1: Based on the spatial distribution range of the multi-source heterogeneous data obtained in Step 1, define the physical boundary of the 3D model; traverse all borehole coordinates. And the coverage area of ​​remote sensing images, calculate spatial extrema:

[0017]

[0018] ; in A boundary buffer distance is used to ensure that all sampling points are included within the reconstruction domain; a three-dimensional spatial reconstruction domain is constructed. ; Step 3-1-2: In the reconstructed domain An initial voxel mesh is built internally, serving as the root node of the octree data structure, and the initial voxel side length is set. The domain will be reconstructed. Discretized into a regularly arranged initial volume element At this point, all elements have not yet been subdivided, and their attribute values ​​are empty. Step 3-1-3: Map the known macroscopic geological and hydrological information to the initial set of voxels, which will serve as hard constraints for subsequent subdivisions; the rules for macroscopic geological constraints are as follows: (1) Terrain surface constraints: DEM elevation data is introduced, and the center of the volume elements is determined. z Voxel elements with coordinates above the ground surface elevation are labeled as "air voxels" and removed in subsequent calculations, retaining only "underground voxels"; among which, , in For the body center z coordinate, DEM surface elevation; (2) Stratigraphic interface constraints: Based on the stratigraphic boundaries revealed by boreholes, a simplified stratigraphic trend surface is constructed for each initial subsurface element. The corresponding stratigraphic position is determined based on the location of its center coordinates, and an initial stratigraphic lithology label is assigned. ; (3) Groundwater level constraint: Introduce groundwater level monitoring data , will the body center The voxel marker is called the "vadose zone". The voxel is labeled as "saturation band".

[0019] Preferably, step 4 specifically includes: constructing a multi-layer database architecture, including a metadata layer, a basic structure layer, a voxel index layer, a high-fidelity attribute voxel storage layer, an uncertainty assessment layer, a Bayesian dynamic update layer, and a Web service interface and visualization interaction layer; wherein, the voxel index layer uses space-filling curves to encode and index the non-uniform grid generated by the adaptive octree; the high-fidelity attribute voxel storage layer is used to store the attribute value of each voxel and its corresponding uncertainty parameter; the Bayesian dynamic update layer is used to store prior probability distribution parameters and provide a posterior update interface, using conjugate prior properties to support model parameter correction and version backtracking based on newly sampled data; and establishing a Web service interface to support remote querying and visualization analysis of uncertainty convergence trends.

[0020] Preferably, step 5 specifically includes: (1) Uncertainty assessment: For each voxel, calculate the data variance, soft classification entropy value, and heterogeneous data conflict degree, and generate a weighted comprehensive uncertainty index; (2) Intelligent decision-making: Based on the generated uncertainty heat field, use adaptive thresholds to divide the region into high-risk blind spot areas, medium-risk areas of concern, and low-risk areas of certainty; for the high-risk blind spot areas, generate a digital sampling task book containing recommended coordinates and adaptive working condition parameters; (3) Dynamic update: After performing on-site sampling, use the Bayesian posterior update formula to correct the attribute mean and variance of local voxels based on the accuracy of the new data, reconstruct the uncertainty field, and generate a new version of the database; the dynamic update specifically includes: The on-site drilling equipment receives sampling instructions, executes sampling operations, and transmits actual drilling parameters back in real time, while also combining this data with newly obtained sample concentration data from laboratory analysis. This constitutes an incremental dataset. ;in, This indicates the observed concentration value at the newly added sampling point; Then, the Bayesian update engine of the database is started to correct the voxel attributes of the affected region. The specific process is as follows: (1) Likelihood probability calculation: Assuming that the observation error follows a normal distribution, based on the newly observed concentration Compared with the current model predictions The difference between them constructs a likelihood function, where The likelihood function represents the predicted concentration value at the corresponding volumetric position, expressed as: ; In the formula, This represents the variance of the measurement error in the newly added observation data; (2) Posterior distribution correction: The distribution of voxel attributes is updated using the conjugate prior property. Let the prior mean of the current voxel attribute be... The prior variance is Then, when introducing new observation data Afterwards, the posterior variance of the updated voxel attributes with posterior mean Calculated separately as follows: ; ; In the formula, This represents the average value of the updated element attributes. This represents the variance of the updated voxel attribute; Based on the updated volume element property distribution, the three-component uncertainty index is recalculated to generate an updated uncertainty thermofield; Generate a new version of the database. For multiple rounds of sampling, repeat the above steps to gradually evolve the database version until the uncertainty index converges or the predetermined sampling completion rate is reached, forming the final version of the high-fidelity sampling database.

[0021] Beneficial effects of this invention: (1) This invention innovatively introduces drilling rig behavior data to achieve "drilling-while sensing". Through deep integration with multi-source heterogeneous data such as remote sensing and samples, it accurately identifies thin interlayers, boulders and abrupt geological change interfaces, significantly improving the construction accuracy of geological models for complex sites. (2) This invention employs adaptive octree subdivision technology to dynamically adjust the grid resolution based on pollution gradient and geological characteristics. This achieves "on-demand densification" of key areas, effectively resolving the contradiction between accuracy and computational efficiency in traditional uniform grid models; (3) The present invention constructs a high-fidelity three-dimensional voxel model based on anisotropic constraints and hybrid simulation, which not only accurately restores the spatial variability of the site, but also realizes true three-dimensional, high-fidelity digital twin reconstruction of the geological structure and pollution plume morphology. (4) This invention establishes a three-component uncertainty assessment system consisting of variance, entropy and conflict degree to quantify the cognitive blind spot of the model, thereby guiding the decision-making of high-priority sampling points and adaptive drilling parameters, effectively avoiding blind construction and cross-contamination. (5) This invention constructs a database containing a Bayesian update layer and uses a posterior probability algorithm to receive new data in real time and correct the model. It realizes the transformation of the database from "static storage" to "dynamic evolution" and provides a convergent digital twin foundation for full lifecycle management.

[0022] (6) This invention solves the problems of insufficient depth of multi-source data fusion, low fidelity due to rigid three-dimensional modeling mesh, and lack of active sampling feedback mechanism based on uncertainty in the prior art, and realizes accurate characterization of complex strata and pollution distribution and intelligent feedback of sampling decision. Attached Figure Description

[0023] Figure 1 This is a flowchart illustrating the overall process framework of an uncertainty-based method for three-dimensional reconstruction and intelligent sampling of underground space provided in this embodiment of the invention. Figure 2 This is a detailed flowchart illustrating the feature extraction and deep fusion process of multi-source heterogeneous data; Figure 3 This is a schematic diagram of the logical flow of constructing a high-fidelity three-dimensional reconstruction volumetric model of underground space based on adaptive octree subdivision; Figure 4 This is a schematic diagram of a high-fidelity underground space 3D reconstruction voxel model and attribute binding structure based on adaptive octree subdivision, where the right side shows a single voxel bound to a multi-dimensional attribute vector; Figure 5 This is a schematic diagram of a high-fidelity intelligent sampling database architecture design; Figure 6 It is a logical block diagram of sampling-assisted decision-making and dynamic database updates based on uncertainty assessment. Detailed Implementation

[0024] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0025] Example 1: This embodiment of the invention provides an uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space, referring to... Figure 1 , Figure 1 This is a flowchart illustrating the first embodiment of an uncertainty-based method for three-dimensional reconstruction and intelligent sampling of underground space.

[0026] In this embodiment, the underground space three-dimensional reconstruction and intelligent sampling method includes the following steps: Step 1: Construct a multi-source heterogeneous data acquisition system and perform data preprocessing and spatiotemporal registration; Step 2: Multimodal heterogeneous feature extraction and unified spatial mapping based on deep learning; Step 3: Construct a high-fidelity 3D reconstruction volumetric model of underground space based on adaptive octree subdivision; Step 4: Construct a high-fidelity intelligent sampling database containing a Bayesian dynamic update layer; Step 5: Perform uncertainty assessment-driven sampling-assisted decision-making to update the database.

[0027] Step 1 includes the following sub-steps: Step 1-1: Acquire borehole sample data, drilling rig behavior data, remote sensing image data, and environmental monitoring data. Borehole sample data includes: soil physicochemical indicators (moisture content, density, pH, etc.), heavy metal concentrations (Pb, Cd, Cu, Zn, etc.), organic pollutant concentrations (PAH, PCB, etc.), borehole depth, borehole diameter, and borehole coordinates. Drilling rig behavior data includes: real-time monitoring data such as borehole pressure curves, rotational torque, drilling speed, borehole depth, and vibration acceleration, stored in a time-series database. Remote sensing image data includes: multispectral raster imagery, digital elevation model (DEM), and land cover type. Environmental monitoring data includes: time-series data such as temperature, relative humidity, groundwater level, precipitation, air pressure, and wind speed.

[0028] Step 1-2: Perform spatial correction on all.

[0029] Steps 1-3: For data with timestamps (drilling rig behavior, environmental monitoring), use basis functions (polynomial basis functions, radial basis functions, or B-spline basis functions) for time interpolation to align data at different time scales (such as daily, hourly, and minute-level data) to a unified time node; for non-time series data (such as sample analysis results), configure them as the timestamp of the sample sampling time.

[0030] Steps 1-4: Handling missing and outlier values.

[0031] Steps 1-5: Convert all data into a unified data format; for numerical data, use a unified unit system; for categorical data, standardize the coding (e.g., use national standard coding for lithological classification).

[0032] Steps 1-6: Use methods such as cross-validation and comparison with authoritative reference datasets to verify the accuracy, completeness and consistency of the data, and ensure that the data quality meets the requirements of subsequent modeling (accuracy ≥ 95%, completeness ≥ 95%).

[0033] In step 1-2, the specific process of spatial correction is as follows: (1) All acquired geospatial data (including borehole opening coordinates, remote sensing image positioning, and environmental monitoring station coordinates) are transformed into a unified reference coordinate system through projection transformation; (2) For multispectral remote sensing images and DEM elevation data, geometric fine correction is performed using ground control points to eliminate geometric distortions during the imaging process; (3) Use bilinear interpolation or cubic convolution to resample the raster data from different sources and at different resolutions to uniformly adjust them to the standard grid size, thereby achieving pixel-level spatial alignment of multi-source images. (4) For non-geographic reference data that does not have direct geographic coordinates (such as soil physicochemical indicators and drilling rig torque logs), use "hole number-sampling depth" or "time stamp-spatial location" as key fields to establish a precise mapping relationship with the corrected geographic coordinates and complete the positioning of attribute information in three-dimensional space.

[0034] In steps 1-4, the specific rules for handling missing and outlier values ​​are as follows: (1) For discrete missing data, the K-nearest neighbor algorithm or a correlation-based prediction model is used to fill in the missing data; (2) For systematic missing data (such as no data for a certain measuring point for a long time), the low-rank tensor decomposition method is used to fill the missing data. (3) For outliers, use the quartile method or the 3σ criterion to identify them, and use the median or neighboring values ​​to replace or remove outliers.

[0035] Step 2 includes the following sub-steps: Step 2-1: Establish multiple parallel feature extraction channels, corresponding to the sample attribute channel, drilling rig behavior channel, remote sensing texture channel, and environmental monitoring channel, respectively.

[0036] Step 2-2: Establish a fusion subsystem to achieve adaptive weighting and cascade fusion of multi-source features.

[0037] Steps 2-3: Perform spatial constraint-based kernel spatial fuzzy clustering Steps 2-4: Fuse the feature vectors All sample points (corresponding to each borehole) are arranged in matrix form, forming a multidimensional attribute mapping matrix. This matrix will serve as the basic input for subsequent 3D voxel modeling and uncertainty assessment.

[0038] Step 2-1 includes the following sub-steps: Step 2-1-1: Extract sample attribute features; Step 2-1-2: Extract drilling rig behavior features; Step 2-1-3: Extract remote sensing texture features; Step 2-1-4: Extract environmental monitoring features.

[0039] In step 2-1-1, the sample attribute feature extraction process is as follows: (1) Detect outliers in the physicochemical properties and contaminant concentrations of borehole samples, and use the IQR method or a modified Z-score method to identify outlier samples and mark or remove them. (2) Vectorize and encode sample attributes. For continuous numerical attributes (such as heavy metal concentration), use them directly as features; for categorical attributes (such as lithology), use One-Hot encoding or sequence encoding; for missing values, use the sample mean or KNN to fill in the missing values ​​before encoding. (3) Perform standardization mapping to map each dimension of features to a distribution with a mean of 0 and a standard deviation of 1, thereby eliminating the influence of dimensional differences; (4) Optionally, principal component analysis (PCA) can be used for dimensionality reduction, retaining principal components with a cumulative explained variance contribution rate ≥ 90%; (5) Output sample attribute feature vector .

[0040] In step 2-1-2, the drilling rig behavior feature extraction process is as follows: (1) Perform gradient analysis on the pressure curve inside the borehole, calculate the pressure change rate and pressure mutation index, and use them to identify changes in formation resistance; (2) Calculate the torque change rate and energy consumption index for the rotational torque time series to reflect the drilling difficulty and formation cohesion; (3) Collect data on drilling speed and location, calculate average drilling speed, drilling acceleration, number of drilling stops, etc., to reflect the smoothness of the drilling process; (4) Calculate the sampling disturbance index and disturbance degree. for:

[0041] In the formula, For the first Soil disturbance index (dimensionless, range 0-1) at each sampling depth location; and These are the drilling pressures at the current and previous moments, respectively (unit: kN). and The drilling rig torque (unit: N·m) is shown for the current moment and the previous moment, respectively. and These are the rated maximum drilling pressure and maximum torque for this drilling rig model (used for normalization to eliminate the influence of dimensions). For the current time window The standard deviation of the internal feed rate (reflects drill bit runout or stuck drill bit phenomenon); This represents the average advance rate within the current time window. These are the weighting coefficients for the pressure, torque, and speed terms, respectively, and they satisfy... .

[0042] (5) Use clustering or decision tree methods to classify drilling status (such as "smooth drilling", "local stuck drill", "forced drilling", etc.), and encode the classification results in One-Hot encoding; (6) Output drilling rig behavior feature vector .

[0043] In step 2-1-3, the remote sensing texture feature extraction process is as follows: (1) Apply the Gray-Level Co-occurrence Matrix (GLCM) method to multispectral remote sensing images to extract five types of texture features: contrast, homogeneity, energy (second moment of angle), correlation, and entropy; specifically, under specified directions (0° / 45° / 90° / 135°) and distances, calculate the gray-level co-occurrence matrix of pixel pairs, thereby deriving the texture features; (2) Boundary detection is performed using the Sobel operator to generate an edge gradient map. Specifically, 3×3 horizontal and vertical gradient operators are convolved with the image to calculate the gradient magnitude. and gradient direction Enhance the accuracy of boundary identification in polluted areas; (3) Calculate spectral indices, such as vegetation indices. Improved soil-adjusted vegetation index (MSAVI) reflects soil cover and moisture content characteristics; (4) Output remote sensing texture feature vector .

[0044] In step 2-1-4, the environmental monitoring feature extraction process is as follows: (1) Use methods such as moving average or wavelet filtering to denoise high-frequency noise and extract trend and periodic features; (2) Perform feature aggregation and calculate the statistical characteristics of each environmental parameter during the sampling period, including mean, standard deviation, maximum value, minimum value, peak value, slope, etc. (3) Use the 3σ criterion or Isolation Forest algorithm to detect anomalies and identify the time periods when extreme weather or groundwater anomalies occur; (4) Use standardized mapping to map the features to the [0,1] interval; (5) Output environmental monitoring feature vector .

[0045] Step 2-2 includes the following sub-steps: Step 2-2-1: Combine the four feature vectors , , and By concatenating the components along the dimensional direction, a preliminary fused vector is obtained. .

[0046] Step 2-2-2: Using cross-validation, the data is divided into training and validation sets. By minimizing the prediction error on the validation set, the optimal fusion weights of the four feature vectors are automatically learned. Specifically, gradient descent or genetic algorithm is used to optimize the weights so that the prediction accuracy of the fusion vector for multiple target variables (such as pollution concentration and formation type) is maximized.

[0047] Step 2-2-3: Considering the potential feature redundancy and high-dimensional sparsity of the concatenated fusion vector, tensor decomposition (such as Tucker decomposition or CP decomposition) is performed on the weighted concatenated fusion vector to further extract higher-order interaction features; finally, L2 norm normalization is applied to obtain the final fusion feature vector. .

[0048] In steps 2-3, for the fused feature vector Given the nonlinear distribution and complex topological structure exhibited in the feature space, a spatially constrained kernel fuzzy C-means clustering algorithm is introduced to map the data to a high-dimensional Hilbert space for enhanced separability. The specific sub-steps are as follows: Step 2-3-1: Introduce Gaussian radial basis functions as kernel functions This transforms the Euclidean distance in the original feature space into a nonlinear distance in the kernel space. The kernel function is defined as:

[0049] In the formula, For the first The fused feature vector of each sampling point For the first Feature vectors of cluster centers Let be the kernel bandwidth parameter. At this point, the distance between the sample point and the cluster center in the kernel space is... .

[0050] Step 2-3-2: To overcome the shortcomings of traditional clustering algorithms, such as sensitivity to noise and neglect of geological spatial continuity, the following objective function is constructed:

[0051] In the formula, This represents the number of cluster categories (corresponding to the number of potential geological / contamination patterns, such as: clean sand, lightly contaminated silt, heavily contaminated clay, etc.). The membership matrix, For the first The sample belongs to the first The probability of a class (soft classification) satisfies ; This is a fuzzy weighted index, typically 2; These are the spatial constraint weighting coefficients; For spatial neighborhood penalty terms, ,in For sample points Spatial neighborhood window.

[0052] Step 2-3-3: Iteratively optimize the objective function using the Lagrange multiplier method until the convergence condition is met. The iterative update formula is as follows: Membership update formula:

[0053] Cluster center update formula:

[0054] Steps 2-3-4: After iterative convergence, output the optimal membership matrix. and cluster center .matrix Each element in This will be directly used to calculate the soft classification entropy value in subsequent steps. Cluster center This represents a typical "geology-pollution-disturbance" coupling pattern within the site and is stored in the database as prior knowledge.

[0055] Step 3 includes the following sub-steps: Step 3-1: Construct the 3D spatial reconstruction domain and load macroscopic constraints; Step 3-2: Adaptive octree subdivision and non-uniform volume element generation; Step 3-3: Construct an anisotropic 3D reconstruction field; Steps 3-4: Assignment and attribute mapping of multi-source data to voxels.

[0056] Step 3-1 includes the following sub-steps: Step 3-1-1: Based on the spatial distribution range of the multi-source heterogeneous data obtained in Step 1, define the physical boundaries of the 3D model. Traverse all borehole coordinates. And the coverage area of ​​remote sensing images, calculate spatial extrema:

[0057]

[0058] ; in The boundary buffer distance is used to ensure that all sampling points are included within the reconstruction domain. Construct the 3D spatial reconstruction domain.

[0059] Step 3-1-2: In the reconstructed domain An initial voxel mesh is built internally, serving as the root node of the octree data structure, and the initial voxel side length is set. The domain will be reconstructed. Discretized into a regularly arranged initial volume element At this point, all elements have not yet been subdivided, and their attribute values ​​are empty.

[0060] Step 3-1-3: Map the known macro-geological and hydrological information to the initial set of voxels as hard constraints for subsequent subdivision.

[0061] In step 3-1-3, the rules for macroscopic geological constraints are as follows: (1) Terrain surface constraints: DEM elevation data is introduced, and the center of the volume elements is determined. z Volume elements whose coordinates are higher than the ground surface elevation ( The elements marked as "air elements" will be removed in subsequent calculations, leaving only "underground elements"; (2) Stratigraphic interface constraints: Based on the stratigraphic boundaries revealed by boreholes, a simplified stratigraphic trend surface is constructed for each initial subsurface element. The corresponding stratigraphic position is determined based on the location of its center coordinates, and an initial stratigraphic lithology label is assigned. ; (3) Groundwater level constraint: Introduce groundwater level monitoring data , will the body center The voxel marker is called the "vadose zone". The voxel is labeled as "saturation band".

[0062] Step 3-2 includes the following sub-steps: Step 3-2-1: Based on the fusion features and macroscopic geological constraints generated in Step 2, calculate the three core indicators driving mesh subdivision, as follows: (1) Pollution concentration gradient norm This indicator is used to capture the diffusion front and core boundary of a contamination plume.

[0063] (2) Curvature of stratigraphic interfaces : Calculates the divergence of the stratigraphic normal vector. This index is used to capture complex geological interfaces such as fold structures and fault fracture zones.

[0064] (3) Characteristics of sudden changes in drilling resistance The second derivative of drilling rig torque and pressure is used to identify anomalous jump points. This indicator is used to detect underground boulders, cavities, or hard interlayers.

[0065] Step 3-2-2: Initialize the volume elements Perform recursive subdivision, with the following subdivision logic: (1) Fine division rules: If volume elements satisfy (High gradient), or (High curvature), or (Sudden change in resistance) will change the current volume element Perform a standard octree split, i.e., along The system is orthogonally divided in three directions to generate eight sub-elements with their size halved. The above judgment is then performed on the generated sub-elements until the preset minimum resolution is reached, so as to achieve high-fidelity capture of pollution boundaries and geological details.

[0066] (2) Coarse division retention rule: If voxels If all indicators are below the threshold and the internal attribute variance is extremely small, then the subdivision stops and the current voxel is retained as a leaf node. Such voxels maintain a large geometric size to reduce the amount of data stored and computational redundancy in homogeneous regions.

[0067] Step 3-2-3: After the above recursive process, the final set of non-uniform volume elements is generated. .

[0068] Step 3-3 includes the following sub-steps: Step 3-3-1: For each volume element The fused feature vector is calculated around it using a sliding window. Spatial autocorrelation coefficients, using Moran's I or Pearson correlation coefficients; Step 3-3-2: Identify the preferred direction of pollution diffusion (such as along groundwater flow or along fault zones) through principal component analysis or variogram analysis. Step 3-3-3: Based on spatial correlation and dominant direction, for each volume element Construct the search neighborhood ellipsoid The major semi-axis of the ellipsoid is set to extend along the main control direction to capture long-range correlations along the flow direction; the minor semi-axis is set to contract perpendicular to the main control direction to shield against interference from lateral irrelevant data. This ellipsoid is used to filter out valid sample points with physical correlations during subsequent data assignment.

[0069] Steps 3-4 include the following sub-steps: Step 3-4-1: For voxels that contain actual sampling points in space (drill hole sampling or drilling rig trajectory), mark them as "hard data voxels" and lock them in the subsequent simulation process without interpolation correction. The assignment rules are as follows: (1) Continuous attributes (such as concentration, pressure): If a volume element contains multiple observation data, the representative value of the volume element is calculated using the volume weighted average method; (2) Discrete attributes (such as lithology): If a volume element contains multiple lithological records, the mode principle is adopted, and the lithology category with the highest frequency is assigned to the volume element.

[0070] Step 3-4-2: For unobserved volume elements, use the anisotropic search neighborhood ellipsoid constructed in Step 4.4. Determine neighborhood samples and employ differentiated interpolation strategies for different attribute types: (1) Use sequential Gaussian simulation with continuous attributes for continuous variables such as pollutant concentration and formation porosity; (2) Sequential indicator simulation using discrete attributes for classification variables such as stratigraphic lithology and adverse geological body labels; (3) For mechanically derived properties such as sampling disturbance degree and drilling specific work, an improved inverse distance weighting method or a neighborhood weighting method based on ellipsoidal distance is adopted.

[0071] Step 3-4-3: Traverse all voxels, complete the assignment, generate a high-fidelity three-dimensional reconstruction voxel model of underground space, and form an attribute mapping table.

[0072] Step 4 includes the following sub-steps: Step 4-1: Based on the voxel attributes and uncertainty assessment results, construct a multi-layer database architecture, including a metadata layer, a basic structure layer, a voxel index layer, a high-fidelity attribute voxel storage layer, an uncertainty assessment layer, a Bayesian dynamic update layer, and a Web service interface and visualization interaction layer. Step 4-2: Implement the database visualization interface to display the 3D reconstructed voxel model, the heat map of pollution concentration distribution, the uncertainty field distribution map, the sampling points and borehole trajectories, and version comparison display, etc. Step 4-3: Establish a Web service interface and data sharing mechanism, provide a RESTful API interface, support remote database query, and support outputting sampling task books and model data in JSON format; provide a data comparison interface between versions to analyze the convergence trend of uncertainty; develop a WebGL-based 3D visualization front-end to provide core views such as 3D volumetric models and uncertain thermal fields.

[0073] In step 4-1, the construction steps of the voxel index layer, the high-fidelity attribute voxel storage layer, and the Bayesian dynamic update layer are as follows: (1) Construct a voxel index layer adapted to non-uniform grids. For the adaptive octree non-uniform grids generated in step 3, the voxel index layer uses Morton code or hash mapping technology to uniquely encode each voxel, and establishes a millisecond-level retrieval index of "spatial coordinates-voxel ID" to solve the addressing problem of massive multi-scale voxels.

[0074] (2) Construct a high-fidelity attribute voxel storage layer. This layer stores the attribute mapping table generated in step 3 into the corresponding voxels. Unlike traditional databases that only store single scalar values, this embodiment also stores the statistical parameters of voxel attributes. Specifically, for each voxel, information such as the mean and variance of its attributes is recorded, along with the data source identifier and timestamp of the voxel, providing a data foundation for subsequent uncertainty calculations.

[0075] (3) Construct a Bayesian dynamic update layer. This layer serves as the computational core of the database and incorporates a built-in Bayesian inference engine. It stores the prior probability distribution of all voxels (i.e., the initial state of modeling in step 3). This layer reserves a posterior update interface. When new sampled data (incremental data) is transmitted back from the field, the conjugate prior property is used to automatically calculate the influence weight of the new data on the surrounding voxels and update the mean and variance parameters of the affected voxels in real time, without having to recalculate the entire model, thereby enabling rapid iteration of the database version (e.g., from v1.0 to v1.1).

[0076] Step 5 includes the following sub-steps: Step 5-1: Three-component uncertainty assessment based on multi-source information fusion; Step 5-2: Uncertainty-driven sampling priority allocation and intelligent decision-making; Step 5-3: On-site sampling execution and Bayesian dynamic update.

[0077] Step 5-1 includes the following sub-steps: Step 5-1-1: Calculate the data variance field. For each 3D volume element... This quantifies the dispersion of multi-source observation data. It searches for the effective set of observations within the statistical volume element and the anisotropic ellipsoid. Calculate the normalized variance index This indicator reflects the physical dispersion of the data. The larger the variance, the more drastic the fluctuations in the data in that region, and the worse the representativeness of the existing observations.

[0078] Step 5-1-2: Calculate the soft clustering entropy value. Based on the membership matrix output by kernel space fuzzy clustering. To quantify the fuzziness of geological model classification. To calculate information entropy. ,in The number of cluster categories. For body element Belongs to the The probability of a geological-like model. When When the time is specified, it indicates that the volume element is at the transition interface between different geological types (such as clay and sand) or different pollution levels.

[0079] Step 5-1-3: Calculate the conflict degree of heterogeneous data. Calculate the volume element. Concentration analysis at the borehole Remote sensing inversion concentration Geostatistical interpolation concentration Maximum normalization deviation between:

[0080] Step 5-1-4: Calculate the final uncertainty index using a weighted summation method, where the weights are adaptively optimized using a cross-validation method, with the goal of minimizing the mean squared error of the predictions on the validation set.

[0081] Step 5-1-5: Generate a three-dimensional continuous thermodynamic field with an uncertain distribution.

[0082] Step 5-2 includes the following sub-steps: Step 5-2-1: Calculate the quartiles based on the statistical histogram of the overall uncertainty index. Define a dynamic segmentation threshold: Low threshold:

[0083] High threshold:

[0084] Step 5-2-2: Divide the region to be sampled into three priority levels based on the threshold, and match the corresponding sampling density and sampling condition parameters, etc. The specific division rules are as follows: (1) Priority I (High-risk blind spot): meets the requirements Encrypted sampling is performed in the designated area, and the feed rate is kept low.

[0085] (2) Priority II (Medium Concern Area) is satisfied The area is sampled using standard grids, with medium-speed advance rates, etc.

[0086] (3) Priority III (Low-risk confidence zone): meets the requirements Sparse verification sampling is performed in the region, and high-speed advance speed is achieved.

[0087] Step 5-2-3: Based on the priority division of the sampling area, generate a digital sampling task book containing recommended point coordinates, predicted stratum lithology, and recommended working condition parameters, and export it as a GIS visualization file and a control instruction file that can be read by the drilling rig.

[0088] Step 5-3 includes the following sub-steps: Step 5-3-1: The on-site drilling equipment receives the instruction to perform sampling, transmits actual drilling parameters (pressure, torque logs, etc.) in real time, and combines them with the new sample concentration data obtained from laboratory analysis. This constitutes an incremental dataset. .

[0089] Step 5-3-2: Start the Bayesian update engine of the database to correct the voxel attributes of the affected region. The specific process is as follows: (1) Likelihood probability calculation: Assuming that the observation error follows a normal distribution, based on the newly observed concentration Compared with the current model predictions The difference between them constructs a likelihood function, where The likelihood function represents the predicted concentration value at the corresponding volumetric position, expressed as: ; In the formula, This represents the variance of the measurement error in the newly added observation data.

[0090] (2) Posterior distribution correction: The distribution of voxel attributes is updated using the conjugate prior property. Let the prior mean of the current voxel attribute be... The prior variance is Then, when introducing new observation data Afterwards, the posterior variance of the updated voxel attributes with posterior mean Calculated separately as follows: ; ; In the formula, This represents the average value of the updated element attributes. This represents the variance of the updated voxel attribute; Step 5-3-3: Based on the updated volume element attribute distribution, recalculate the three-component uncertainty index to generate an updated uncertainty thermofield.

[0091] Step 5-3-4: Generate a new version of the database (e.g., upgrade from v1.0 to v1.1). For multi-round sampling, repeat the above steps, gradually evolving the database version (v1.1 → v1.2 → ...) until the uncertainty index converges. (or until the predetermined sampling completion rate is reached, forming the final version of the high-fidelity sampling database.)

[0092] Example 2: This example selects a factory contaminated site with complex geological conditions as the application object. This site has an interlayered sand-clay structure and a heterogeneously distributed trichloroethylene contamination plume. The specific implementation process of the intelligent sampling and data reconstruction uncertainty method for underground space described in this invention in actual engineering is as follows: Step 1: Construct a multi-source heterogeneous data acquisition system and perform data preprocessing. First, exploratory boreholes (labeled ZK-01 to ZK-05) were drilled at the site to acquire borehole sample data, including TCE concentration, soil pH, and moisture content sampled every 0.5 meters. Simultaneously, sensors mounted on a hydraulic push drill rig recorded real-time drilling rig behavior data at a frequency of 10Hz, including borehole pressure curves, rotational torque, and drilling speed. In addition, multispectral remote sensing imagery of the site acquired by UAV and recent groundwater level monitoring data were combined. For the above multi-source data, basis functions were used for time interpolation alignment to unify high-frequency drilling rig behavior data and low-frequency environmental monitoring data to a standard time axis. Projection transformation was used to unify the spatial coordinate reference. Missing values ​​in the borehole data were restored using a low-rank tensor completion algorithm to ensure the integrity of the input data. The data format units were standardized, and cross-validation was used to verify the accuracy and completeness of the data, ensuring the quality of the modeling input.

[0093] Step 2: Feature Fusion and Spatial Mapping Based on Deep Learning. Parallel feature extraction channels are established to extract the physicochemical properties of the samples, the spatial gradient of pressure and torque change rate of the drilling rig, and the texture spectral features of remote sensing images and environmental monitoring data. A deep learning fusion subsystem is used to adaptively weight and concatenate these features, extracting and normalizing higher-order features based on tensor decomposition. A spatially constrained kernel spatial fuzzy clustering algorithm is employed, introducing a Gaussian radial basis kernel function and a spatial neighborhood penalty term to generate a multidimensional attribute mapping matrix containing stratigraphic structure, pollution distribution, sampling disturbance information, and environmental monitoring data. For example, at a borehole depth of 6.5 meters in ZK-03, the algorithm identifies the feature combination of "high torque change rate + low drilling speed," which is mapped to a geological pattern of "calcareous nodule aquitard" through cluster analysis, thus generating a multidimensional attribute mapping matrix containing stratigraphic structure and pollution distribution information.

[0094] Step 3: Construct a high-fidelity 3D reconstruction voxel model of underground space based on adaptive octree subdivision. A 3D spatial reconstruction domain covering the underground space of the site is constructed, with DEM topography and groundwater level loaded as macroscopic constraints. Adaptive octree subdivision is performed: In the homogeneous fill area where ZK-01 is located, because both the pollution concentration gradient norm and the curvature of the stratum interface are below the preset threshold, the algorithm retains large-size coarse voxels with a side length of 4 meters to save computational power; while in the 8-12 meter depth range near ZK-02, the TCE concentration gradient norm is detected and the drilling resistance mutation index exceeds the standard, so the algorithm automatically recursively splits the voxels to the minimum resolution (e.g., 0.125 meters), generating a non-uniform high-density mesh to accurately characterize the pollution plume boundary and hard interlayer details. Subsequently, the groundwater flow direction is identified as the main controlling direction, an anisotropic search ellipsoid is constructed, and sequential Gaussian simulation and sequential indicator simulation are performed on unsampled voxels to generate a high-fidelity 3D numerical model.

[0095] Step 4: Construct a high-fidelity sampling database with a Bayesian dynamic update layer. Based on the generated 3D reconstruction model, the system goes beyond static storage and constructs a multi-layered database architecture that supports dynamic evolution. Specifically, this database includes a voxel index layer for storing adaptive octree non-uniform grid codes, enabling efficient retrieval of massive micro-voxels; a high-fidelity attribute voxel storage layer for recording the lithology category and TCE concentration prediction value of each voxel; and the core Bayesian dynamic update layer and uncertainty assessment layer. The Bayesian dynamic update layer abandons the traditional approach of storing only a single numerical value and innovatively stores the prior probability distribution parameters of each voxel attribute, reserving a mathematical interface for subsequently accepting new data and performing posterior corrections. Furthermore, the system deploys a web service interface and a visualization interaction layer, supporting remote users to view 3D cross-sections of underground space and render uncertainty fields in real time.

[0096] Step 5: Execute uncertainty assessment-driven sampling-assisted decision-making and dynamic database updates. The system first calculates the data variance (physical discreteness), soft classification entropy (geological ambiguity), and heterogeneous data conflict degree of all voxels based on parameters in the database, generating a weighted uncertainty index and further generating an uncertainty thermal field. The calculation reveals that the uncertainty index of the central area of ​​the site (coordinates X:120, Y:85) exceeds the high-risk threshold. The system classifies it as sampling "Priority I (High-Risk Blind Spot Filling)" and automatically generates a digital sampling task sheet, recommending the placement of verification borehole ZK-06 at this location. It also predicts a hard interlayer at a depth of 5 meters in this area and instructs the drilling rig to reduce its rotation speed by 30%. After on-site drilling is completed, the newly acquired sample and drilling rig data are transmitted back via an interface, triggering the database's Bayesian update engine. The engine performs a convolution operation using the likelihood function of the new data and the prior distribution stored in the database to calculate the posterior probability, thereby correcting the attribute mean of the voxels in this area and converging the variance. Ultimately, the system automatically generated a new version of the database (v1.1), realizing a closed-loop evolution from a static model to a growable digital twin.

[0097] The above embodiments are merely preferred technical solutions of the present invention and should not be considered as limitations on the present invention. The scope of protection of the present invention should be limited to the technical solutions described in the claims, including equivalent substitutions of the technical features described in the claims. That is, equivalent substitutions and improvements within this scope are also within the scope of protection of the present invention.

Claims

1. An uncertainty method for 3D reconstruction and intelligent sampling of underground space, characterized in that: Includes the following steps: Step 1: Construct a multi-source heterogeneous data acquisition system and perform data preprocessing and spatiotemporal registration; Step 2: Multimodal heterogeneous feature extraction and unified spatial mapping based on deep learning; Step 3: Construct a high-fidelity 3D reconstruction volumetric model of underground space based on adaptive octree subdivision; Step 4: Construct a high-fidelity intelligent sampling database containing a Bayesian dynamic update layer; Step 5: Perform uncertainty assessment-driven sampling-assisted decision-making to update the database.

2. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 1, characterized in that: In step 1, data preprocessing and spatiotemporal registration specifically include: acquiring borehole sample data, drilling rig behavior data, remote sensing image data, and environmental monitoring data; for drilling rig behavior and environmental monitoring data with timestamps, using basis functions for time interpolation alignment; for spatial data, unifying the coordinate reference through projection transformation and geometric fine correction; for missing values ​​in the data, constructing a three-dimensional data tensor and using a low-rank tensor completion algorithm based on kernel norm minimization for data recovery; standardizing the data format units and using cross-validation to verify the accuracy and completeness of the data to ensure the quality of modeling input.

3. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 2, characterized in that: For spatial data, unifying the coordinate datum through projection transformation and geometric calibration specifically includes: (1) The borehole opening coordinates, remote sensing image positioning, and environmental monitoring station coordinates are transformed to a unified reference coordinate system through projection transformation; (2) Perform geometric fine correction on multispectral remote sensing images and DEM elevation data using ground control points to eliminate imaging geometric distortion; (3) Use bilinear interpolation or cubic convolution to resample the raster and unify the standard grid size of raster data from different sources and with different resolutions; (4) For non-geographic reference data, establish a precise mapping relationship with the corrected geographic coordinates using "bore number-sampling depth" or "time stamp-spatial location" as key fields.

4. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 1, characterized in that: In step 2, the multimodal heterogeneous feature extraction and unified spatial mapping specifically include: establishing parallel feature extraction channels: extracting physicochemical property features from sample data, corresponding to sample attribute channels, drilling rig behavior channels, remote sensing texture channels, and environmental monitoring channels respectively; calculating the spatial gradient of pressure and the rate of change of torque from drilling rig behavior data to quantify abrupt changes in formation mechanics; extracting texture and spectral features from remote sensing data; performing feature aggregation on environmental monitoring data; constructing a deep learning fusion subsystem to perform cascaded fusion of multi-source features using an adaptive weighting mechanism; extracting and normalizing high-order features based on tensor decomposition; and using a spatially constrained kernel spatial fuzzy clustering algorithm, introducing a Gaussian radial basis kernel function and a spatial neighborhood penalty term to generate a multidimensional attribute mapping matrix containing information on formation structure, pollution distribution, and sampling perturbation.

5. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 4, characterized in that: The process of extracting drilling rig behavioral features is as follows: (1) Perform gradient analysis on the pressure curve inside the borehole, calculate the pressure change rate and pressure mutation index, and use them to identify changes in formation resistance; (2) Calculate the torque change rate and energy consumption index for the rotational torque time series to reflect the drilling difficulty and formation cohesion; (3) Collect data on drilling speed and location, calculate average drilling speed, drilling acceleration, and number of drilling stops to reflect the smoothness of the drilling process; (4) Calculate the sampling disturbance index and disturbance degree. for: ； In the formula, For the first Soil disturbance index at each sampling depth location; and These are the drilling pressures at the current moment and the previous moment, respectively; and These are the drilling rig torques at the current moment and the previous moment, respectively; and These are the rated maximum drilling pressure and maximum torque for this drilling rig model, respectively. For the current time window Standard deviation of inward advance speed; This represents the average advance rate within the current time window. These are the weighting coefficients for the pressure, torque, and speed terms, respectively, and they satisfy... ; (5) Use clustering or decision tree methods to classify the drilling status and encode the classification results in One-Hot encoding; (6) Output drilling rig behavior feature vector .

6. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 4, characterized in that: The remote sensing texture feature extraction process is as follows: (1) For multispectral remote sensing images, the gray-level co-occurrence matrix (GLCM) method is applied to extract five types of texture features: contrast, homogeneity, energy, correlation and entropy. Under specified directions, i.e. 0° / 45° / 90° / 135° and distance, the gray-level co-occurrence matrix is ​​calculated by pixel pairs, and the corresponding texture features are derived from it. (2) Use the Sobel operator to perform boundary detection, generate an edge gradient map, and calculate the gradient magnitude and gradient direction; (3) Calculate spectral indices based on multispectral images to reflect soil cover and moisture content characteristics; (4) Integrate the features obtained in steps (1) to (3) to output the remote sensing texture feature vector. .

7. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 1, characterized in that: Step 3 specifically includes: constructing a three-dimensional spatial reconstruction domain and loading DEM topography, stratigraphic interfaces, and groundwater levels as macroscopic constraints; performing adaptive octree subdivision: calculating the pollution concentration gradient norm, stratigraphic interface curvature, and drilling resistance mutation index within the voxels; if the above indicators exceed a preset threshold, recursively splitting the voxels into sub-voxels until the minimum resolution is reached, thereby generating a non-uniform voxel mesh; constructing an anisotropic three-dimensional reconstruction field: calculating the experimental variogram to identify the main control direction of pollutant diffusion, and constructing an anisotropic search ellipsoid with the major axis along the main control direction; performing voxel assignment: following the hard data locking principle, sequential Gaussian simulation is used for continuous attributes, and sequential indicator simulation is used for discrete attributes to generate a high-fidelity underground space three-dimensional reconstruction voxel model.

8. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 7, characterized in that: Constructing a three-dimensional spatial reconstruction domain and loading DEM topography, stratigraphic interfaces, and groundwater level as macroscopic constraints specifically includes: Step 3-1-1: Based on the spatial distribution range of the multi-source heterogeneous data obtained in Step 1, define the physical boundary of the 3D model; traverse all borehole coordinates. And the coverage area of ​​remote sensing images, calculate spatial extrema: ； in A boundary buffer distance is used to ensure that all sampling points are included within the reconstruction domain; a three-dimensional spatial reconstruction domain is constructed. ; Step 3-1-2: In the reconstructed domain An initial voxel mesh is built internally, serving as the root node of the octree data structure, and the initial voxel side length is set. , will reconstruct the domain Discretized into a regularly arranged initial volume element At this point, all elements have not yet been subdivided, and their attribute values ​​are empty. Step 3-1-3: Map the known macroscopic geological and hydrological information to the initial set of voxels, serving as hard constraints for subsequent subdivision; the rules for macroscopic geological constraints are as follows: (1) Terrain surface constraints: DEM elevation data is introduced, and the center of the volume elements is determined. z Voxel elements with coordinates above the ground surface elevation are labeled "air voxels" and removed in subsequent calculations, retaining only "underground voxels"; among which, , in For the body center z coordinate, DEM surface elevation; (2) Stratigraphic interface constraints: Based on the stratigraphic boundaries revealed by boreholes, a simplified stratigraphic trend surface is constructed for each initial subsurface element. The corresponding stratigraphic position is determined based on the location of its center coordinates, and an initial stratigraphic lithology label is assigned. ; (3) Groundwater level constraint: Introduce groundwater level monitoring data , will the body center The voxel marker is called the "vadose zone". The voxel is labeled as "saturation band".

9. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 1, characterized in that: Step 4 specifically includes: constructing a multi-layered database architecture, including a metadata layer, a basic structure layer, a voxel index layer, a high-fidelity attribute voxel storage layer, an uncertainty assessment layer, a Bayesian dynamic update layer, and a Web service interface and visualization interaction layer; wherein, the voxel index layer uses space-filling curves to encode and index the non-uniform grid generated by the adaptive octree; the high-fidelity attribute voxel storage layer is used to store the attribute values ​​of each voxel and their corresponding uncertainty parameters; the Bayesian dynamic update layer is used to store prior probability distribution parameters and provide a posterior update interface, using conjugate prior properties to support model parameter correction and version backtracking based on newly sampled data; and establishing a Web service interface to support remote querying and visualization analysis of uncertainty convergence trends.

10. The uncertainty method for three-dimensional reconstruction and intelligent sampling of underground space according to claim 1, characterized in that: Step 5 specifically includes: (1) Uncertainty assessment: For each voxel, calculate the data variance, soft classification entropy value and heterogeneous data conflict degree, and generate a weighted comprehensive uncertainty index; (2) Intelligent decision-making: Based on the generated uncertainty heat field, use adaptive thresholds to divide the region into high-risk blind spot, medium-risk concern and low-risk confidence areas; for the high-risk blind spot, generate a digital sampling task book containing recommended coordinates and adaptive working condition parameters; (3) Dynamic update: After performing on-site sampling, use the Bayesian posterior update formula to correct the attribute mean and variance of local voxels based on the accuracy of the new data, reconstruct the uncertainty field and generate a new version of the database; Dynamic update specifically includes: The on-site drilling equipment receives sampling instructions, executes sampling operations, and transmits actual drilling parameters back in real time, while also combining this data with newly obtained sample concentration data from laboratory analysis. This constitutes an incremental dataset. ;in, This indicates the observed concentration value at the newly added sampling point; Then, the Bayesian update engine of the database is started to correct the voxel attributes of the affected region. The specific process is as follows: (1) Likelihood probability calculation: Assuming that the observation error follows a normal distribution, based on the newly observed concentration Compared with the current model predictions The difference between them constructs a likelihood function, where The likelihood function represents the predicted concentration value at the corresponding volumetric position, expressed as: ； In the formula, This represents the variance of the measurement error in the newly added observation data; (2) Posterior distribution correction: The distribution of voxel attributes is updated using the conjugate prior property. Let the prior mean of the current voxel attribute be... The prior variance is Then, when introducing new observation data Afterwards, the posterior variance of the updated voxel attributes with posterior mean Calculated separately as follows: ； ； In the formula, This represents the average value of the updated voxel attribute. This represents the updated variance of the voxel attribute; Based on the updated volume element property distribution, the three-component uncertainty index is recalculated to generate an updated uncertainty thermofield; Generate a new version of the database. For multiple rounds of sampling, repeat the above steps to gradually evolve the database version until the uncertainty index converges or the predetermined sampling completion rate is reached, forming the final version of the high-fidelity sampling database.

Images