Three-dimensional digital core reconstruction method for vesicular basalt

Abstract

The present application relates to the field of digital core modeling, and particularly relates to a three-dimensional digital core reconstruction method for vesicular basalt, comprising the following steps: extracting target characteristic parameters from CT scanning gray body data of a vesicular basalt sample, wherein the target characteristic parameters include porosity, cluster number, cluster size statistics, spatial uniformity index and simplified compactness; generating an initial three-dimensional digital core model according to the target characteristic parameters; iteratively optimizing the initial three-dimensional digital core model by using an adaptive Markov chain-Monte Carlo algorithm, so that the cluster characteristics of the optimized three-dimensional digital core model approximate the target characteristic parameters; and performing curvature smoothing post-processing on the optimized three-dimensional digital core model to obtain a three-dimensional digital core model for the vesicular basalt. The model not only highly matches the real vesicular basalt in macroscopic statistical characteristics, but also presents natural smooth curved surface characteristics in the micro-pore boundary morphology, thereby providing a reliable digital basis for subsequent rock physical property analysis.

Description

Three-dimensional digital core reconstruction method for vesicular basalt Technical Field

[0001] This invention relates to the field of digital core modeling technology, and in particular to a method for reconstructing three-dimensional digital cores of vesicular basalt. Background Technology

[0002] Basalt is the most widely distributed volcanic rock in the Earth's crust, and its prevalent vesicles give it excellent reservoir properties. In the deep layers of oil and gas basins, well-porosity basalt can not only serve as reservoirs for geothermal and oil and gas resources, but has also been proven to have enormous potential for carbon dioxide geological sequestration. Vesicle structure is the foundation of basalt's good reservoir and permeability capacity; the pore structure characteristics reflected by the morphology, size, number, and connectivity of vesicles within basalt are important evaluation indicators for basalt reservoirs.

[0003] Digital core technology acquires high-resolution images of rock samples, processes, segments, and reconstructs these images to create a three-dimensional digital core model with realistic pore structure. Subsequently, numerical simulations of various physical properties, such as seepage, electrical conductivity, elasticity, and thermal conductivity, are performed on this model to study the rock's reservoir properties and physical characteristics.

[0004] In recent years, three-dimensional digital core modeling has become a key research tool in rock physics and oil and gas exploration. Constructing digital core models of basalt can more intuitively and accurately depict the microscopic pore structure characteristics of basalt reservoirs and provide rich digital samples for studying the seepage and hydrocarbon storage characteristics of basalt reservoirs with different pore characteristics. Existing technologies mainly include physical experimental methods and numerical simulation methods. Physical experimental methods use high-precision experimental equipment to scan and image real cores to establish three-dimensional digital cores, such as X-ray CT scanning, laser scanning confocal microscopy, and focused ion beam scanning electron microscopy. Numerical reconstruction methods are based on two-dimensional images or statistical information and apply mathematical methods to reconstruct digital models equivalent to real cores, including Gaussian field methods, simulated annealing methods, multi-point geostatistical methods, Markov chain-Monte Carlo methods, process methods, deep learning methods, and hybrid methods.

[0005] However, existing technologies have significant shortcomings: physical experimental modeling is costly, has a limited number of samples, and the pore structure parameters are uncontrollable; existing numerical modeling methods have not been optimized for the spatial distribution of basalt pores, and their simulation of the connectivity of complex pore structures in basalt is poor. These problems limit the widespread application of basalt digital core models and make it difficult to meet the needs of basalt reservoir evaluation and carbon dioxide geological storage research. Summary of the Invention

[0006] In view of this, the present invention aims to provide a three-dimensional digital core reconstruction method for vesicular basalt. By optimizing and constraining the Markov chain-Monte Carlo modeling process, it achieves high-precision restoration of the spatial distribution law of basalt vesicles, improves the model's simulation effect on complex pore structures, and provides rich digital samples with controllable pore structure parameters for seepage simulation, pore structure evaluation, etc.

[0007] To achieve the above objectives, the technical solution created by this invention is implemented as follows:

[0008] A method for three-dimensional digital core reconstruction of vesicular basalt includes the following steps:

[0009] S1: Extract target feature parameters from the CT scan grayscale data of the vesicular basalt sample. The target feature parameters include porosity, number of clusters, cluster size statistics, spatial uniformity index and simplified compactness.

[0010] S2: Generate an initial three-dimensional digital core model based on the target feature parameters;

[0011] S3: The initial three-dimensional digital core model is iteratively optimized using the adaptive Markov chain-Monte Carlo algorithm, so that the cluster features of the optimized three-dimensional digital core model approximate the target feature parameters.

[0012] S4: Perform curvature smoothing post-processing on the optimized three-dimensional digital core model to obtain a three-dimensional digital core model of vesicular basalt.

[0013] Furthermore, step S1 specifically includes the following steps:

[0014] S11: Use a custom threshold to segment the pore space of the CT scan grayscale data to obtain a binary pore model;

[0015] S12: Perform connected component analysis on the binary pore model, identify connected pore clusters, calculate the number of voxels for each connected pore cluster as the cluster size, and collect cluster size statistics, including the maximum cluster volume, minimum cluster volume, and average cluster volume.

[0016] S13: Calculate the pairwise distances between the centroids of the connected pore clusters to obtain the distance set, and use the ratio of the standard deviation of the distance set to the mean of the distance set as the spatial uniformity index.

[0017] S14: For each connected pore cluster, calculate the axial boundary dimensions to obtain the bounding box volume, and use the ratio of the cluster volume prime number to the bounding box volume as the simplified compactness.

[0018] S15: The proportion of pore voxels in the binary pore model is used as the porosity to obtain the target feature parameters.

[0019] Furthermore, step S11 specifically includes the following steps:

[0020] S111: Read the grayscale data from the CT scan, and determine the grayscale dynamic range by statistically analyzing the minimum and maximum grayscale values; at the same time, use the interactive threshold module of Avizo software to select the porosity phase and obtain a fixed segmentation threshold.

[0021] S112: Binarize the grayscale volume data of CT scans with a fixed segmentation threshold to obtain preliminary pore segmentation results: voxels with grayscale values ​​less than the fixed segmentation threshold are identified as pore voxels, and voxels with grayscale values ​​greater than the fixed segmentation threshold are identified as matrix voxels, thus obtaining a binary pore model.

[0022] S113: Calculate the proportion of pore voxels in the binary pore model and use it as a reference porosity for the segmentation result. It is only used to indicate the porosity level of the current threshold segmentation to the user.

[0023] S114: Receives the target porosity input and uses it as the target parameter;

[0024] S115: If the pore morphology, connectivity, and noise level of the current binary pore model meet the requirements, solidify and output the binary pore model; if the pore morphology, connectivity, and noise level of the current binary pore model do not meet the requirements, use the interactive threshold module of Avizo software to recalibrate the fixed segmentation threshold, and repeat steps S112 to S114 until a binary pore model that meets the requirements is obtained.

[0025] Furthermore, step S2 specifically includes the following steps:

[0026] S21: Arrange seed points in three-dimensional space according to the preset minimum spacing;

[0027] S22: Centered on the seed point, derive the Gaussian kernel standard deviation parameter based on the target cluster size to generate a Gaussian field. By performing threshold binarization segmentation on the Gaussian field, a preliminary model of the natural morphology cluster is obtained.

[0028] S23: Apply spherical opening and closing operations to the prototype model of natural morphological clusters to smooth the pore boundaries;

[0029] S24: Calculate the difference between the porosity of the prototype model after step S23 and the target porosity, and select candidate voxels at the pore boundary to perform growth or erosion to adjust the number of pore voxels to the target porosity.

[0030] S25: Adjust the axial scaling factor of the Gaussian kernel to generate an anisotropic Gaussian field. Through threshold cutting operation and boundary fine-tuning operation based on statistical quantile, an initial three-dimensional digital core model that conforms to the target characteristic parameters is obtained.

[0031] Furthermore, step S25 specifically includes the following steps:

[0032] S251: Determine the scaling factor for each axis based on the anisotropy factor;

[0033] S252: Construct an anisotropic Gaussian kernel function centered on the seed point. The kernel function value is the negative exponent of the sum of the squares of each coordinate in the exponential decay term divided by twice the square of the corresponding axis scaling factor.

[0034] S253: In the frequency domain, the spectrum of the anisotropic Gaussian kernel is multiplied by the spectrum of the white noise field, and then an inverse transformation is performed to generate a continuous Gaussian random field.

[0035] S254: Calculate the cumulative distribution function of all voxel values ​​in the Gaussian random field, use the target porosity as the cut threshold to find the quantile, and set the voxels above the cut threshold as pores to obtain the initial three-dimensional digital core model.

[0036] Furthermore, step S3 specifically includes the following steps:

[0037] S31: Calculate the deviation between the cluster size and the target feature parameters of the current 3D digital core model, and dynamically adjust the weights of splitting and connectivity operations;

[0038] S32: If a super-large cluster is detected, a splitting operation is triggered. A cutting plane is randomly selected within the super-large cluster to divide the super-large cluster into two clusters.

[0039] S33: If adjacent clusters are detected, the matrix voxels are converted into porosiform voxels along the shortest path between the adjacent clusters to build a connected path and merge the adjacent clusters.

[0040] S34: Generate a new 3D digital core model state using an adaptive Markov chain-Monte Carlo algorithm, and calculate the energy difference between the new 3D digital core model state and the current 3D digital core model state using an energy function;

[0041] S35: If the energy difference is less than zero, accept the new three-dimensional digital core model state; if the energy difference is greater than zero, accept the new three-dimensional digital core model state according to the probability acceptance criterion; through iterative iteration, make the cluster features of the optimized three-dimensional digital core model approximate the target feature parameters.

[0042] Furthermore, step S34 specifically includes the following steps:

[0043] S341: Calculate the deviation between the porosity of the new three-dimensional digital core model and the target porosity, and use it as the porosity energy component;

[0044] S342: Calculate the difference between the number of clusters in the new 3D digital core model state and the number of clusters in the target feature parameters, as well as the difference between the cluster size statistics in the new 3D digital core model state and the cluster size statistics in the target feature parameters, as the cluster distribution energy component.

[0045] S343: Calculate the difference between the bounding box volume ratio and the simplified compactness of each connected pore cluster in the state of the new 3D digital core model as a morphological energy component;

[0046] S344: Calculate the difference between the cluster centroid distance variation coefficient and the spatial homogeneity index of the target characteristic parameters in the state of the new three-dimensional digital core model, as the spatial matching energy component;

[0047] S345: The energy components of porosity, cluster distribution, morphology and spatial matching are weighted and summed to obtain the total energy of the new three-dimensional digital core model state. The difference between the total energy of the new three-dimensional digital core model state and the total energy of the current three-dimensional digital core model state is calculated as the energy difference.

[0048] Furthermore, step S4 specifically includes the following steps:

[0049] S41: Extract the pore boundary voxels of the optimized 3D digital core model, calculate the difference between the local neighborhood average and the center value of each boundary voxel, and use it as the curvature symbol;

[0050] S42: If the curvature sign is negative, the corresponding matrix voxel is converted into a pore voxel for growth; if the curvature sign is positive, the corresponding pore voxel is converted into a matrix voxel for erosion.

[0051] S43: Apply sphere opening and closing operations to the three-dimensional digital core model processed in step S42 to smooth the pore boundaries;

[0052] S44: Convert the three-dimensional digital core model processed in step S43 into a floating-point model, perform Gaussian blur operation on the floating-point model, and then perform threshold binarization segmentation to obtain a smooth boundary model.

[0053] S45: If there are residual super-large clusters, calculate the distance field for the super-large clusters, identify the local maxima of the distance field as seed points, and obtain the three-dimensional digital core model of vesicular basalt by segmenting the super-large clusters using the watershed algorithm.

[0054] Furthermore, step S45 specifically includes the following steps:

[0055] S451: Calculate the shortest distance from each porosity voxel within the supercluster to the cluster boundary, and use it as the distance field;

[0056] S452: Identify local maxima in the distance field as seed markers, with each seed marker corresponding to a sub-cluster center;

[0057] S453: Invert the distance field into an elevation field and diffuse it from various sub-markers to the neighborhood simultaneously. When the diffusion areas of different seed marks meet, a dividing boundary is established at the meeting point.

[0058] S454: Divide the super-large cluster into different sub-clusters according to the dividing boundary, update the model cluster list, and obtain the three-dimensional digital core model of vesicular basalt.

[0059] Compared with the prior art, the present invention can achieve the following technical effects:

[0060] This invention first extracts target porosity feature parameters precisely from CT scan grayscale data using user-defined thresholds and connected component analysis. Then, based on seed points and Gaussian kernels, it generates an initial digital core model with a natural morphology and dynamically adjusts it to the target porosity. Next, it introduces an adaptive Markov chain-Monte Carlo algorithm to dynamically adjust the operation weights of cluster splitting and connected path construction based on the deviation between the current 3D digital core model and the target features. Iterative optimization is achieved using a multi-component energy function combined with a probability acceptance criterion. Finally, post-processing techniques such as curvature boundary smoothing and watershed segmentation are used to eliminate residual defects. Ultimately, this method achieves a precise match between the number, size distribution, and spatial uniformity of clusters and the target features. This high-fidelity matching of multiple features significantly improves the realism and representativeness of the digital core model, providing a more accurate digital foundation for rock physics simulation and reservoir evaluation. Attached Figure Description

[0061] Figure 1 is a schematic flowchart of the three-dimensional digital core reconstruction method for vesicular basalt according to an embodiment of the present invention. Detailed Implementation

[0062] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not constitute a limitation thereof.

[0063] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0064] The invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0065] As shown in Figure 1, an embodiment of the present invention provides a method for reconstructing a three-dimensional digital core of porous basalt, comprising the following steps:

[0066] S1: Extract target feature parameters from the CT scan grayscale data of the vesicular basalt sample. The target feature parameters include porosity, number of clusters, cluster size statistics, spatial uniformity index, and simplified compactness.

[0067] In one embodiment, the process of extracting target feature parameters from CT scan grayscale data of a vesicular basalt sample includes the following steps:

[0068] S11: Use a custom threshold to segment the pore space of the CT scan grayscale data to obtain a binary pore model.

[0069] Step S11 specifically includes the following steps:

[0070] S111: Read the grayscale data from the CT scan, and determine the grayscale dynamic range by statistically analyzing the minimum and maximum grayscale values; at the same time, use the interactive threshold module of Avizo software to select the porous phase and obtain a fixed segmentation threshold.

[0071] S112: Binarize the grayscale volume data of CT scans with a fixed segmentation threshold to obtain preliminary pore segmentation results: voxels with grayscale values ​​less than the fixed segmentation threshold are identified as pore voxels, and voxels with grayscale values ​​greater than the fixed segmentation threshold are identified as matrix voxels, thus obtaining a binary pore model.

[0072] S113: Calculate the proportion of pore voxels in the binary pore model and use it as a reference porosity for the segmentation result. This reference porosity is only used to indicate the porosity level of the current threshold segmentation and is not used as a hard constraint target for subsequent optimization. At the same time, save the binary model for use in subsequent processes.

[0073] S114: Receives the input target porosity and uses it as the target parameter for subsequent initial three-dimensional digital core model generation, structural optimization, and porosity hard constraint correction.

[0074] S115: If the pore morphology, connectivity, and noise level of the current binary pore model meet the requirements, solidify and output the binary pore model as input for subsequent digital core modeling and optimization post-processing; if the pore morphology, connectivity, and noise level of the current binary pore model do not meet the requirements, use the interactive threshold module of Avizo software to recalibrate the fixed segmentation threshold, and repeat steps S112 to S114 until a binary pore model that meets the requirements is obtained.

[0075] Specifically, the process begins by reading the grayscale data from the CT scan and calculating its minimum and maximum grayscale values ​​to characterize the dynamic range of grayscale under the current scanning conditions. Based on this, the basic settings for grayscale scale unification and visualization window width / level are established. Subsequently, the CT scan grayscale data is imported into Avizo software. In the Interactive Threshold module, combining the grayscale histogram distribution and multi-section slice / volume rendering visualization results, the grayscale intervals of the porosity phase are manually interpreted, selected, and interactively verified to obtain the global fixed segmentation threshold corresponding to the pore-matrix phase separation (i.e., the pore threshold parameter obtained by Avizo interactive threshold calibration). This fixed segmentation threshold is then written into the segmentation program as the sole criterion for global threshold segmentation. Following this, a one-time global threshold binarization is performed on the CT scan grayscale data using the fixed segmentation threshold: for any voxel, voxels with grayscale values ​​less than or equal to the fixed segmentation threshold are marked as pore voxels, and voxels with grayscale values ​​greater than the fixed segmentation threshold are marked as matrix voxels, thus obtaining a preliminary binary pore model. Furthermore, to suppress noise and isolated discrete voxels introduced by threshold segmentation and improve the topological consistency of pore boundaries, a three-dimensional morphological opening operation is performed on the preliminary binary results. Small cubic structuring elements are used to implement an erosion-dilation sequence to remove disconnected small-scale pseudo-pores and weaken discrete noise points, resulting in a cleaned binary pore model. Based on this, the number of pore voxels and the overall prime number are statistically analyzed, and the porosity is calculated for quantitative characterization and consistency checks of the segmentation results. Since the threshold is determined by interactive Avizo calibration and used as a fixed parameter in the program, the above process does not involve threshold iterative updates and convergence searches constrained by the target porosity. The final morphologically cleaned binary pore model is then output as the pore space segmentation result.

[0076] For example, in processing grayscale data from basalt CT scans with a resolution of approximately 3 μm, the grayscale range is 0–65535. First, in the interactive thresholding module of Avizo software, combined with grayscale histogram peak-valley characteristics and multi-slice comparison, the porous phase is manually selected and verified to determine a globally fixed segmentation threshold. This threshold is then used as the input parameter for one-time threshold binarization. Three-dimensional morphological opening operations are then used for curvature smoothing post-processing to denoise the data, ultimately outputting a binary pore model. Simultaneously, porosity is calculated and output as a statistical indicator to verify the consistency between the segmentation results and the overall pore structure characteristics of the sample. This method, based on manually calibrated thresholds, avoids threshold drift and unstable convergence that may be introduced by automatic threshold iteration strategies under conditions of noise, beam hardening effects, or local grayscale drift-induced grayscale non-uniformity. This improves the robustness and repeatability of subsequent feature extraction, such as pore connectivity analysis, pore size distribution statistics, and pore cluster identification.

[0077] In another implementation, to further enhance the suppression of threshold jaggedness and microcavity artifacts, a small-scale closing operation (dilation-erosion sequence) can be added after the morphological opening operation to fill the small pores formed during binarization and smooth the pore boundaries, while maintaining the stability of the pore voxel statistics as much as possible. This "open-close" joint morphological cleaning strategy is particularly suitable for basalt samples with uneven pore spatial distribution, significant gray-scale fluctuations, or low local contrast, and can significantly reduce the interference of pseudo-isolated pores, boundary gaps, and noisy connectivity bridges on the subsequent topological connectivity network and cluster identification results.

[0078] S12: Perform connected component analysis on the binary pore model, identify connected pore clusters, calculate the number of voxels for each connected pore cluster as the cluster size, and collect cluster size statistics, including the maximum cluster volume, minimum cluster volume, and average cluster volume.

[0079] Specifically, the 26-adjacency criterion is used to traverse the pore voxels in the binary pore model, and all unconnected pore regions are marked as independent clusters. For each identified cluster, the number of voxels it contains is counted as the cluster size. After traversing all clusters, the maximum value, minimum value, and arithmetic mean of all cluster sizes are recorded as the cluster size statistics.

[0080] For example, in a typical binary porosity model of vesicular basalt, connected domain analysis may identify hundreds of clusters, ranging in size from tens of voxels to tens of thousands of voxels, with the average cluster size reflecting the main vesicular scale. This statistic directly provides a size distribution benchmark for subsequent digital core reconstruction.

[0081] S13: Calculate the pairwise distances between the centroids of the connected pore clusters to obtain the distance set, and use the ratio of the standard deviation of the distance set to the mean of the distance set as the spatial uniformity index.

[0082] Specifically, for each connected pore cluster, the centroid coordinates are calculated, which is the average of the coordinates of all voxels within the cluster. All cluster centroids are collected, and the Euclidean distance between any two centroids is calculated, forming a distance set. If the number of clusters is greater than one, the standard deviation and mean of this set are calculated; the ratio of these two values ​​is the spatial uniformity index. If there is only one cluster, the index is defined as 0.

[0083] For example, in basalt samples with relatively uniform vesicle distribution, the standard deviation of the centroid distance set is small, the average value is large, and the ratio is usually below 0.2, indicating high uniformity in cluster spacing. However, in samples with locally clustered vesicles, this ratio may exceed 0.5, reflecting uneven distribution. The introduction of this index helps to quantify the spatial arrangement of vesicles in basalt reservoirs, provides a quantitative basis for imposing uniformity constraints on reconstruction models, and effectively reduces the spatial distribution error of traditional methods.

[0084] In one embodiment, when the number of clusters is large, the centroid set can be sampled first to reduce the computational load while maintaining the statistical representativeness of the indicators. This approach can significantly shorten the computation time when processing large-size digital cores.

[0085] S14: For each connected pore cluster, calculate the axial boundary dimensions to obtain the bounding box volume, and use the ratio of the cluster volume prime number to the bounding box volume as the simplified compactness.

[0086] Specifically, for each cluster, its minimum and maximum coordinates in the X, Y, and Z directions are determined, and the bounding box volume is calculated by multiplying the axial ranges. The simplified compactness of the cluster is obtained by dividing the number of voxels by the bounding box volume. The average or distribution of the calculated compactness for all clusters is then used as the overall compactness feature.

[0087] For example, in basalts dominated by near-spherical vesicles, the compactness is typically above 0.6, while in samples with more elongated or irregular vesicles, this value may drop below 0.4. This index directly reflects the compactness of the vesicle shape and is of great significance for assessing reservoir performance. By incorporating it into the reconstruction objective, the pore morphology of the generated model can be made closer to that of the real sample.

[0088] S15: The proportion of pore voxels in the binary pore model is used as the porosity to obtain the target feature parameters.

[0089] Specifically, the porosity is obtained by dividing the total number of voxels labeled as pores in the binary porosity model by the total number of primes in the model. Porosity, cluster size statistics, spatial uniformity index, simplified compactness, and the number of related clusters are all included as the target feature parameter set.

[0090] For example, after the complete extraction process, a typical basalt sample may have a porosity of 0.25, a cluster count of approximately 350, an average cluster size of 5000 voxels, a spatial homogeneity index of 0.28, and an average compactness of 0.52. These parameters provide target feature benchmarks for the subsequent construction of a three-dimensional digital core model, ensuring that the generated digital core faithfully reproduces the real vesicular structure in multiple dimensions, and significantly improving the accuracy of numerical simulations of basalt reservoir seepage and storage performance.

[0091] S2: Generate an initial three-dimensional digital core model based on the target feature parameters.

[0092] In one embodiment, the process of generating an initial three-dimensional digital core model based on target feature parameters includes the following steps:

[0093] S21: Arrange seed points in three-dimensional space according to the preset minimum spacing.

[0094] Specifically, the average pore diameter is obtained from the target feature parameters. This value is then used as the minimum interval to generate a series of seed point coordinates uniformly or randomly within the 3D modeling space. These seed points will serve as the center positions for subsequently generated pore clusters. The spacing of these seed points ensures that the pore clusters in the initial 3D digital core model do not overlap excessively, providing a basis for simulating the spatial distribution of pores in real basalt.

[0095] S22: Using the seed point as the center, derive the Gaussian kernel standard deviation parameter based on the target cluster size to generate a Gaussian field. By performing threshold binarization segmentation on the Gaussian field, a preliminary model of the natural morphology cluster is obtained.

[0096] Step S22 specifically includes the following steps:

[0097] Step S221: Determine the Gaussian kernel standard deviation parameter.

[0098] The average volume or equivalent diameter of the target cluster is extracted from the target feature parameters. Assuming the target cluster is approximately spherical, its equivalent radius R can be calculated using the volume formula. The standard deviation σ of the Gaussian kernel is directly related to the cluster size. In one embodiment, σ is set to R divided by an empirical coefficient, typically ranging from 2 to 3, for example, σ = R / 2.5. This setting ensures that the attenuation range of the Gaussian function in three-dimensional space matches the size of the target cluster.

[0099] Step S222: Generate a Gaussian field.

[0100] For each seed point, using its three-dimensional coordinates as the center and the standard deviation σ determined in step S221, the value of the three-dimensional Gaussian function is calculated for each voxel position in the neighborhood surrounding that point. The form of the Gaussian function is: By superimposing the Gaussian fields generated by all seed points, a continuous three-dimensional scalar field with values ​​between 0 and 1 is obtained. Regions with higher values ​​in the field correspond to the spaces where future pore clusters may appear.

[0101] Step S223: Threshold binarization yields natural morphology clusters.

[0102] A global threshold T is set for the superimposed continuous Gaussian field. Voxels with values ​​greater than T are marked as pores, and voxels with values ​​less than or equal to T are marked as matrix. The selection of the threshold T is crucial, as it initially determines the total volume of generated pores. In one embodiment, the initial threshold T can be set to 0.5, and then dynamically adjusted based on the difference between the porosity of the initial binarization result and the target porosity. Through threshold cutting, the originally smoothly transitioned Gaussian field is transformed into a set of binary pore clusters with continuous, soft edges and natural morphology. These clusters are approximately spherical or ellipsoidal in shape, avoiding artificially generated sharp edges.

[0103] S23: Apply spherical opening and closing operations to the prototype model of natural morphological clusters to smooth the pore boundaries.

[0104] For the prototype model of the natural morphological clusters obtained in step S22, morphological operations are performed using a spherical structuring element with radius r. First, an opening operation is performed, i.e., erosion followed by dilation. This operation eliminates small spikes and isolated noise points on the pore surface. Then, a closing operation is performed, i.e., dilation followed by erosion. This operation fills the tiny pits and holes inside the pores. The spherical structuring element ensures consistent smoothness in all directions. After the spherical opening and closing operations, the pore boundaries of the model become smoother and more continuous, removing any artificial traces that may have been introduced by discretization and initial generation.

[0105] S24: Calculate the difference between the porosity of the prototype model after step S23 and the target porosity, and select candidate voxels at the pore boundary to perform growth or erosion to adjust the number of pore voxels to the target porosity.

[0106] Step S24 specifically includes the following steps:

[0107] Step S241: Calculate the porosity difference.

[0108] The total number of voxels labeled as pores in the current prototype model is denoted as N. current Based on the target porosity and the total volume of the model, calculate the target pore volume number N. target Calculate the difference ΔN = N target -N current .

[0109] Step S242: Identify boundary candidate voxels.

[0110] The prototype model is traversed to identify all voxels located at the interface between pores and matrix. Specifically, they are divided into two categories: one category is currently matrix, but has at least one directly adjacent voxel that is a pore, and these voxels are considered as growth candidate voxels; the other category is currently pore, but has at least one directly adjacent voxel that is matrix, and these voxels are considered as erosion candidate voxels.

[0111] Step S243: Perform growth or erosion operations.

[0112] If ΔN > 0, it indicates insufficient porosity voxels, requiring a growth operation. From the candidate voxels for growth, prioritize points highly surrounded by porosity voxels, invert their state from matrix to pore, and repeat this process until the actual increase in voxels approaches ΔN. If ΔN < 0, it indicates an excess of porosity voxels, requiring an erosion operation. From the candidate voxels for erosion, prioritize points located at the pore edge protruding position and highly surrounded by matrix, invert their state from pore to matrix, and repeat this process until the actual decrease in voxels approaches |ΔN|. Through this local adjustment at the boundary, the model porosity can be accurately adjusted to the target value while minimizing changes to the overall pore structure.

[0113] S25: Adjust the axial scaling factor of the Gaussian kernel to generate an anisotropic Gaussian field. Through threshold cutting operation and boundary fine-tuning operation based on statistical quantile, an initial three-dimensional digital core model that conforms to the target characteristic parameters is obtained.

[0114] In step S25, the process of adjusting the axial scaling factor of the Gaussian kernel to generate an anisotropic Gaussian field, and obtaining an initial three-dimensional digital core model that conforms to the target characteristic parameters through threshold cutting operation based on statistical quantiles and boundary fine-tuning operation, specifically includes the following steps:

[0115] S251: Determine the scaling factor for each axis based on the anisotropy factor.

[0116] The principal axis scaling factor is 1 plus an anisotropy factor, and the secondary axis scaling factor is 1 plus a portion of the anisotropy factor. Specifically, the anisotropy factor quantifies the stretching of basalt pores in a certain direction, typically ranging from 0 to 1, where 0 represents complete isotropy and a larger value indicates more pronounced stretching. Setting the principal axis scaling factor to 1 plus the anisotropy factor significantly extends the relevant length in the principal direction. For example, when the anisotropy factor is 0.5, the principal axis scaling factor is 1.5, ensuring slower Gaussian decay in that direction and resulting in elongated pore morphology. The secondary axis scaling factor uses a portion of the anisotropy factor, for example, 0.3 times, making the secondary axis scaling factor 1 plus 0.3 times the anisotropy factor. When the anisotropy factor is 0.5, the secondary axis scaling factor is 1.15. This asymmetric scaling method avoids excessive deformation of the secondary axis, maintaining the overall natural synergy of the pores. In one possible implementation, the main axis direction is pre-determined as the main extension direction of pores based on principal component analysis of the actual core, thereby making the generated initial morphology more consistent with the directional characteristics of the original pores in basalt.

[0117] S252: Construct an anisotropic Gaussian kernel function centered on the seed point. The kernel function value is the negative exponent of the sum of the squares of each coordinate in the exponential decay term divided by twice the square of the corresponding axis scaling factor.

[0118] For example, the seed point is taken as the center of the stomatal cluster, and an anisotropic Gaussian kernel function is constructed in three-dimensional space with this point as the origin. The expression for the anisotropic Gaussian kernel function is: ,in, Three-dimensional spatial coordinates The kernel function value at that location, Based on the relevant length, derived from the target average pore diameter, , , The scaling factors for each axis are determined in step S251. This construction method causes the anisotropic Gaussian kernel function to decay slowly along the principal axis, forming an ellipsoidal distribution, while decaying relatively quickly along the secondary axis, resulting in an irregular but natural stretched shape. It should be noted that the denominator in the exponent term is twice the square of the corresponding axis scaling factor, ensuring that the standard deviation of the Gaussian distribution is proportional to the scaling factor, thereby precisely controlling the aspect ratio of the pores. In one embodiment, when the target pores exhibit vertical layering stretching, the Z-axis is set as the principal axis, and the scaling factor is increased accordingly. This allows the generated clusters of pores to extend slenderly along the vertical direction, simulating the pore morphology in deep porous basalt affected by gravity, thereby improving the initial model's preliminary approximation of the real structure.

[0119] S253: In the frequency domain, the spectrum of the anisotropic Gaussian kernel is multiplied by the spectrum of the white noise field, and then an inverse transformation is performed to generate a continuous Gaussian random field.

[0120] Specifically, firstly, a three-dimensional Fourier transform is performed on the anisotropic Gaussian kernel to obtain its spectrum. Then, a three-dimensional white noise field with the same size as the model mesh is generated and subjected to a Fourier transform. The point product operation is performed in the frequency domain, that is, the Gaussian kernel spectrum is directly multiplied by the white noise spectrum as a filter to achieve an equivalent spatial domain convolution. Then, an inverse Fourier transform is performed to obtain a continuous Gaussian random field. This Gaussian random field presents multiple randomly distributed clouds in space, and the shape of each cloud is controlled by the Gaussian kernel function and has a specified anisotropy. This frequency domain implementation significantly improves computational efficiency, especially with large-size three-dimensional meshes, avoiding the huge overhead of point-by-point convolution in the spatial domain. At the same time, the introduction of white noise ensures that each generation has randomness, so that the morphology of different vesicle clusters varies slightly, which is closer to the natural variability of real basalt vesicles. In one possible implementation, candidate fields with different random seeds are generated in parallel multiple times, and the field whose spatial statistical characteristics best match the target is selected as the subsequent basis, thereby further improving the diversity and fidelity of the initial morphology.

[0121] S254: Calculate the cumulative distribution function of all voxel values ​​in the Gaussian random field, use the target porosity as the cut threshold to find the quantile, and set the voxels above the cut threshold as pores to obtain the initial three-dimensional digital core model.

[0122] For example, firstly, all voxel values ​​in the Gaussian random field are sorted, or a cumulative distribution function is constructed using a histogram method. Then, the corresponding quantile is found based on the target porosity, ensuring that the proportion of voxels above this threshold is precisely equal to the target value. This statistical thresholding ensures that the porosity of the initial model is strictly matched, avoiding volume deviations caused by traditional fixed segmentation thresholds. After segmentation, high-value regions form connected or isolated pore clusters, exhibiting natural and smooth boundaries. In one embodiment, if the target porosity is 0.18, the value corresponding to the 0.82 quantile in the cumulative distribution function is used as the threshold, directly setting voxels above this value to 1 and the rest to 0, thus obtaining an initial three-dimensional digital core model with accurate volume. This approach, combined with the aforementioned anisotropic field generation, ensures that the initial morphology of the initial three-dimensional digital core model simultaneously meets the requirements for porosity and shape stretching, significantly reducing the number of subsequent optimization iterations and improving overall reconstruction efficiency.

[0123] S3: The initial 3D digital core model is iteratively optimized using the adaptive Markov chain-Monte Carlo algorithm, so that the cluster features of the optimized 3D digital core model approximate the target feature parameters.

[0124] In one embodiment, the process of iteratively optimizing the initial three-dimensional digital core model using an adaptive Markov chain-Monte Carlo algorithm to make the cluster features of the optimized three-dimensional digital core model approximate the target feature parameters specifically includes the following steps:

[0125] S31: Calculate the deviation between the cluster size and the target feature parameters of the current 3D digital core model, and dynamically adjust the weights of splitting and connectivity operations.

[0126] In one embodiment, the process of calculating the deviation between the cluster size of the current 3D digital core model and the target feature parameters, and dynamically adjusting the weights of splitting and connectivity operations, specifically includes the following steps:

[0127] S311: Extract all pore clusters in the current 3D digital core model.

[0128] Connectivity domains are marked on the current three-dimensional digital core model to identify all independent regions that are interconnected and composed of pore voxels. Each region is a pore cluster.

[0129] S312: Calculate the cluster size distribution characteristics of the current three-dimensional digital core model.

[0130] The number of voxels contained in all pore clusters is counted to obtain the cluster size set of the current three-dimensional digital core model, and then the average value, maximum value and size distribution histogram of the set are calculated.

[0131] S313: Calculate cluster size deviation.

[0132] The current cluster size distribution characteristics obtained in step S312 are compared with the target feature parameters extracted from the real vesicular basalt sample. The deviation calculation includes multiple dimensions, such as calculating the absolute difference between the average cluster volume of the current 3D digital core model and the average cluster volume of the target feature parameter, calculating the ratio of the maximum cluster volume of the current 3D digital core model to the maximum cluster volume of the target feature parameter, and calculating the L1 norm distance between the current cluster size histogram and the target histogram. These deviation values ​​collectively quantify the difference in cluster size between the current 3D digital core model and the real vesicular basalt sample.

[0133] S314: Dynamically adjust operation weights based on deviation.

[0134] Set the baseline selection probabilities for splitting and connecting operations. After each iteration or every few iterations, temporarily adjust the actual selection probabilities of these two operations based on the deviation calculated in step S313. Specifically, if a very large cluster with a volume significantly exceeding the maximum cluster volume of the target feature parameter is detected in the current 3D digital core model, and the number of clusters is less than the target value, then the cluster size deviation is mainly due to the clusters being too large and too few. In this case, the weight of the splitting operation is increased, raising its probability of being selected in the next iteration to a high value, such as 70%, while correspondingly reducing the weight of other operations (such as simple voxel flipping). Conversely, if a large number of small clusters with a volume much smaller than the average cluster volume of the target feature parameter are detected in the current 3D digital core model, and the number of clusters far exceeds the target value, then the deviation is mainly due to the clusters being too fragmented and too numerous. In this case, the weight of the connecting operation is increased, for example, to 60%, to promote the merging of small clusters.

[0135] S32: If a super-large cluster is detected, a split operation is triggered. A cutting plane is randomly selected within the super-large cluster to divide the super-large cluster into two clusters.

[0136] In one embodiment, step S32 specifically includes the following steps:

[0137] S321: Ultra-large cluster detection.

[0138] Traverse all pore clusters in the current model and mark clusters with a volume greater than a preset threshold (which is usually set to 1.2 to 1.5 times the target maximum cluster volume) as super-large clusters to be split.

[0139] S322: Randomly select a cutting reference point.

[0140] Within the marked super-large cluster, a pore voxel is randomly selected as the point through which the cutting plane passes. To ensure the effectiveness of the cutting, this random point must be at least a certain distance from the cluster boundary to avoid producing tiny fragments at the edge.

[0141] S323: Determine the direction of the cutting plane.

[0142] A random 3D unit vector is generated as the normal vector of the cutting plane. The direction of this normal vector is sampled uniformly with equal probability in 3D space to ensure the isotropy of the cutting direction in long-term iterations and avoid introducing artificial direction bias.

[0143] S324: Perform cluster splitting.

[0144] Using the points selected in step S322 and the normal vector determined in step S323, an infinitely extending virtual plane is constructed. All porosity voxels in this supercluster are classified according to their positional relationship with this plane: voxels located on the positive side of the plane's normal vector are assigned to sub-cluster A, and voxels located on the negative side are assigned to sub-cluster B. Subsequently, connectivity analysis is performed on sub-clusters A and B respectively, confirming that they are internally connected and separated by the plane, with no porosity voxels connected between them. Thus, the original supercluster has been successfully divided into two independent porosity clusters.

[0145] S325: Update model data.

[0146] The two newly generated clusters after segmentation are added to the model's pore cluster list, and the original super-large cluster is removed from the list. The porosity, number of clusters, and cluster size distribution of the model are recalculated.

[0147] S33: If an adjacent cluster is detected, the matrix voxel is converted into a pore voxel along the shortest path between the adjacent clusters to build a connected path and merge the adjacent clusters.

[0148] In one embodiment, step S33 specifically includes the following steps:

[0149] S331: Neighboring small cluster detection.

[0150] Calculate the minimum Euclidean distance between all pairs of pore clusters. This distance is defined as the shortest distance between any pair of boundary voxels in two clusters. If the minimum distance between two clusters is less than or equal to a preset adjacency threshold (e.g., the side length of two voxels), then the two clusters are determined to be spatially adjacent candidate connected pairs.

[0151] S332: Filter cluster pairs to be connected.

[0152] From all adjacent cluster pairs, select "small cluster pairs" where both cluster volumes are smaller than the target average cluster volume. Prioritize cluster pairs that are closest in distance and whose sum of volumes is still within the target cluster size range.

[0153] S333: Determine the shortest connected path.

[0154] For a selected pair of subclusters, find the pair of boundary voxels that are closest to each other in both clusters, denoted as voxel P (belonging to cluster A) and voxel Q (belonging to cluster B). Then, in a 3D voxel mesh, calculate the shortest path from voxel P to voxel Q. This path is obtained by searching in a space containing only matrix voxels using either Bressnum's line algorithm or a 3D Dijkstra algorithm, ensuring that the path consists of a series of adjacent matrix voxels.

[0155] S334: Construct a connected path.

[0156] The state of all matrix voxels on the shortest path obtained in step S333 is modified to that of porosity voxels. This operation is equivalent to "digging" a narrow channel in the matrix barrier that originally isolated the two small clusters.

[0157] S335: Verify and complete the merge.

[0158] Due to the establishment of the connecting path, the originally independent clusters A and B are connected into a single connected whole through the new pore channels. The newly formed unified region is marked as a connected domain, confirming it as a single pore cluster. The model cluster list is updated, adding the merged new cluster and deleting the original clusters A and B.

[0159] S34: Generate a new 3D digital core model state using an adaptive Markov chain-Monte Carlo algorithm, and calculate the energy difference between the new 3D digital core model state and the current 3D digital core model state using an energy function.

[0160] In one embodiment, the process of generating a new 3D digital core model state using an adaptive Markov chain-Monte Carlo algorithm specifically includes the following steps:

[0161] S341a: Generate candidate new model states.

[0162] Based on the dynamically adjusted operation weights, an operation (such as splitting, connecting, or flipping the state of a basic voxel) is randomly selected and executed on the current model, thereby generating a candidate new 3D digital core model.

[0163] S342a: Calculate the energy function value.

[0164] The energy function E is used to comprehensively evaluate the similarity between the 3D digital core model and the target feature parameters. Its specific form is a weighted sum of multiple feature deviations, for example: .in, This indicates the porosity of the new model. Porosity represents the target characteristic parameter. This indicates the differences in the distribution of cluster quantity and size. This represents the difference in spatial statistical characteristics such as the correlation function between two points. Calculate the energy value of the current model state for each point. Energy values ​​of candidate new model states .

[0165] S343a: Calculate the energy difference.

[0166] calculate . The size and sign of the candidate variation reflect the degree to which the three-dimensional digital core model deviates from or approaches the target feature parameter.

[0167] In one embodiment, the process of calculating the energy difference between the new state of the three-dimensional digital core model and the current state of the three-dimensional digital core model using an energy function specifically includes the following steps:

[0168] S341b: Calculate the deviation between the porosity of the new three-dimensional digital core model and the target porosity, and use it as the porosity energy component.

[0169] Specifically, voxel statistics are performed on the newly generated 3D digital core model to obtain the proportion of the current total number of pore voxels to the total number of voxels, which is taken as the porosity of the new 3D digital core model. This porosity is then subtracted from the porosity of the target feature parameters extracted from the original CT scan grayscale data, and the absolute value or the difference of squares is taken as the porosity energy component. This component directly reflects the volume ratio deviation and is usually expressed in square form to amplify larger errors.

[0170] For example, if the target porosity of a basalt sample is 0.28, and the new model calculates 0.26, the deviation is 0.02. The porosity energy component can be taken as a penalty term of 0.0004, so that the volume deviation can be corrected first in the subsequent acceptance criteria.

[0171] S342b: Calculate the difference between the number of clusters in the new 3D digital core model state and the number of clusters in the target feature parameters, as well as the difference between the cluster size statistics in the new 3D digital core model state and the cluster size statistics in the target feature parameters, as the cluster distribution energy component.

[0172] In one embodiment, a 26-adjacent connected component analysis is first performed on the new model to identify all independent pore clusters, count the total number of clusters, and calculate the voxel count for each cluster to obtain a cluster size set, including the maximum cluster size, minimum cluster size, and average cluster size. The absolute difference is obtained by subtracting the number of clusters in the new model from the number of clusters in the target feature parameter. The L1 norm difference or chi-square distance is calculated between the histogram of the new model's cluster size statistics and the target histogram. The two are then weighted and combined to form a cluster distribution energy component. This component imposes a high penalty for cases where there are too many or too few pore clusters or mismatched size distributions in basalt.

[0173] Specifically, when the target number of clusters is 45 and the average cluster size corresponds to a volume of 8500 voxels, and the number of clusters in the new model is 32 and there is an ultra-large cluster with a volume of more than 20000 voxels, the cluster number deviation contributes a large penalty. At the same time, the maximum cluster size deviation further increases this component value, prompting subsequent iterations to trigger splitting operations to balance the cluster distribution.

[0174] In another implementation, if the basalt sample shows a cluster size distribution that is biased towards small, dense clusters, the weight of the cluster number deviation can be increased to make the energy component more sensitive to the case of too many clusters, thereby quickly driving the connectivity path construction operation, reducing the proportion of isolated small clusters, and improving overall connectivity.

[0175] Preferably, for the few cases of large clusters dominating deep basalt reservoir samples, an additional term of the square difference between the maximum cluster size and the target maximum cluster size can be added to the energy component of the cluster distribution to further enhance the penalty effect on ultra-large clusters and avoid the appearance of overly connected regions in the model that do not conform to geological reality.

[0176] For example, in one embodiment, the target maximum cluster size is 15,000 voxels. The new model has a cluster of 28,000 voxels. This component is significantly increased by the squared difference, causing most harmful variations to be rejected until the split operation splits the cluster into two sub-clusters of approximately 12,000 voxels, which reduces the energy and thus achieves cluster distribution convergence toward the target.

[0177] S343b: Calculate the difference between the bounding box volume ratio and the simplified compactness of each connected pore cluster in the state of the new 3D digital core model as a morphological energy component.

[0178] In one embodiment, the axial boundary range of each cluster in the new model is calculated to obtain the bounding box volume. The compactness of a single cluster is then obtained by dividing the cluster voxels by the bounding box volume. The average value or distribution statistics of all clusters are taken, and the difference between this and the average compactness of the target feature parameter is calculated. The squared or absolute values ​​are then summed to obtain the morphological energy component. This component measures the compactness of the stomatal cluster shape, penalizing elongated or highly irregular cluster morphologies.

[0179] Specifically, a target simplified compactness of 0.62 indicates that the pores are relatively rounded. If the average compactness of the new model drops to 0.48, it indicates that the cluster shape is overstretched. This component increases the curvature smoothing or erosion operation in post-processing, making the cluster boundaries fuller.

[0180] In one possible implementation, to accommodate different basalt types, the compactness deviations of large and small clusters can be calculated separately and summed with weights, with the compactness deviation of small clusters having a higher weight, in order to prevent small pores from being stretched into unnatural tubular structures.

[0181] For example, in amygdaloidal basalt formed by rapid cooling after volcanic eruption, the target compactness is relatively high, close to 0.75. After the new model is initially generated, the compactness is 0.59. This component produces a significant penalty through differential accumulation. Subsequently, by adjusting the Gaussian kernel deformation parameters and performing sphere opening and closing operations, the compactness is gradually increased to 0.71, significantly improving the morphological fidelity.

[0182] Preferably, when the vesicular basalt sample shows that some clusters are ellipsoidal, a variance matching term for compactness distribution can be added to the morphological energy component to further constrain the diversity of cluster shapes to be consistent with the target.

[0183] S344b: Calculate the difference between the cluster centroid distance variation coefficient and the spatial homogeneity index of the target characteristic parameters in the state of the new three-dimensional digital core model, as the spatial matching energy component.

[0184] In one embodiment, the centroid coordinates of all clusters in the new model are calculated, the Euclidean distances between all centroid pairs are collected, and the coefficient of variation of these distances is calculated (i.e., standard deviation divided by the mean plus a small positive number to prevent division by zero). The absolute difference or squared difference between this coefficient of variation and the target spatial uniformity index is used as the spatial matching energy component. This component reflects the uniformity of the stomatal clusters in three-dimensional space; the smaller the coefficient of variation, the more uniform the distribution. A larger deviation results in a higher penalty to correct cluster aggregation or sparseness.

[0185] Specifically, the target space uniformity index is 0.35, which indicates that the centroid distance is relatively consistent. If the coefficient of variation of the new model reaches 0.68, it means that some regions are densely clustered while others are sparse. The increase of this component prompts the adjustment of seed point arrangement or connectivity operation in the optimization to make the distribution more uniform.

[0186] In another implementation, for basalt with slightly clustered pore distribution under deep high-pressure environment, the weight of this component can be appropriately relaxed, but it is still retained to maintain the overall spatial regularity.

[0187] For example, in one embodiment, the target index is 0.42, and the coefficient of variation in the initial stage of the new model is 0.79. This component dominates the increase in total energy, which leads to a decrease in acceptance rate until the centroid distance is normalized through boundary growth and inter-cluster path construction. After the coefficient of variation drops to 0.45, the energy decreases significantly, thereby obtaining a digital core that is more consistent with the geological distribution.

[0188] Preferably, a logarithmic transformation can be performed on the distance set before calculating the coefficient of variation to better handle situations with large scale differences and improve sensitivity to spatial inhomogeneity.

[0189] S345b: The energy components of porosity, cluster distribution, morphology and spatial matching are weighted and summed to obtain the total energy of the new three-dimensional digital core model state. The difference between the total energy of the new three-dimensional digital core model state and the total energy of the current three-dimensional digital core model state is calculated as the energy difference.

[0190] In one embodiment, weights are preset for each component, such as porosity energy weight of 0.4, cluster distribution weight of 0.3, morphology weight of 0.15, and spatial matching weight of 0.15. The total energy is obtained by multiplying each component by its weight and summing the results. The total energy is then calculated for the new model and the current model, and the difference between them is obtained. This is used for subsequent judgments based on the Metropolis acceptance criterion.

[0191] For example, when the new model significantly improves cluster distribution and spatial matching but slightly deviates from the porosity, if the weighted total energy decreases by 0.012, the energy difference will be significant. If the value is negative, the state is accepted directly, and the model is driven to converge toward the optimal convergence of multiple features.

[0192] In one possible implementation, the weights can be adaptively adjusted with each iteration stage, with higher weights for porosity in the early stages and increased weights for morphology and space in the later stages, to achieve phased optimization and improve the final model's overall reflection of the basalt pore structure.

[0193] S35: If energy difference If the energy difference is less than zero, the new 3D digital core model state is accepted; if the energy difference is less than zero, the new state is accepted. If the value is greater than zero, the new three-dimensional digital core model state is accepted according to the probability acceptance criterion; through iterative iteration, the cluster features of the optimized three-dimensional digital core model are made to approximate the target feature parameters.

[0194] In one embodiment, step S35 specifically includes the following steps:

[0195] S51: Accept or reject the new state.

[0196] Determine the energy difference The value of. If If the candidate model is better, then the candidate model is unconditionally accepted as the current model for the next iteration. If the new model indicates a worse state, then the probability is used. Accept the candidate model, where T is the simulated annealing temperature parameter. Generate a uniformly random number in the interval [0,1]. If the random number is less than p, accept the worse new model state; otherwise, reject it and keep the current model state unchanged.

[0197] S52: Update temperature and cycle.

[0198] As the number of iterations increases, the temperature parameter T is gradually reduced according to a predetermined cooling plan, which gradually decreases the probability of the algorithm accepting different solutions in the later stages, thus achieving stable convergence. Steps S31 to S35 are repeated to form an iterative optimization loop. In each loop, the model gradually adjusts the number, size, and spatial relationship of its internal pore clusters by accepting operations that are beneficial to reducing energy (such as splitting super-large clusters and connecting small clusters). After a sufficient number of iterations, the characteristic parameters of the current model, such as cluster size distribution and spatial uniformity, will continuously approach the preset target characteristic parameters, ultimately obtaining a high-fidelity three-dimensional digital core model of vesicular basalt.

[0199] In one embodiment, the process of dynamically adjusting the weights in step S31 can be illustrated by a specific scenario. For example, the target feature parameters show that the average cluster volume of the real core is 1000 voxels, and the maximum cluster volume does not exceed 5000 voxels. In a certain optimization iteration, the current model is found to have a maximum cluster volume of 8000 voxels and an average cluster volume of 1500 voxels, with the number of clusters only 70% of the target value. At this point, the deviation calculated in step S313 clearly indicates the problem of "too large clusters and too few clusters". Step S314 then responds by significantly increasing the weight of the splitting operation from the baseline 20% to 70%, while decreasing the weight of the connectivity operation from 30% to 10%. This makes it highly likely that in subsequent iterations, the algorithm will execute the splitting operation in step S2 to cut those ultra-large clusters with a volume exceeding 5000 voxels, thereby quickly and effectively correcting the current deviation.

[0200] Understandably, the random selection of the cutting plane in step S32 is crucial. For example, for an irregularly shaped, large cluster, a fixed cut along the coordinate axis might not effectively separate its core. Using a random-direction plane cut ensures that, over long-term iterations, regardless of the cluster's shape, it has a chance to be separated from its thicker internal parts, resulting in more natural-looking sub-clusters that more closely resemble the vesicular separation of real basalt. This avoids artificial anisotropy in the 3D digital core model caused by a fixed cutting direction.

[0201] It should be noted that the shortest path principle used in step S33 to construct connected paths has clear technical benefits. Excavating straight or near-straight channels between adjacent clusters minimizes the increase in pore volume while maintaining connectivity, thus maximizing the stability of the overall porosity of the 3D digital core model. Simultaneously, straight channels conform to the physical phenomenon in real rocks where pore throats typically extend along the direction of least resistance, enhancing the physical realism of the 3D digital core model.

[0202] Specifically, the exponential probability acceptance criterion in step S35 is crucial for the Markov chain-Monte Carlo algorithm to escape local optima. For example, in the later stages of optimization, the 3D digital core model may get stuck in a low-energy local state, where any single operation may temporarily increase the energy difference. If you completely refuse Changes will halt optimization. By using... The probability of accepting such upward movement gives it the opportunity to escape the current local depression, explore a lower-energy, better model state region, and ultimately find a better global solution, making the final three-dimensional digital core model more closely approximate the target features.

[0203] S4: Perform curvature smoothing post-processing on the optimized three-dimensional digital core model to obtain a three-dimensional digital core model of vesicular basalt.

[0204] In one embodiment, the process of performing curvature smoothing post-processing on the optimized three-dimensional digital core model specifically includes the following steps:

[0205] S41: Extract the pore boundary voxels of the optimized 3D digital core model, and calculate the difference between the local neighborhood average value and the center value of each boundary voxel as the curvature symbol.

[0206] In one embodiment, step S41 specifically includes the following steps:

[0207] S411: Traverse all voxels in the 3D binary model and identify the set of voxels at the interface between pores and matrix.

[0208] Specifically, if a voxel is a pore and there is at least one matrix voxel in its 26 neighborhood, or if a voxel is a matrix and there is at least one pore voxel in its 26 neighborhood, then the voxel is marked as a boundary voxel.

[0209] S412: For each boundary voxel, define a local neighborhood window, such as a cube window centered on the voxel with a side length of L.

[0210] Calculate the arithmetic mean of the gray values ​​of all voxels within the window, where the gray value of the pore voxels is 1 and the gray value of the matrix voxels is 0.

[0211] S413: Obtain the gray value of the current boundary voxel as the center value.

[0212] S414: Calculate the curvature symbol.

[0213] The curvature sign is equal to the center value minus the local neighborhood average. Understandably, if the center value is 1 (pore) and the neighborhood average is low (surrounded by matrix), the curvature sign is positive, indicating that the pore voxel is located in a locally convex position. Conversely, if the center value is 0 (matrix) and the neighborhood average is high (surrounded by pores), the curvature sign is negative, indicating that the matrix voxel is located in a locally concave position.

[0214] S42: If the curvature sign is negative, the corresponding matrix voxel is converted into a pore voxel for growth; if the curvature sign is positive, the corresponding pore voxel is converted into a matrix voxel for erosion.

[0215] In one embodiment, step S42 specifically includes the following steps:

[0216] S421: Classify all boundary voxels according to the curvature sign calculated in step S41.

[0217] Set a threshold ε close to zero, for example, ε=0.01. If the curvature sign is less than -ε, the voxel is determined to be located in a significantly concave region; if the curvature sign is greater than ε, the voxel is determined to be located in a significantly convex region.

[0218] S422: Perform a growth operation on matrix voxels that are identified as significantly depressed areas.

[0219] Changing the value of this voxel from 0 to 1 transforms it from a matrix into a pore. This operation is equivalent to filling the small pits or grooves at the pore boundaries, making the pore morphology fuller.

[0220] S423: For porosity voxels identified as significantly raised areas, perform an erosion operation.

[0221] Changing the value of this voxel from 1 to 0 transforms it from a pore to a matrix. This operation is equivalent to smoothing out the spikes or protrusions on the pore boundary, making the boundary smoother.

[0222] In one embodiment, to control the smoothing intensity and avoid excessive changes in porosity, a processing scale parameter ρ can be set. For example, the aforementioned growth or erosion operation can be performed only on the top 20% of boundary voxels with the highest absolute values ​​of curvature signs. By adjusting ρ, the boundary smoothing effect can be balanced with the degree of preservation of the overall model characteristics.

[0223] S43: Apply sphere opening and closing operations to the three-dimensional digital core model processed in step S42 to smooth the pore boundaries.

[0224] In one embodiment, step S43 specifically includes the following steps:

[0225] S431: Select a spherical structural element with radius r, for example, r = 2 voxel units.

[0226] S432: First, perform morphological opening operation on the three-dimensional digital core model, that is, corrosion followed by expansion.

[0227] The erosion operation removes all boundary pore voxels that cannot be fully contained within the structural elements, thereby severing and eliminating fine pore-connecting branches and isolated noise points. The subsequent dilation operation restores the size of the main pore structure, but the removed fine branches are not restored.

[0228] S433: Next, perform a morphological closing operation on the result of the opening operation, that is, first dilate and then erode.

[0229] The expansion operation expands the pore region, filling the tiny depressions on the pore surface and the internal micropores. The subsequent corrosion operation shrinks the pore boundaries, but the already filled micro-cavities are preserved. After processing with the sphere opening and closing operation sequence, high-frequency noise and micro-structural defects in the three-dimensional digital core model are effectively suppressed.

[0230] S44: Convert the three-dimensional digital core model processed in step S43 into a floating-point model, perform Gaussian blurring on the floating-point model, and then perform threshold binarization segmentation to obtain a smooth boundary model.

[0231] In one embodiment, step S44 specifically includes the following steps:

[0232] S441: Convert the data type of the three-dimensional digital core model processed in step S43 from integer to floating-point, with each voxel value represented as 1.0 or 0.0.

[0233] S442: Perform three-dimensional Gaussian filtering on the floating-point model.

[0234] Construct a 3D Gaussian convolution kernel, with its standard deviation σ controlling the smoothness. For example, let's take σ = 1.5. Convolve the Gaussian kernel with the floating-point model, so that the originally sharp 0-1 boundary transitions into a gradual grayscale field. At the boundary, the voxel values ​​are between 0 and 1.

[0235] S443: Set a binarization threshold.

[0236] The binarization threshold is typically set to 0.5. In the Gaussian-filtered floating-point model, locations with voxel values ​​greater than or equal to 0.5 are identified as pores and assigned a value of 1; locations with voxel values ​​less than 0.5 are identified as matrix and assigned a value of 0. After this process, the resulting binary model has smoother, more natural pore boundaries, eliminating the step-like artifacts caused by voxelization.

[0237] S45: If there are residual super-large clusters, calculate the distance field for the super-large clusters, identify the local maxima of the distance field as seed points, and obtain the three-dimensional digital core model of vesicular basalt by segmenting the super-large clusters using the watershed algorithm.

[0238] In one embodiment, step S45 specifically includes the following steps:

[0239] S451: Calculate the shortest distance from each pore voxel within the supercluster to the cluster boundary, and use it as the distance field.

[0240] Specifically, the process first identifies ultra-large pore clusters in the current model whose volume exceeds the target maximum cluster size threshold. Then, distance transformation calculations are performed only on all pore voxels within these clusters. The distance transformation uses Euclidean distance metric, finding the straight-line distance from each pore voxel to the nearest cluster boundary voxel and assigning this distance to the corresponding location to form a distance field. Larger values ​​in this distance field indicate proximity to the cluster center, while values ​​close to 0 indicate cluster edges or narrow connections. The distance field obtained through this step provides a basic terrain description for subsequent seed identification and segmentation.

[0241] S452: Identify local maxima in the distance field as seed markers, with each seed marker corresponding to a sub-cluster center.

[0242] For example, by traversing all voxel positions in the distance field, the distance value of each voxel is compared with the distance values ​​of other voxels in its 26 neighborhoods. If the distance value of the current voxel is strictly greater than that of all neighboring voxels, it is determined to be a local maximum point. To avoid the influence of noise, a small-amplitude Gaussian smoothing filter can be applied to the distance field first to make the distance value transition more continuous, and then local maximum detection is performed. Preferably, a minimum distance threshold is set, and only local maxima with distance values ​​greater than a certain value are retained as valid seeds to prevent too many invalid markers from being generated at the cluster edges. Each selected seed marker represents a potential sub-cluster core region, because the peak of the distance field usually corresponds to a thicker independent part within the cluster. This automatic identification method ensures that the seed position objectively reflects the internal topology of the cluster and avoids the segmentation bias that may be caused by random or uniform point scattering.

[0243] In one possible implementation, when a super-large cluster exhibits multiple significantly enlarged regions, the distance field naturally forms multiple peaks, correspondingly generating multiple seed markers. For example, if the cluster consists of three spherical pores connected by narrow channels, the distance field will exhibit local maxima at the center of each spherical portion, thus automatically generating three seeds. This multi-sub-case approach helps to accurately divide the super-large cluster into three sub-clusters, significantly improving the consistency between the cluster size distribution and the actual basalt vesicle statistics.

[0244] S453: Invert the distance field into an elevation field and diffuse it from various sub-markers to the neighborhood simultaneously. When the diffusion areas of different seed marks meet, a dividing boundary is established at the meeting point.

[0245] Specifically, the original distance field is negatively evaluated, or the negative distance transformation result is used directly to create an inverted elevation field. The original distance maxima become the lowest valleys, and the original edge locations become highlands. This simulates the flow of rainwater from highlands to valleys, but in actual segmentation, flooding begins from the marked seed points upwards. The algorithm uses a queue or priority queue approach, initiating diffusion simultaneously from all seed marks and gradually assigning neighboring voxels to the nearest seeds in ascending order of distance value. When diffusion fronts from different seeds meet at a voxel location, a segmentation boundary is set at that location, typically marked as a watershed line or used directly as a barrier between subclusters. This synchronous diffusion mechanism ensures that the segmentation boundary falls on the natural weak connections of clusters, such as narrow necks, thus achieving reasonable morphological separation.

[0246] For example, in a super-large cluster formed by two large pores connected by a narrow channel, after the distance field is reversed, two troughs will be formed corresponding to two seeds. The diffusion process extends upwards from the two troughs simultaneously, eventually meeting at the middle of the channel and establishing a dividing boundary. The boundary obtained in this way is located precisely in the narrowest region of the cluster. The volumes of the two sub-clusters formed after the division are closer to the target average cluster size, while maintaining the original morphological characteristics of the channel and avoiding unnatural cross-sections that may result from random planar cutting.

[0247] It should be noted that the segmentation sensitivity of this diffusion process can be adjusted by controlling the diffusion step size or introducing weights. For example, diffusion can be appropriately suppressed in high curvature regions to preserve certain fine structures. This approach makes the segmentation more adaptable to the irregular characteristics of basalt pores, improving the fidelity of the final digital core model.

[0248] S454: Divide the super-large cluster into different sub-clusters according to the dividing boundary, update the model cluster list, and obtain the three-dimensional digital core model of vesicular basalt.

[0249] After the segmentation boundaries are determined, the voxels within the original super-large clusters are regrouped according to their respective seed tags, with each group forming an independent sub-cluster. Subsequently, the connected component information of the entire 3D digital core model is updated, and the number of clusters and the volume statistics of each cluster in the cluster list are refreshed. Preferably, a slight morphological closure operation is performed on the newly formed segmentation boundaries to smooth them, eliminating possible jagged boundaries and making the sub-cluster edges more consistent with the smooth characteristics of real vesicles. The final output 3D digital core model shows a significant improvement in cluster size distribution; after the super-large clusters are effectively segmented, the overall pore structure is closer to the target features extracted from the original CT scan.

[0250] For example, if an excessively large cluster, more than twice the average size, remains after optimization of the 3D digital core model, the aforementioned watershed segmentation algorithm can automatically divide it into three appropriately sized sub-clusters, thereby reducing the maximum cluster size error from 150% to less than 20%. This automatic segmentation significantly improves the 3D digital core model's ability to reflect the spatial distribution and connectivity of basalt vesicular clusters, providing higher-fidelity digital samples for subsequent seepage simulations.

[0251] Based on the above overall process, the three-dimensional digital core reconstruction method for vesicular basalt provided in this invention extracts target feature benchmarks through adaptive threshold segmentation, and combines multi-stage processing such as Gaussian kernel initial generation, boundary adjustment, adaptive Markov chain optimization, and curvature-driven post-processing to finally generate a high-fidelity three-dimensional digital core model that highly matches the real rock in terms of multi-scale pore structure features. This three-dimensional digital core model can provide a reliable digital foundation for permeability prediction, residual oil distribution simulation, and carbon reservoir sequestration safety assessment in oil and gas development, avoiding the statistical bias in cluster distribution of traditional stochastic reconstruction methods and improving the credibility of numerical simulation results.

[0252] Compared with existing three-dimensional digital core modeling, the present invention has the following advantages:

[0253] 1. This invention requires only a single high-resolution CT scan of a single vesicular basalt sample to generate a large number of three-dimensional digital samples with controllable porosity and adaptive pore structure parameters (cluster size, connectivity, coarseness, etc.) through numerical reconstruction. Compared to the high cost of physical experimental methods, which rely on multiple samples and high-precision scanning equipment, and the additional overhead of existing numerical methods requiring multiple experimental calibrations, this invention significantly reduces the economic and time costs of modeling. It also overcomes the limitations of physical experimental methods, such as the limited number of samples and difficulty in covering diverse pore characteristics, providing ample digital sample support for basalt reservoir research.

[0254] 2. This invention employs a Markov Chain Monte Carlo (MCMC) modeling method with spatial feature constraints on pore clusters. Key features extracted from real core samples, such as the number of clusters, spatial uniformity, and compactness, are used as optimization targets to dynamically adjust the model's iteration direction. Compared to existing numerical methods that do not optimize for the vesicle distribution patterns in basalt, this method significantly reduces the deviation between the spatial distribution of vesicles in digital cores and real cores. It achieves a high-precision reconstruction of the natural morphology, size distribution, and spatial arrangement of original basalt pores, providing a more reliable model foundation for subsequent pore structure evaluation.

[0255] 3. This invention utilizes a post-processing scheme combining curvature feature recognition and multi-scale morphological filtering, along with a strategy for constructing connected paths between adjacent small clusters during the iteration process, to specifically address the core deficiency of existing numerical methods in simulating the complex pore connectivity of basalt. Through directional growth / erosion of pore boundary concave-convex regions, spherical opening and closing operations, and Gaussian smoothing, the model connectivity error is significantly reduced, making the pore network connectivity state of the digital core more closely resemble that of real basalt reservoirs. This provides accurate microscopic pore connectivity support for core applications such as seepage simulation and carbon dioxide sequestration potential assessment.

[0256] 4. This invention, through the comprehensive application of Gaussian kernel deformation to simulate irregular vesicular morphology, energy function-driven iterative optimization, and a refined post-processing workflow, ultimately generates a high-level three-dimensional digital core model of vesicular basalt in terms of morphological realism, statistical feature matching, and spatial structure reproduction. This model not only matches macroscopic porosity but also faithfully reflects the complex geometric topology and spatial configuration of basaltic vesicles at the microscopic scale, thus providing a more reliable and geologically representative three-dimensional computational model for rock physics simulations (such as seepage, electrical conductivity, and mechanical analysis) based on digital cores.

[0257] It should be understood that the various forms of processes shown above can be used to reorder, add, or delete steps. For example, the steps described in this invention disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this invention can be achieved, and this is not limited herein.

[0258] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A method for reconstructing three-dimensional digital cores of vesicular basalt, characterized in that, The process includes the following steps: S1: Extracting target feature parameters from the CT scan grayscale data of the vesicular basalt sample. Target feature parameters include porosity, cluster number, cluster size statistics, spatial uniformity index, and simplified compactness; S2: Generating an initial 3D digital core model based on the target feature parameters; S3: Iteratively optimizing the initial 3D digital core model using an adaptive Markov chain-Monte Carlo algorithm to make the cluster features of the optimized 3D digital core model approximate the target feature parameters. Step S3 specifically includes the following steps: S31: Calculating the deviation between the cluster size of the current 3D digital core model and the target feature parameters, and dynamically adjusting the weights of splitting and connectivity operations; S32... S33: If a super-large cluster is detected, a splitting operation is triggered. A cutting plane is randomly selected within the super-large cluster to divide it into two clusters. S4: If adjacent small clusters are detected, matrix voxels are converted into porosity voxels along the shortest path between adjacent small clusters to construct a connected path and merge the adjacent small clusters. S5: A new 3D digital core model state is generated using an adaptive Markov chain-Monte Carlo algorithm. The energy difference between the new 3D digital core model state and the current 3D digital core model state is calculated using an energy function. S6: If the energy difference is less than zero, the new 3D digital core model state is accepted; if the energy difference is greater than zero, the new 3D digital core model state is accepted according to a probability acceptance criterion. Through iterative iteration, the cluster features of the optimized 3D digital core model are made to approximate the target feature parameters. In step S34, the energy difference between the new 3D digital core model state and the current 3D digital core model state is calculated using an energy function, specifically including the following steps: S341: Calculate the deviation between the porosity of the new 3D digital core model state and the target porosity, as the porosity energy component; S342: Calculate the difference between the number of clusters in the new 3D digital core model state and the number of clusters in the target feature parameters, as well as the difference between the cluster size statistics of the new 3D digital core model state and the cluster size statistics of the target feature parameters, as the cluster distribution energy component. S343: Calculate the difference between the bounding box volume ratio and the simplified compactness of each connected pore cluster in the new 3D digital core model state, as the morphological energy component; S344: Calculate the difference between the cluster centroid distance variation coefficient and the spatial uniformity index of the target feature parameter in the new 3D digital core model state, as the spatial matching energy component; S345: Perform a weighted summation of the porosity energy component, cluster distribution energy component, morphological energy component, and spatial matching energy component to obtain the total energy of the new 3D digital core model state, and calculate the difference between the total energy of the new 3D digital core model state and the total energy of the current 3D digital core model state as the energy difference; S4: Perform curvature smoothing post-processing on the optimized three-dimensional digital core model to obtain a three-dimensional digital core model of vesicular basalt.

2. The method for reconstructing three-dimensional digital cores of vesicular basalt according to claim 1, characterized in that, Step S1 specifically includes the following steps: S11: Perform pore space segmentation on the CT scan grayscale volume data using a custom threshold to obtain a binary pore model; S12: Perform connected component analysis on the binary pore model to identify connected pore clusters, calculate the number of voxels in each connected pore cluster as the cluster size, and collect cluster size statistics, including the maximum cluster volume, minimum cluster volume, and average cluster volume; S13: Calculate the pairwise distances between the centroids of connected pore clusters to obtain a distance set, and use the ratio of the standard deviation of the distance set to the mean of the distance set as a spatial uniformity index; S14: For each connected pore cluster, calculate the axial boundary dimension to obtain the bounding box volume, and use the ratio of the number of cluster voxels to the bounding box volume as the simplified compactness; S15: Use the proportion of pore voxels in the binary pore model as the porosity to obtain the target feature parameter.

3. The method for reconstructing three-dimensional digital cores of vesicular basalt according to claim 2, characterized in that, Step S11 specifically includes the following steps: S111: Read the grayscale data from the CT scan, and determine the dynamic range of grayscale by statistically analyzing the minimum and maximum grayscale values; simultaneously, select the pore phase using the interactive threshold module of Avizo software to obtain a fixed segmentation threshold; S112: Binarize the grayscale data from the CT scan using the fixed segmentation threshold to obtain preliminary pore segmentation results: voxels with grayscale values ​​less than the fixed segmentation threshold are identified as pore voxels, and voxels with grayscale values ​​greater than the fixed segmentation threshold are identified as matrix voxels, thus obtaining a binary pore model; S113: Calculate the pore size of the binary pore model. The percentage of void elements is used as a reference porosity for the segmentation result, and is only used to indicate the porosity level of the current threshold segmentation to the user; S114: Receive the input target porosity as the target parameter; S115: If the pore morphology, connectivity and noise level of the current binary pore model meet the requirements, solidify and output the binary pore model; If the pore morphology, connectivity and noise level of the current binary pore model do not meet the requirements, recalibrate the fixed segmentation threshold using the interactive threshold module of Avizo software, and repeat steps S112 to S114 until a binary pore model that meets the requirements is obtained.

4. The method for reconstructing three-dimensional digital cores of vesicular basalt according to claim 1, characterized in that, Step S2 specifically includes the following steps: S21: Arrange seed points in three-dimensional space according to the preset minimum spacing; S22: Using the seed points as the center, derive the Gaussian kernel standard deviation parameter based on the target cluster size to generate a Gaussian field, and obtain a preliminary model of the natural morphology cluster by threshold binarization segmentation of the Gaussian field; S23: Apply spherical opening and closing operations to the preliminary model of the natural morphology cluster to smooth the pore boundaries; S24: Calculate the difference between the porosity of the preliminary model processed in step S23 and the target porosity, and select candidate voxels at the pore boundaries to perform growth or erosion to adjust the number of pore voxels to the target porosity; S25: Adjust the axial scaling factor of the Gaussian kernel to generate an anisotropic Gaussian field, and obtain an initial three-dimensional digital core model that conforms to the target characteristic parameters through threshold cutting operation based on statistical quantiles and boundary fine-tuning operation.

5. The method for reconstructing three-dimensional digital cores of vesicular basalt according to claim 4, characterized in that, Step S25 specifically includes the following steps: S251: Determine the scaling factor of each axis based on the anisotropy factor; S252: Construct an anisotropic Gaussian kernel function centered on the seed point. The kernel function value is the negative exponent of the sum of the squares of each coordinate in the exponential decay term divided by twice the square of the corresponding axis scaling factor; S253: In the frequency domain, perform a dot product operation between the spectrum of the anisotropic Gaussian kernel and the spectrum of the white noise field, and then generate a continuous Gaussian random field through inverse transformation; S254: Calculate the cumulative distribution function of all voxel values ​​in the Gaussian random field, and use the quantile point obtained by reverse lookup based on the target porosity as the cutting threshold. Voxels higher than the cutting threshold are set as pores to obtain the initial three-dimensional digital core model.

6. The method for reconstructing three-dimensional digital cores of vesicular basalt according to claim 1, characterized in that, Step S4 specifically includes the following steps: S41: Extract the pore boundary voxels of the optimized three-dimensional digital core model, and calculate the difference between the local neighborhood average value and the center value of each boundary voxel as the curvature sign; S42: If the curvature sign is negative, convert the corresponding matrix voxel into a pore voxel for growth; if the curvature sign is positive, convert the corresponding pore voxel into a matrix voxel for erosion. S43: Apply sphere opening and closing operations to the 3D digital core model processed in step S42 to smooth the pore boundaries; S44: Convert the 3D digital core model processed in step S43 into a floating-point model, perform Gaussian blurring on the floating-point model, and then perform threshold binarization segmentation to obtain a smooth boundary model; S45: If there are residual super-large clusters, calculate the distance field for the super-large clusters, identify the local maxima of the distance field as seed points, and segment the super-large clusters using the watershed algorithm to obtain a 3D digital core model of vesicular basalt.

7. The method for reconstructing three-dimensional digital cores of vesicular basalt according to claim 6, characterized in that, Step S45 specifically includes the following steps: S451: Calculate the shortest distance from each pore voxel within the super-large cluster to the cluster boundary, which is used as the distance field; S452: Identify local maxima in the distance field as seed markers, with each seed marker corresponding to a sub-cluster center; S453: Invert the distance field into an elevation field and diffuse it from various sub-markers to the neighborhood simultaneously. When different seed marker diffusion areas meet, establish a dividing boundary at the meeting point; S454: Divide the super-large cluster into different sub-clusters according to the dividing boundary, update the model cluster list, and obtain the three-dimensional digital core model of vesicular basalt.

Images