Method for multi-scale imaging and 3d modeling of digital core in coalbed methane exploration and development

Abstract

This invention provides a method for multi-scale imaging and 3D modeling of digital cores in coalbed methane exploration and development, belonging to the technical field of digital core multi-scale imaging and 3D modeling. The method includes: S1. Multi-scale image acquisition; S2. Image preprocessing; S3. Image feature analysis and optimization; S4. Image segmentation; S5. 3D reconstruction. This invention integrates K-means pre-annotation and U-Net deep learning, combined with surface fitting and MRF global optimization, to capture the morphological features of pores, fractures, and minerals, correct edge jaggedness, solve the problem of misjudgment due to near-identical grayscale values, and achieve accurate reservoir structure segmentation.

Description

Technical Field

[0001] This invention relates to the fields of digital core multi-scale imaging and 3D modeling technology, and in particular to digital core multi-scale imaging and 3D modeling methods in coalbed methane exploration and development. Background Technology

[0002] Traditional coal reservoir image analysis techniques have significant limitations when dealing with complex geological structures:

[0003] Insufficient segmentation accuracy: Relying on a single threshold method or basic clustering algorithm, classifying categories solely by grayscale values ​​makes it difficult to distinguish between micro-fractures and coal matrix, or mineral edges and matrix with similar grayscale values. This leads to missed detection of pores and fractures and misjudgment of minerals, making it impossible to accurately characterize the internal structure of the reservoir.

[0004] High dependence on manual intervention: It requires a large amount of manually labeled training data, and the labeling results are affected by subjective experience, resulting in poor consistency and difficulty in adapting to large-scale sample analysis.

[0005] Model decoupling from reality: When two-dimensional segmentation results are directly used for three-dimensional modeling, the three-dimensional model is easily distorted due to jagged edges and discontinuous pores; the extraction of pore networks ignores the constraints of the natural morphology of pore throats, resulting in large deviations between permeability calculations, fracturing design and actual field conditions.

[0006] Lack of multi-technology collaboration: The segmentation, modeling, and application stages are carried out independently, lacking a data feedback mechanism. As errors accumulate in the preceding steps, the reliability of the final engineering guidance decreases. Summary of the Invention

[0007] This invention provides a method for digital core multi-scale imaging and three-dimensional modeling in coalbed methane exploration and development to solve one or more of the problems mentioned above.

[0008] To achieve the above objectives, the present invention adopts the following technical solution:

[0009] Digital core multi-scale imaging and 3D modeling methods for coalbed methane exploration and development include:

[0010] S1. Multi-scale image acquisition; acquire micron-scale CT images, perform scanning parameter adaptation processing to obtain micron-scale images containing minerals and matrix; acquire nano-scale CT images, perform scanning parameter adaptation processing to obtain nano-scale images containing microcracks and nanopores;

[0011] S2. Image preprocessing; Denoising and grayscale enhancement are performed on the micrometer-level image to obtain the micrometer-level preprocessed image; Denoising and grayscale enhancement are performed on the nanometer-level image to obtain the nanometer-level preprocessed image; Spatial alignment and fusion are performed on the micrometer-level and nanometer-level preprocessed images to obtain the fused preprocessed image;

[0012] S3. Image Feature Analysis and Optimization: Perform grayscale distribution and morphological feature analysis on the fused preprocessed image to obtain feature analysis results; Based on the feature analysis results, normalize the fused preprocessed image to obtain an optimized image;

[0013] S4. Image segmentation; K-means clustering is used for initial segmentation of the optimized image to obtain pre-labeled labels; manually corrected labels are obtained, and a U-Net model is trained based on the pre-labeled labels and manually corrected labels to perform semantic segmentation on the optimized image to obtain a class probability map; surface fitting is performed on the edge pixels in the class probability map to obtain smooth edge curves; combining the class probability map and the smooth edge curves, a Markov random field is used for global optimization to obtain the final segmentation result;

[0014] S5. 3D reconstruction; image sequence alignment is performed on the final segmentation results to obtain an aligned segmentation sequence; the aligned segmentation sequence is stacked along the depth direction and combined with a smooth edge curve for edge optimization to obtain a 3D digital core model.

[0015] This manual also includes the following methods for digital core multi-scale imaging and 3D modeling in coalbed methane exploration and development:

[0016] S6. Pore network extraction; extract the effective pore space from the three-dimensional digital core model to obtain the effective pores; identify the pores and throats and calculate the parameters of the effective pores to obtain the pore throat parameters; construct a pore network model based on the pore throat parameters.

[0017] This manual also includes the following methods for digital core multi-scale imaging and 3D modeling in coalbed methane exploration and development:

[0018] S7. Model Validation and Application: The porosity and permeability of the pore network model are compared and validated with laboratory measurement data to obtain a validated pore network model; Based on the validated pore network model, reservoir evaluation and fracturing scheme design are carried out to obtain reservoir stimulation schemes.

[0019] In this specification, the training of the U-Net model based on pre-labeled labels and manually corrected labels in S4 includes: constructing a training set with pre-labeled labels accounting for 80% and manually corrected labels accounting for 20%, training using a weighted cross-entropy loss function, and stopping training when the cross-union ratio of the validation set is not less than 0.92.

[0020] In this specification, S4 combines the class probability map and the smooth edge curve, and uses a Markov random field for global optimization, including: using the class probability map as a data constraint, the smooth edge curve as an edge constraint, and the class consistency of adjacent pixels as a spatial smoothing term, constructing an energy function, and iteratively optimizing until the energy function converges.

[0021] In this specification, S4, which uses K-means clustering to perform initial segmentation on the optimized image, includes: determining the number of cluster centers to be 3 based on the gray-level distribution features obtained in S3, corresponding to pores, matrix, and minerals respectively; iterating the clustering until the change in the number of cluster centers is less than 5, and then stopping to obtain the pre-labeled labels for the initial segmentation.

[0022] In this specification, S4 involves surface fitting of edge pixels in the category probability map, which includes: extracting edge pixels of pore-matrix and mineral-matrix in the category probability map, fitting the spatial distribution of edge pixels with a quadratic polynomial surface to obtain a continuous and smooth edge curve, so that the edge direction deviation is less than 2°.

[0023] In this specification, the quadratic polynomial surface fitting uses the least squares method to fit the edge pixels and calculates the average distance between the fitted curve and the original edge pixels. When the average distance is less than 3 pixels, it is determined to be a smooth edge curve.

[0024] In this specification, the gradient descent method is used for iterative optimization of the energy function. The pixel category label is updated in each iteration until the difference between the energy function values ​​of two adjacent iterations is less than 0.01, and the global optimal final segmentation result is obtained.

[0025] In this specification, S5 involves stacking the aligned segmentation sequences along the depth direction and combining them with a smooth edge curve for edge optimization. This includes: stacking the aligned segmentation sequences along the depth direction to form an initial three-dimensional model; and correcting the stepped edges of pores and cracks in the three-dimensional model based on the smooth edge curve obtained in S4, so that the fitting error between the edge voxels and the surface is less than 1 voxel.

[0026] In summary, the present invention has at least the following beneficial effects:

[0027] Improved segmentation accuracy: By integrating K-means pre-labeling and U-Net deep learning, and combining surface fitting and MRF global optimization, it can capture the morphological features of pores, fractures and minerals, correct edge jaggedness, solve the problem of misjudgment when grayscale is close to the target, and achieve accurate segmentation of reservoir structure.

[0028] Reduce reliance on manual labor: K-means is used to automatically generate pre-labeled data, reducing the workload of manual labeling. At the same time, algorithmic constraints are used to ensure labeling consistency, adapting to the needs of large-scale sample analysis.

[0029] Enhanced model reliability: The 2D segmentation results are smoothed at the edges and then used for 3D modeling to ensure the continuity of the spatial distribution of pores and minerals; the pore network extraction is combined with natural morphology constraints to make the model more closely match the actual geological characteristics of the reservoir.

[0030] The entire process is collaborative and closed-loop: the segmentation, modeling and application stages are linked through data transmission and feedback mechanisms. Optimization of the preceding steps directly improves the quality of subsequent results, and the final reservoir stimulation scheme is more in line with the field requirements, providing reliable technical support for the efficient development of coal reservoirs. Attached Figure Description

[0031] Figure 1 This is a schematic diagram of the digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development involved in this invention.

[0032] Figure 2 This is a schematic diagram of the three-dimensional reconstructed pore model involved in this invention.

[0033] Figure 3 This is a schematic diagram of the three-dimensional reconstructed mineral model involved in this invention.

[0034] Figure 4 This is a schematic diagram of the three-dimensional reconstructed matrix model involved in this invention. Detailed Implementation

[0035] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0036] like Figure 1 As shown, this embodiment provides a method for multi-scale imaging and 3D modeling of digital cores in coalbed methane exploration and development, including:

[0037] S1. Multi-scale image acquisition; acquire micron-scale CT images, perform scanning parameter adaptation processing to obtain micron-scale images containing minerals and matrix; acquire nano-scale CT images, perform scanning parameter adaptation processing to obtain nano-scale images containing microcracks and nanopores;

[0038] S2. Image preprocessing; Denoising and grayscale enhancement are performed on the micrometer-level image to obtain the micrometer-level preprocessed image; Denoising and grayscale enhancement are performed on the nanometer-level image to obtain the nanometer-level preprocessed image; Spatial alignment and fusion are performed on the micrometer-level and nanometer-level preprocessed images to obtain the fused preprocessed image;

[0039] S3. Image Feature Analysis and Optimization: Perform grayscale distribution and morphological feature analysis on the fused preprocessed image to obtain feature analysis results; Based on the feature analysis results, normalize the fused preprocessed image to obtain an optimized image;

[0040] S4. Image segmentation; K-means clustering is used for initial segmentation of the optimized image to obtain pre-labeled labels; manually corrected labels are obtained, and a U-Net model is trained based on the pre-labeled labels and manually corrected labels to perform semantic segmentation on the optimized image to obtain a class probability map; surface fitting is performed on the edge pixels in the class probability map to obtain smooth edge curves; combining the class probability map and the smooth edge curves, a Markov random field is used for global optimization to obtain the final segmentation result;

[0041] S5. 3D reconstruction; image sequence alignment is performed on the final segmentation results to obtain an aligned segmentation sequence; the aligned segmentation sequence is stacked along the depth direction and combined with a smooth edge curve for edge optimization to obtain a 3D digital core model.

[0042] In some embodiments, acquiring micron-level CT images in S1 and performing scan parameter adaptation processing includes setting the scan voltage to 120kV, the current to 100μA, and the slice thickness to 40μm, so that the resolution of the micron-level CT image reaches 35μm / pixel. Acquiring nanon-level CT images in S1 and performing scan parameter adaptation processing includes setting the scan voltage to 80kV, the current to 50μA, and the slice thickness to 5μm, so that the resolution of the nanon-level CT image reaches 5μm / pixel.

[0043] In some embodiments, S2, the noise reduction and grayscale enhancement processing of the micron-level image includes: first, using nonlocal mean filtering to remove Gaussian noise, then using median filtering to remove salt-and-pepper noise, and finally using adaptive histogram equalization to enhance grayscale differences. S5, the image sequence alignment of the final segmentation result includes: based on the smooth edge curve obtained in S4, performing edge feature matching on adjacent final segmentation results, correcting the scan displacement, and making the alignment error less than 1 pixel. S6, the identification and parameter calculation of effective pores and throats includes: using the maximum sphere method to identify pores and throats, calculating the pore radius, throat length, and coordination number, and removing invalid throats blocked by minerals. S2, the spatial alignment and fusion of the micron-level preprocessed image and the nano-level preprocessed image includes: establishing a coordinate mapping relationship through metal markers on the core surface, superimposing the micro-fracture features of the nano-level preprocessed image onto the corresponding positions of the micron-level preprocessed image, forming a fused preprocessed image containing multi-scale features.

[0044] In some embodiments, the verification of porosity and permeability of the pore network model in S7 by comparing it with laboratory measurement data includes: calculating the porosity and permeability of the pore network model, comparing it with the porosity measured by the helium method and the permeability measured by the pulse decay method, and determining that the verification is passed when the porosity deviation is less than 3% and the permeability deviation is less than 15%. The reservoir evaluation and fracturing scheme design based on the verified pore network model in S7 includes: evaluating the reservoir based on the pore volume, throat radius, and coordination number of the pore network model, determining the fracturing target, construction pressure, and proppant parameters based on the evaluation results, and obtaining the reservoir stimulation scheme.

[0045] In some embodiments, the grayscale distribution and morphological feature analysis of the fused preprocessed image in S3 includes: verifying the three-peak distribution characteristics of pores, matrix, and minerals by statistically analyzing grayscale histograms; extracting pore and fissure edges and analyzing their orientation and length distribution using the Canny operator; and identifying mineral clumps and statistically analyzing their area and morphological parameters using a region growing algorithm. The normalization processing of the fused preprocessed image based on the feature analysis results in S3 includes: determining the minimum and maximum grayscale values ​​according to the grayscale distribution characteristics; mapping the image grayscale values ​​to the [0,1] interval using linear normalization to eliminate grayscale deviations between different scanning batches.

[0046] In some embodiments, extracting effective pore space from the three-dimensional digital core model in S6 includes: screening pores with a volume of not less than 5 voxels and connected paths, and removing isolated pores completely enclosed by minerals to obtain effective pore space. Reservoir evaluation based on the validated pore network model in S7 includes: evaluating reservoir capacity based on pore volume and pore radius, evaluating seepage capacity based on throat radius and coordination number, and evaluating stimulation potential based on mineral distribution characteristics to form a comprehensive reservoir evaluation result. Determining fracturing targets based on the evaluation results in S7 includes: selecting areas with a pore volume ratio greater than 25%, an average coordination number greater than 3.5, and a mineral volume ratio less than 15% as fracturing targets, and verifying the reliability of the targets using well logging curves.

[0047] In some embodiments, a feedback correction step is also included: if the model verification in S7 fails, the corresponding step is backtracked for correction, wherein when the porosity deviation exceeds the limit, the effective pore extraction in S6 is corrected, and when the permeability deviation exceeds the limit, the throat identification in S6 is corrected, until the verification passes and a closed loop is formed.

[0048] In some embodiments, the scanning parameter adaptation process in S1 further includes: performing slice thickness calibration for micron- and nanometer-scale CT images to ensure that the physical depth difference error between adjacent image slices is less than 2 μm, providing a basis for the spatial continuity of subsequent 3D reconstruction. The parameters for nonlocal mean filtering are set as follows: window size 7×7, similarity threshold 0.1, and weighted averaging of pixels with grayscale differences less than the threshold, preserving pore and fissure edge features while removing Gaussian noise. Adaptive histogram equalization uses an 8×8 local window and a contrast limit threshold of 2.0, increasing the grayscale difference between pores and the matrix from 10-20 to 20-30, enhancing the separability of K-means clustering. During edge optimization, voxels deviating from smooth edge curves in the initial 3D model are corrected by category, prioritizing the retention of pore voxels along the surface direction, improving the continuity of pore and fissure direction in 3D space by 40%.

[0049] In some embodiments, the maximum sphere method is used to identify pores and throats. The pore radius is the shortest distance from the center of the sphere to the nearest matrix or mineral voxel, and the throat is the smallest sphere in the path connecting adjacent pores. The throat length is the distance between the centers of the two pores minus the radii of the two pores. In reservoir capacity evaluation, the standard for a high-quality reservoir is: a total pore volume greater than 5 × 10⁻⁶. 7 μm 3 The average pore radius is greater than 50 μm; in the evaluation of seepage capacity, the standard for a high-quality seepage channel is: an average throat radius greater than 30 μm and an average coordination number greater than 3. Specifically, the logging curve verification is as follows: the logging curve corresponding to the fracturing target point must meet the following requirements: a gamma value lower than 60 API and a resistivity higher than 50 Ω·m, ensuring that the target point is a high-quality coal reservoir area.

[0050] In some embodiments, during the retrospective correction process, if the porosity deviation exceeds the limit, the effective porosity screening threshold in S6 needs to be readjusted (e.g., the minimum volume is reduced from 5 voxels to 3 voxels); if the permeability deviation exceeds the limit, micro-throats with a diameter less than 5 μm need to be identified and the throat parameters recalculated. In S1, when acquiring micron-level and nanon-level CT images, the same core sample is scanned in sections. The micron-level scan covers the entire core, while the nanon-level scan focuses on the high-porosity development areas identified in the micron-level scan, achieving targeted optimization of data acquisition. Metal markers are 0.1 mm diameter tungsten beads, with at least three non-collinear markers evenly arranged on the core surface. Alignment is achieved by matching the coordinates of the markers in the two scales of the image and calculating the spatial transformation matrix, with an alignment error of less than 2 pixels. When the region growth algorithm identifies mineral agglomerates, pixels with a gray value greater than 180 are used as seed points. The growth criterion is that the gray value difference between adjacent pixels is less than 10, and the area of ​​the grown region is greater than 100 pixels, indicating a valid mineral agglomerate.

[0051] In some embodiments, the initial centers of K-means clustering are set as follows: the median of the gray-level range of 50-80 for pore centers, the median of the gray-level range of 80-180 for matrix centers, and the median of the gray-level range of 180-255 for mineral centers, accelerating the clustering convergence speed. The weighted cross-entropy loss function assigns a 1.5x weight to mineral categories and a 1.2x weight to pore categories, addressing the issue of a small proportion of mineral and pore samples in the training data. Edge feature matching uses the Hausdorff distance to measure the similarity of edge curves in adjacent layers. When the distance is greater than 3 pixels, alignment is achieved through translation and rotation, ensuring continuous tracking of pores and fissures in the depth direction. Connectivity path determination employs a three-dimensional region labeling algorithm, marking each pore as a connected region and retaining only pores with connected paths to the model boundary, ensuring that the extracted effective pores participate in actual seepage.

[0052] In some embodiments, when using the maximum sphere method to traverse the pore space, the sphere center interval is set to 5 μm to ensure that pores and throats with diameters greater than 10 μm are not missed, while keeping the computational load within an acceptable range. Porosity is calculated as the ratio of the sum of the volumes of all pores in the pore network model to the total volume of the three-dimensional digital core model. Permeability is calculated based on Darcy's law, obtained through fluid flow simulation of the pore-throat network. In determining the proppant parameters, the proppant particle size is less than 0.8 times the minimum throat radius, and the compressive strength is greater than 1.2 times the reservoir closure pressure, ensuring that the proppant can effectively support fractures without clogging throats.

[0053] In some embodiments, during micrometer-level CT image scanning, a 180° rotation scanning mode is used, acquiring one frame of image every 0.5° rotation. Two-dimensional slice images are generated through a tomographic reconstruction algorithm to improve the spatial resolution of the image. During nanometer-level CT image scanning, a high-resolution detector is used, and the exposure time is set to 500ms to reduce image noise and improve the grayscale contrast of the nanopores by 20%.

[0054] In some embodiments, the weight calculation of nonlocal mean filtering is based on the gray-level similarity of the pixel neighborhood. The higher the similarity, the greater the weight, so that the gray-level gradient retention rate of the pore and crack edges after filtering is greater than 80%. The spatial smoothing term of the Markov random field sets the penalty coefficient when adjacent pixels are of different categories. The penalty coefficient is 1.0 when pores are adjacent to the matrix and 1.5 when minerals are adjacent to the matrix, which enhances the spatial consistency of the edges.

[0055] In some embodiments, during the evaluation of fracturing potential, when the proportion of minerals filling the fractures is greater than 30%, the area is identified as having high fracturing difficulty, requiring an increase in fracturing pressure of 10%–15%. When the coefficient of variation of the throat radius is greater than 0.5, a segmented fracturing method is adopted. When supplementing the identification of micro-throats, an adaptive spherical radius algorithm is used to reduce the spherical center interval to 2 μm in densely porous areas, ensuring that throats with a diameter of 3–5 μm are effectively identified.

[0056] In some embodiments, after the 3D model is optimized, a volume consistency check is performed. The deviation between the total pore volume of the 3D model and the sum of the pore volumes of each layer in the 2D segmented image is calculated. If the deviation is greater than 2%, edge optimization is performed again. If the logging curve of the fracturing target does not meet the gamma value and resistivity conditions, the target position is adjusted within an adjacent 5cm range based on the lateral distribution characteristics of the 3D digital core model to ensure that the target falls in a high-quality reservoir area. When removing invalid throats blocked by minerals, a 3D spatial query is used to determine whether the throat passes through a mineral area. If the overlap length between the throat and the mineral voxel is greater than 30% of the total throat length, it is determined to be an invalid throat. After the image sequence alignment is completed, an alignment quality report is generated, including the edge matching degree of adjacent layers and the cumulative displacement of the overall sequence. If the cumulative displacement is greater than 10 pixels, the entire sequence alignment is performed again.

[0057] The technical concept of this invention is as follows:

[0058] Coal Sample Collection: Coal samples were drilled from target coal seams (such as C7 and C8 seams) in the Longtan / Xuanwei Formation in southern Sichuan. Sampling depths were determined based on the exploration targets (e.g., 760.32–760.62 m in Well Wei-1 and 686.50–686.55 m in Well FB-1). Sampling specifications included: full-diameter cores: 5–10 cm in diameter and 5–10 cm in length (for macroscopic structural observation); plunger samples: 2.5 cm in diameter and 2.5–5 cm in length (for micron-scale CT scanning); small sample pieces: 1–2.5 cm in diameter (for nanometer-scale CT scanning). Pretreatment: Cleaning: Removing mud and impurities from the coal sample surface to avoid contaminating imaging; Drying: Drying in a 60℃ oven to constant weight (to avoid moisture affecting grayscale values); Marking: Marking depth and direction (e.g., axial and radial) on the coal sample surface to ensure consistent imaging direction. Processing Results: Clean, dry, and fully marked coal samples were obtained, providing a basic sample for subsequent multi-scale imaging.

[0059] S1. Multi-scale image acquisition

[0060] Core objective: To provide "low-noise, high-fidelity" raw data for subsequent high-precision segmentation and modeling of S4-S7, and to solve the problems of "grayscale blurring, artifact interference, and large scale differences" in CT images of coal reservoirs in southern Sichuan.

[0061] S1.1 Data Acquisition (Adapting to Subsequent Segmentation Requirements)

[0062] To address the segmentation challenge of "large differences in pore and crack scale (nanometer to millimeter)" in S4, a "multi-scale scanning scheme" is required:

[0063] Micron-level CT scanning: Covering millimeter-scale minerals and matrix (scanning parameters: voltage 120kV, current 100μA, slice thickness 40μm, resolution 35μm / pixel), ensuring complete imaging of mineral blocks (>100μm); Equipment: MicroXCT-400 μ-CT scanner; Parameters: voltage 170kV, power 17W, resolution 35μm, slice thickness 40-50μm;

[0064] NanoCT scanning: Focusing on microcracks and nanopores (scanning parameters: voltage 80kV, current 50μA, slice thickness 5μm, resolution 5μm / pixel), capturing microcracks <50μm (key features that U-Net needs to learn in S4); Equipment: Phoenix Nanotom nanoCT; Parameters: resolution 65-150nm, scanning range for small samples with a diameter <1mm;

[0065] Data alignment: Micron / nano CT images are fused by core markers (e.g., pasted metal markers) to ensure consistent spatial coordinates of multi-scale data in the same area (error < 2 pixels), providing a basis for cross-scale segmentation of S4.

[0066] S1.2 Data Quality Check (Predicting Subsequent Interference Factors)

[0067] After data collection, three types of invalid data need to be removed (to avoid affecting the accuracy of S4 segmentation):

[0068] Noisy images: If the grayscale standard deviation is >30 (normal range 10-20), it is judged as noise pollution (such as caused by vibration during scanning) and needs to be marked for noise reduction;

[0069] Artifact images: Check for "ring artifacts" (CT detector error) or "blurred edges" (insufficient scan dose). Images with artifacts accounting for more than 5% need to be rescanned.

[0070] Scale mismatch: When the deviation of the same aperture edge between nano-CT and micro-CT is greater than 5 pixels, it is necessary to realign by marking points.

[0071] S2. Image preprocessing (to reduce interference for S4 segmentation)

[0072] Core objective: To remove image noise and artifacts, enhance the gray-level difference between pores and the matrix, and provide high-quality input for K-means clustering (dependent on gray-level distribution) and U-Net segmentation (dependent on edge features) in S4.

[0073] S2.1 Noise Removal (Adapted to S4 Edge Preservation Requirements)

[0074] The coal reservoir images in southern Sichuan contain "Gaussian noise (electronic noise) + salt-and-pepper noise (scanning interference)," requiring a "tiered noise reduction strategy."

[0075] Step 1: Nonlocal mean filtering (Gaussian noise removal)

[0076] Preserve the edges of the pores and cracks (S4 surface fitting requires continuous edges), parameter settings: window size 7×7 (balancing noise reduction and edge preservation), similarity threshold 0.1 (only smoothing areas with similar gray levels), reducing the noise standard deviation from 25 to 8;

[0077] Step 2: Median filtering (to remove salt-and-pepper noise)

[0078] For isolated noise points (which are easily misidentified as pores in S4 by K-means), a 3×3 median filter is used to remove pixels with gray-level abrupt changes >50 (such as salt and pepper noise points), while avoiding excessive blurring of edges (edge ​​gray-level changes are retained >30).

[0079] S2.2 Artifact Correction (Improves Grayscale Consistency)

[0080] Ring artifact correction: The detector response bias is corrected by polynomial fitting—the radial gray mean of the image is calculated, and the fitted curve eliminates the "concentric circle gray fluctuations", reducing the gray standard deviation of the same matrix region from 15 to 5 (ensuring the gray stability of K-means clustering in S4).

[0081] Edge artifact removal: crop the image edges by 10 pixels (including edge bright bands caused by scanning equipment) to prevent U-Net in S4 from misidentifying artifacts as mineral edges.

[0082] S2.3 Gray-scale enhancement (enhancing the identification of pores and cracks)

[0083] To address the issue of "mineral and matrix grayscale overlap" in S4, "Adaptive Histogram Equalization (CLAHE)" is employed:

[0084] The local window size is 8×8 (to adapt to the micro-crack scale), and the contrast limit threshold is 2.0 (to avoid excessive noise enhancement), which increases the gray level difference between the pores (gray level 50-80) and the matrix (80-180) from 10-20 to 20-30, providing a clearer classification boundary for gray level clustering of K-means in S4.

[0085] S3. Image Feature Analysis and Optimization

[0086] Core objective: Based on the S2 preprocessing results, extract key features (grayscale, morphology, distribution) of pores, fissures, and minerals to provide a basis for the selection of the S4 algorithm and further optimize image quality.

[0087] S3.1 Feature Analysis (Guiding S4 Algorithm Parameters)

[0088] 1. Gray-scale distribution analysis

[0089] Statistically analyze the gray-level histogram of the preprocessed image to verify whether it conforms to the "tri-peak hypothesis" of K-means clustering in S4: if the peak spacing of the pore peak (50-80), matrix peak (80-180), and mineral peak (180-255) is >30 (e.g., 60 in actual measurement), then K-means can be directly clustered; if the peaks overlap (spacing <20), it is necessary to return to S2 to enhance the gray-level difference (e.g., adjust the CLAHE parameter).

[0090] 2. Morphological Feature Analysis

[0091] Pore ​​and fracture: Potential edges are extracted through edge detection (Canny operator), and the length (>5 pixels is considered valid) and orientation (e.g., NE30° is the main orientation in well FB-1) are statistically analyzed to provide a basis for the design of convolution kernels in U-Net in S4 (3×3 convolution is used to capture linear extension features).

[0092] Minerals: Mark the area (>100 pixels is considered valid) and shape (blocky / stripy) of mineral clumps to ensure that U-Net in S4 can learn their blocky distribution characteristics.

[0093] S3.2 Preprocessing Optimization (Connecting to S4 Input)

[0094] Multi-scale image fusion: Nanoscale CT (microcracks) and microscale CT (minerals) images are aligned by marker points to generate a "multi-resolution superimposed image"—the microscale image preserves the minerals and matrix, while the nanoscale image supplements the microcrack details, providing complete input for multi-scale feature extraction of U-Net in S4;

[0095] Normalization: Normalize the image grayscale to [0,1] ( This eliminates grayscale deviations between different scanning batches, ensuring the consistency of U-Net training data in S4.

[0096] S3.3 Connection with subsequent steps

[0097] Output: Optimized image (Noise removal, grayscale enhancement, normalization) and feature analysis report (grayscale distribution, pore and crack morphology);

[0098] Support for S4: Directly used as input to K-means and U-Net, the grayscale-enhanced image increases the spacing between K-means cluster centers in S4 by 40%, reducing pre-labeling errors; the feature analysis results guide the adjustment of U-Net network parameters (such as adapting the convolution kernel size to the micro-crack morphology).

[0099] Coordination logic between S1-S3 and S4-S7

[0100] Data-driven chain: S1’s multi-scale scanning provides full-scale targets (nanometer to millimeter) for S4 segmentation; S2’s noise reduction enhancement ensures grayscale stability of K-means clustering in S4 (cluster center deviation < 5); S3’s feature analysis guides S4 algorithm parameters (such as U-Net convolution kernel size).

[0101] Feedback mechanism: If edge recognition is blurred in S4, backtrack to S3 to check the grayscale enhancement effect (e.g., whether the CLAHE parameter needs to be adjusted); if the 3D reconstruction in S5 finds discontinuous pores, return to S2 to verify whether noise removal is excessive (edges are blurred).

[0102] Targeted optimizations in S1-S3 lay the foundation for high-precision segmentation, modeling, and application in subsequent S4-S7, forming a complete closed loop of "data acquisition-preprocessing-segmentation-modeling-application".

[0103] S4. Image segmentation (fusion of U-Net, K-means, Markov random field, and surface fitting)

[0104] Core Logic and Objectives

[0105] CT images of coal reservoirs in southern Sichuan present three major challenges: ① Microfractures are extremely small (width < 50 μm) and have only a slight difference in grayscale value from the coal matrix (grayscale difference < 20); ② Minerals (such as pyrite) overlap with the matrix grayscale (some minerals have a grayscale value of 80-150, which coincides with the matrix); ③ The edges of pores and fractures appear jagged due to imaging noise. Traditional thresholding methods rely solely on grayscale values, which easily leads to "missed detection of microfractures and misjudgment of mineral boundaries."

[0106] This solution achieves precise segmentation through a "four-level collaborative process":

[0107] 1. K-means clustering: Utilizes the characteristics of gray-scale distribution to initially classify categories and generate low-cost pre-labeled data (reducing manual labeling by 80%).

[0108] 2. U-Net Semantic Segmentation: Based on pre-labeled data, it learns the morphological features of pores and fractures (such as the linear extension of microfractures and the blocky distribution of minerals) and outputs pixel-level class probabilities.

[0109] 3. Surface Fitting: Smooth the edges identified by U-Net to solve the jagged edge problem and generate continuous boundary curves;

[0110] 4. Markov Random Field (MRF): It combines the probabilities of U-Net (to ensure local correctness) with the boundary constraints of surface fitting (to ensure global morphological rationality) to finally output high-precision segmentation results.

[0111] S4.1 K-means Clustering: Initial Segmentation and Pre-labeling

[0112] Function: The differences in the physical properties of "porosity (low density, gray level 50-80), coal matrix (medium density, gray level 80-180), and minerals (high density, gray level 180-255)" in coal reservoirs result in a "trimodal gray level distribution." K-means automatically classifies these three types of targets through clustering, providing pre-labeled labels for U-Net and reducing the cost of manual labeling.

[0113] 1.1 Model Construction and Input / Output

[0114] Input: S3 preprocessed grayscale image (size H×W, pixel grayscale values) (Noise and artifacts have been removed). For example, a CT slice image is 512×512 pixels and contains 100 microcracks and 20 mineral grains.

[0115] Core assumption: Image grayscale distribution can be divided into three significant clusters (pores, matrix, and minerals), and the grayscale values ​​of targets of the same type are concentrated in a certain range.

[0116] Output: Cluster centers: (Typical gray level of pores, such as 60) (Typical gray level of the matrix, such as 130) (Typical mineral gray level, such as 210); Pre-labeled tags: (0 = pores, 1 = matrix, 2 = minerals), with dimensions consistent with the input image.

[0117] 1.2 Specific Process

[0118] Step 1: Initialize cluster centers (based on prior knowledge)

[0119] Based on statistical data on coal reservoirs in southern Sichuan, the initial center was set as follows:

[0120] (Initial center of pores) (Initial center of matrix) (Initial mineral centers) The initial values ​​are based on the actual data distribution to avoid clustering converging to a local optimum.

[0121] Step 2: Pixel allocation (categorized by grayscale distance)

[0122] For each pixel (i,j), calculate its absolute distance to the three centers (the smaller the distance, the more likely the pixel belongs to that class):

[0123] (Distance to the center of the pore);

[0124] (Distance to the center of the matrix);

[0125] (Distance to the mineral center);

[0126] t is the iteration number, initially t=0;

[0127] Pixel attribution rules:

[0128] like and ,but (Porosity);

[0129] like and ,but (Matrix);

[0130] like and ,but (mineral).

[0131] For example: If the grayscale value of a pixel is g=70, the calculation is as follows: =20, =50, =130, then it is classified as pore (L=0).

[0132] Step 3: Update cluster centers (adjust within-cluster means)

[0133] Calculate the center (average grayscale value of all pixels within the class) in the (t+1)th iteration:

[0134] ; ;

[0135] ;

[0136] In the formula: The total number of pixels in class k in the t-th iteration (e.g.) (Number of pixels with pores).

[0137] For example: if the pore class contains 1000 pixels and the total grayscale value is 60000, then .

[0138] Step 4: Convergence Judgment and Pre-label Generation

[0139] Repeat steps 2-3 until the change in all centers is less than 1. (k=1,2,3), at which point the center is stable (denoted as). ).

[0140] Based on the center, the segmentation threshold is determined, and the final pre-label is generated:

[0141] ,like ; ,like ; ,like .

[0142] 1.3 Interaction with U-Net (Data Transfer)

[0143] Pre-labeling As one of the training labels for U-Net (accounting for 80%), it constrains the model's learning through U-Net's loss function—for example, when U-Net predicts a pixel as a hole (probability 0.9)... When labeled as matrix, the loss function increases significantly, forcing the model to adjust its parameters to reduce such errors.

[0144] S4.2U-Net Semantic Segmentation: High-Precision Class Probability Prediction

[0145] Function: K-means relies solely on grayscale values, resulting in a high misclassification rate (approximately 30%) for "micro-cracks (grayscale close to the matrix) and mineral edges (grayscale gradients)". U-Net learns the spatial morphological features of the image (such as the linear extension of micro-cracks and the blocky boundaries of minerals) and outputs the probability (0-1) of each pixel belonging to three target classes, providing a precise basis for subsequent optimization.

[0146] 2.1 Model Construction (Network Structure and Parameters)

[0147] Input: S3 preprocessed grayscale image (normalized to [0,1], denoted as...) +K-means pre-labeling (80%) + Manually corrected labels (20%).

[0148] Network structure: U-Net (Encoder-Decoder symmetric structure, adapted to multi-scale features of coal reservoir images)

[0149] Encoder (feature extraction): Contains 4 convolutional blocks, each block consisting of "2 layers of 3×3 convolution (stride 1, padding 1) + ReLU activation + 2×2 max pooling (stride 2)".

[0150] Input layer: 256×256×1 (grayscale image); Convolution block 1: Output 128×128×64 (extracting shallow features such as micro-fracture edges); Convolution block 2: Output 64×64×128 (extracting mid-layer features such as pore distribution); Convolution block 3: Output 32×32×256 (extracting deep features such as mineral blocks); Convolution block 4: Output 16×16×512 (extracting overall lithological distribution features).

[0151] Decoder (resolution restoration): Contains 4 deconvolution blocks, each block consists of "1 layer of 2×2 deconvolution (stride 2) + fused Encoder corresponding layer features + 2 layers of 3×3 convolution + ReLU activation".

[0152] Deconvolution block 1: Output 32×32×256 (fusing deep and mid-level features); Deconvolution block 2: Output 64×64×128 (fusing mid-level and shallow features); Deconvolution block 3: Output 128×128×64 (fusing shallow and edge features); Deconvolution block 4: Output 256×256×3 (3 channels corresponding to the probabilities of three types of objects). Output layer: 1×1 convolution + Softmax activation, output probability map. (k=0 / 1 / 2 corresponds to pores / matrix / minerals, and) ).

[0153] 2.2 Model Training (Integrating K-means Pre-labeling)

[0154] Dataset Construction: 2000 CT slices from three wells in southern Sichuan (FB-1, L1-1, and Wei-1) were selected, covering typical types such as primary structured coal and fractured coal. 80% of the samples were pre-labeled using K-means. (Reducing labor costs); 20% of the samples had labels manually corrected. (Manual correction of error-prone areas such as micro-cracks and mineral edges to ensure accuracy in critical areas).

[0155] Loss function (weighted cross-entropy, distinguishing between pre-labeled and human-labeled labels): To balance the weights of the two types of labels (human-labeled labels are more reliable), the loss function is defined as follows:

[0156] ;

[0157] In the formula: (Manual labels have higher weight); If pixel (i,j) actually belongs to category k (0 / 1 / 2), then 0; If the K-means pre-labeled pixel (i,j) belongs to category k, then 0; Predict the probability (output value) that pixel (i,j) belongs to class k for U-Net.

[0158] For example: the artificial label of a micro-crack pixel (i,j) (Porosity), if U-Net predicts (Low probability), then the loss is (The punishment is severe).

[0159] Training process: Optimizer: Adam (learning rate = 1e-4, decay rate = 0.9, to avoid overfitting); Iteration stopping condition: Intersection over Union (IoU) of the validation set ≥ 0.92 (IoU = overlap area between predicted and true labels / total coverage area); Output: The trained U-Net model can directly output a probability map from a new image. .

[0160] 2.3 Interaction with preceding and subsequent algorithms

[0161] Input dependency: K-means pre-labeling As training data, through the loss function The constrained model learns the morphological characteristics of pores, fractures, and minerals;

[0162] Output propagation: Probabilistic graph It serves as input for surface fitting (extracting edge pixels) and also provides the basis for the "data items" of the MRF (ensuring correct local classification).

[0163] S4.3 Surface Fitting: Edge Smoothing and Boundary Constraints

[0164] Function: The edges of holes and cracks output by U-Net may appear "jagged" due to pixel discreteness (such as pixel jumps at the edge of micro-cracks). By using surface fitting, a continuous and smooth edge curve is generated, which provides geometric constraints for MRF (ensuring that the edges conform to natural shapes and avoiding jaggedness).

[0165] 3.1 Model Construction and Input / Output

[0166] Input: U-Net probabilistic graph "Pore-matrix edge pixel set" Specific extraction rules: Filter all pixels (i,j) that satisfy... (A porosity probability close to 0.5 indicates that the pixel is located at the transition edge between the pore and the matrix.) This set contains the discrete pixel coordinates of the edge. (m is the number of edge pixels, usually ≥100, to ensure the reliability of the fit).

[0167] Core assumption: The edges of natural pores and fractures in coal reservoirs are continuous and smooth (not discrete serrated), and can be fitted with low-order polynomial surfaces (such as quadratic polynomials), which can capture the edge trend and avoid overfitting by high-order polynomials.

[0168] Output: Equation of smooth edge surface and coefficient matrix (Edge constraints used in MRF). Wherein, The theoretical boundary between pores and matrix ( For pores, (as a matrix).

[0169] 3.2 Detailed Calculation Process (Quadratic Polynomial Fitting)

[0170] Step 1: Extracting edge pixel coordinates

[0171] from Extract the two-dimensional coordinates of the pixels (x is the horizontal pixel index, y is the vertical pixel index) to obtain the coordinate set. For example, if the edge of a microcrack contains 150 pixels, then m = 150, and the coordinate range is typically... (H and W represent the image height and width).

[0172] Step 2: Define the surface equation (quadratic polynomial)

[0173] Choose a quadratic polynomial (to balance fitting accuracy and smoothness, and avoid oscillations caused by noise in higher-order polynomials): ;

[0174] In the formula: The coefficients to be determined (collectively referred to as the coefficient matrix) ).

[0175] Step 3: Solve for coefficients using least squares (to make edge pixels fit the curved surface)

[0176] Objective: To make edge pixels appear on the surface The sum of squared distances is minimized (i.e., the fitting error is minimized). For each edge pixel... Ideally Therefore, we construct a system of equations:

[0177] Represented as a matrix: ;

[0178] in: for Design matrix (each row corresponds to one edge pixel): ; for Zero vector (ideal margin error).

[0179] Solve the coefficient matrix using least squares (to make...) Minimum):

[0180] ;

[0181] Step 4: Generate smooth edge constraints

[0182] The solution obtained Substituting into the surface equation, we get For any pixel If it belongs to the edge region (i.e. ), then through The categories of sign judgment theory: if Theoretically, it is porosity (category 0); if Theoretically, it is a matrix (Category 1). Minor deviations are permissible. (pixels), if If the deviation is not corrected, it is considered an edge deviation and needs to be corrected in the subsequent MRF.

[0183] 3.3 Interaction with preceding and subsequent algorithms

[0184] Input Dependency: Probabilistic Graph of U-Net Determine the edge pixel set (Pixels with a probability close to 0.5) directly affect the object of surface fitting;

[0185] Output transfer: Surface equations and coefficient As an input to the "edge constraint term" of the MRF (see S4.4), it forces the edge pixel class to conform to the smooth surface.

[0186] S4.4 Markov Random Field (MRF): Global Optimization and Final Segmentation

[0187] Function: Combining the local probabilities of U-Net (ensuring class correctness) and the global edge constraints of surface fitting (ensuring morphological rationality), it corrects two types of problems: 1. Isolated noise points (such as single pixels misclassified as pores by K-means / U-Net); 2. Jagged edges (such as edge pixel jumps in the output of U-Net).

[0188] The final output is a segmentation result that is spatially continuous and has a reasonable shape.

[0189] 4.1 Model Construction (Energy Function Definition)

[0190] Core variables: The final category of pixel (i,j) is represented by (0=pore, 1=matrix, 2=mineral).

[0191] Energy function: The lower the energy, the more reasonable the segmentation result (satisfying both data consistency and spatial continuity): ;

[0192] In the formula: (Data item): Constrained segmentation is consistent with U-Net probabilities; (Spatial smoothing term): Constrains adjacent pixels to have the same category; (Edge Constraint Term): Constraint edge pixel fitting surface fitting results; (Spatial smoothing weights) (Edge constraints have higher weights, prioritizing edge shape).

[0193] 4.2 Calculation of each energy term (integration of U-Net and surface fitting)

[0194] 1. Data Items (Depending on U-Net probabilities)

[0195] Ensure the segmentation results are consistent with the high-probability categories predicted by U-Net:

[0196] ;

[0197] For example: if U-Net predicts the probability that pixel (i,j) is a hole. Then when hour, (Low energy); if misjudged as matrix (c=1), then Possibly (Energy levels increased significantly).

[0198] 2. Spatial smoothing term (Adjacent pixel constraints)

[0199] Based on the assumption that "adjacent pixels have a higher probability of belonging to the same class," the penalty for class mutation is:

[0200] ;

[0201] In the formula: It represents the 8-neighborhood (top, bottom, left, right, and 4 diagonal directions) of pixel (i,j). For indicator functions: if but (Energy increases), otherwise .

[0202] For example: if a pixel is a pore (c=0), but its 8-neighborhood contains 5 matrix elements (c=1), then this contribution... (High energy level, requires correction).

[0203] 3. Edge constraint terms (Depending on surface fitting results)

[0204] Force edge pixel category to fit the surface :

[0205] ;

[0206] In the formula Edge indicator function: if s < 0 (theoretical porosity) and ,but If s > 0 (theoretical matrix) and ,but ;like (Allowable deviation), then .

[0207] 4.3 Optimization Process (Iterative Conditional Mode ICM)

[0208] Step 1: Initialize the category

[0209] Initial class settings based on U-Net probabilistics (high probability priority):

[0210] ;

[0211] Step 2: Pixel-by-pixel iterative optimization

[0212] For each pixel (i,j), iterate through all possible classes (0,1,2) and calculate the energy when the class is assumed to be k. Update the category with the lowest energy:

[0213] ;

[0214] For example: the U-Net aperture probability of a certain edge pixel is 0.5, and the surface fitting... (Theoretical porosity), if the current category is matrix (c=1), then (High energy), after iteration it will be corrected to porosity (c=0).

[0215] Step 3: Convergence and Output

[0216] Repeat step 2 until all pixel categories no longer change. For all (i,j) true, output the final segmentation result. .

[0217] 4.4 Interaction with Pre-algorithms

[0218] Probabilistic graphs that directly rely on U-Net This ensures that the segmentation results are consistent with the deep learning predictions;

[0219] Dependence on surface fitting Force edge pixels to fit the smooth curve;

[0220] Output The data is transferred to S5 3D reconstruction as the basis for subsequent modeling.

[0221] Summary of Algorithm Fusion Relationships

[0222] 1. K-means → U-Net: Through pre-labeling Provide training data, and the formula is related to the U-Net loss function. 2. U-Net → Surface Fitting: Through Probability Plots Define edge pixel set 3. Surface fitting → MRF: through definition 4. U-Net → MRF: through constraining edge shape; definition This ensures that the local categories are correct.

[0223] Through four-level collaboration, the micro-fracture recognition rate of coal reservoir images in southern Sichuan was improved by 40%, and the edge smoothness was improved by 50%, providing high-precision input for subsequent 3D modeling and reservoir evaluation.

[0224] S5. 3D Reconstruction (Based on S4 Optimized Segmentation Results)

[0225] Core objective: To convert the high-precision 2D segmented images output by S4 into 3D digital core models, accurately reconstructing the spatial distribution characteristics of pores, matrix, and minerals in coal reservoirs (including the orientation of pores and fractures, and the aggregation patterns of minerals), providing a "three-dimensional data source" for subsequent pore network extraction and reservoir evaluation. The 3D reconstructed pore model, mineral model, and matrix model are shown below. Figure 2 , Figure 3 , Figure 4 As shown.

[0226] S5.1 Input Data and Preparatory Steps

[0227] Input data: 1. The final segmentation result output by S4 (Sequence of 2D images): N layers in total (e.g., 1539 layers in a CT scan), each layer is H×W (e.g., 2024×2024 pixels), pixel values ​​0=pores, 1=matrix, 2=minerals; 2. Scan metadata: layer thickness (e.g., 40μm, depth difference between two adjacent layers), scanning direction (along the core axis); 3. Edge parameters of S4 surface fitting (Used for 3D edge optimization).

[0228] Pre-check: Verification The integrity of the sequence (no missing layers, no misordering) and the presence of slight core displacement during the scanning process are preliminarily determined by the interlayer edge features (such as the micro-crack edges identified by S4). (If the misalignment between the edges of two adjacent layers is greater than 3 pixels, it needs to be marked for correction.)

[0229] S5.2 Specific Implementation Steps

[0230] 1. Alignment of two-dimensional segmented image sequences

[0231] Because core samples may undergo slight rotation or translation during CT scanning, interlayer misalignment needs to be corrected through "feature matching": Feature extraction: for each layer Extract the edge curves of the S4 surface fitting. (Edge of layer z); Matching and alignment: Based on layer 1, calculate the overlap between the edge of layer z and the edge of layer z-1 (using curve similarity algorithm). If the overlap is <80%, align the edges by translating / rotating layer z (maximum adjustment ≤5 pixels); Effect verification: After alignment, the pores and mineral boundaries of adjacent layers are continuous (e.g., if a microcrack extends along the x-axis in layer z, layer z+1 should continue this direction with a deviation ≤1 pixel).

[0232] 2. Generation of 3D volume data

[0233] The aligned 2D images are stacked along the depth direction (z-axis) to generate 3D volume data. The voxel value is defined as: ; This is the two-dimensional segmentation result of the z-th layer;

[0234] Spatial scale transformation: Converting pixel coordinates to actual physical coordinates (e.g., 1 pixel = 35μm, then the actual position in the x-direction is...). Actual depth in the z-direction Ensure that the model size is consistent with the actual core (e.g., if the core length is 50mm, the model z-axis length error is <0.1mm).

[0235] 3. 3D model optimization (based on S4 edge constraints)

[0236] After stacking two-dimensional segmented images, the edges of pores and cracks in the three-dimensional model may appear "stepped" due to small interlayer errors (such as discontinuous orientation of microcracks in the depth direction). Edge constraint correction using S4 surface fitting is required.

[0237] Edge voxel recognition: from 3D volume data Extract voxels from the pore-matrix edge (i.e., voxels whose adjacent voxels are pores and matrix respectively), and label them as... ;

[0238] Smoothing Correction: Surface Equations Based on S4 ,right Adjust the voxels in the middle - if a certain edge voxel deviates If the corresponding three-dimensional extension direction (such as a sudden change of >5° along the z-axis) is corrected, its category is modified (prioritizing the fit to the surface trend).

[0239] Mineral distribution verification: Examining the spatial distribution of minerals in the 3D model ( (The region) should be consistent with the two-dimensional mineral distribution trend of the S4 segment (e.g., if a mineral cluster extends in a "northwest-southeast" direction in the two-dimensional image, the direction should be continued in the three-dimensional model to avoid "isolated mineral clusters").

[0240] 4. 3D model visualization and quality inspection

[0241] Visualization: The optimized 3D model was visualized using Avizo software. Render the image and observe the connectivity of pores and fractures (such as whether a pore extends continuously along the depth direction) and the aggregation morphology of minerals (such as whether they fill along the fractures) through "slice viewing" (along the xz and yz planes).

[0242] Quality Indicators: Edge Continuity: Deviation of pore-matrix edge orientation in 3D space ≤ 2° (based on S4 surface constraints); Volume Consistency: Deviation between the total pore volume of the 3D model and the total volume calculated from the 2D segmented image < 2%. ).

[0243] S5.3 Output and Pre-processing Interaction

[0244] Output: Optimized 3D digital core model (Including physical coordinate information), and model quality report (indicators such as edge continuity and volume consistency).

[0245] Interaction with S4: S4's Provide basic segmentation data and surface fitting for edges. Used for 3D edge optimization, it improves the continuity of pores and cracks in the depth direction by 40% (compared to unconstrained reconstruction).

[0246] Support for subsequent steps: It provides accurate three-dimensional pore space for S6 pore network extraction, avoiding "missed throats" caused by model errors (such as the 5% throat identification error in traditional reconstruction, which can be reduced to less than 1%).

[0247] S6. Pore Network Extraction (Based on Optimized 3D Model)

[0248] Core objective: To extract quantitative parameters (such as radius, length, and connectivity) of pores and throats from the three-dimensional digital core model of S5, and to construct a pore network model (PNM) that reflects the true seepage characteristics of coal reservoirs, providing a "digital twin" for subsequent permeability calculation and reservoir evaluation.

[0249] S6.1 Input Data and Pre-analysis

[0250] Input data: 1. S5 optimized 3D model 2. Mineral distribution results of S4 (Used to identify ineffective pores blocked by minerals); 3. Porosity measured in the laboratory. (Used to verify the extraction results).

[0251] Preliminary analysis: By observing slices of a 3D model, the dominant direction of porosity development can be predicted.

[0252] The dominant direction of pore development (such as extending along the direction of tectonic stress) directly affects the connectivity path and seepage characteristics of pore throats. It needs to be determined by multi-directional slicing of a three-dimensional model to provide directional guidance for the subsequent extraction of pore throats in S6.2 (avoiding blind searching).

[0253] Specific implementation methods

[0254] 1. Multi-directional slice extraction

[0255] 3D model optimized from S5 Three types of slices were extracted: horizontal slices (xy plane): one slice was extracted every 50 layers along the depth direction (z-axis) (a total of 20 slices were extracted) to observe the distribution of pores on the horizontal plane (e.g., whether they are densely distributed along the NE30° direction); vertical slices (xz plane): one slice was extracted every 50 pixels along the y-axis to observe the extension of pores in the vertical-horizontal combination direction (e.g., whether they are continuous along the z-axis depth direction); lateral slices (yz plane): one slice was extracted every 50 pixels along the x-axis to verify whether pores have the characteristic of "preferential development along a certain oblique direction (e.g., x+z direction)".

[0256] 2. Indicators for Determining Advantageous Directions

[0257] For the pores in each slice ( We conducted statistics on the region and used three indicators to determine the advantageous directions:

[0258] Elongation length: The maximum extension distance of pores in different directions (0°, 30°, 45°, 60°, 90°) is measured (e.g., in a horizontal slice, the average extension length of pores in the NE30° direction is 800 μm, which is significantly greater than 500 μm in other directions).

[0259] Number of pores: Statistically calculate the percentage of pores within a 10° range in each direction (e.g., 35% of pores within a 30°±5° range is the highest value);

[0260] Edge conformation consistency: combining edge parameters fitted with S4 surface (The slope of the edge curve reflects the direction). If 80% of the edge direction is concentrated between 25° and 35°, it further verifies that the direction is the dominant direction.

[0261] 3. Verification based on geological background

[0262] The coal reservoirs in southern Sichuan are influenced by Yanshanian tectonic stress, and pores (especially tectonic fractures) often develop along the "NE30°-45°" direction (consistent with the direction of the regional principal stress). If the dominant direction observed in the slice (such as NE32°) matches the geological background, the judgment is reliable; if the deviation is >10° (such as the observation is EW direction), the slice needs to be re-examined (it may be due to misjudgment caused by local mineral obstruction).

[0263] Output and Function: Output a report on the dominant pore direction, clearly stating the core conclusions (e.g., "The pores in the C8 coal seam of well FB-1 mainly extend along the NE32° direction, with weak connectivity in the vertical direction (z-axis)"). This conclusion is used to guide throat extraction in S6.2.

[0264] When searching for connecting paths between pores, prioritize expansion along the dominant direction (NE32°) to reduce invalid searches by 50% (e.g., avoid blindly searching for throats along the vertical direction).

[0265] When calculating throat length, a "connectivity weight" is assigned to throats along the dominant direction (e.g., length error is allowed to be ±5μm, while non-dominant directions are strictly controlled to ±3μm) to ensure the extraction accuracy of critical seepage paths.

[0266] S6.2 Specific Implementation Steps

[0267] 1. Extraction of effective pore space

[0268] From the 3D model of S5 Separate pores, but remove "ineffective pores" (such as isolated pores encased in minerals, or noisy pores that are too small):

[0269] Pore ​​screening criteria:

[0270] Volume ≥ 5 voxels (corresponding to actual volume ≥ (To avoid noise)

[0271] There are connecting paths with other pores (isolated pores completely surrounded by matrix / minerals are eliminated through three-dimensional connectivity analysis).

[0272] Operation method: Use Avizo's "Connected Components" tool to mark all connected pore clusters, retain the top 95% of pore clusters by volume (remove tiny isolated pores), and generate effective pore space. .

[0273] 2. Maximum sphere method for extracting pores and throats (core step)

[0274] Based on the principle of "maximum inscribed sphere", pores (large spheres) and throats (small spheres between pores) are identified, adapting to the pore and fracture morphology of the 3D model:

[0275] Pore ​​identification:

[0276] exist Iterate through all possible sphere centers and calculate the distance from each sphere center to the nearest matrix / mineral boundary (i.e., the maximum sphere radius). The "locally largest sphere" (a sphere not contained within other larger spheres) is retained as the porosity representative:

[0277] ;

[0278] Let the coordinates of the sphere center of the i-th pore be... It is a three-dimensional Euclidean distance.

[0279] For example: if the shortest distance from the center of a pore to the matrix boundary is 58 μm, then its radius... Corresponding pore volume .

[0280] Throat identification:

[0281] In the interconnected paths between pores, the sphere with the smallest radius is defined as the throat (the "bottleneck" controlling fluid flow). For connected pores i and j, the throat radius... The minimum sphere radius on the path, and the throat length. The distance between the centers of the two pores minus their radii:

[0282] ;

[0283] For example: pore i ( ) and porosity j ( If the distance between the centers of the spheres is 210 μm, then the throat length is... .

[0284] 3. Phrevical throat parameter statistics and network construction

[0285] Core parameter calculation:

[0286] Pore ​​parameters: Total number of pores (For example, the C8 coal seam in well FB-1 has 139,885 pores), average pore radius (e.g., 58.58μm), pore volume distribution (divided into intervals such as 40~60μm and 60~80μm according to radius, and the volume percentage of each interval is statistically analyzed);

[0287] Laryngeal parameters: Total number of larynxes (Determined by pore connectivity, e.g., the number of throats in well FB-1 is approximately 120,000), average throat radius. (e.g., 35μm), average laryngeal length (e.g., 87.14μm);

[0288] Connectivity parameter: coordination number (Number of throats connected to a single pore, such as the average coordination number of 3.46 in well FB-1), percentage of connected volume (the proportion of the total volume of connected pores to the total pore volume, such as 65.14%).

[0289] Pore ​​network model construction:

[0290] The pore-throat relationship is represented by a "node-edge" structure—nodes represent pores (with the attribute being...). , , The edge represents the throat (with the attribute of...). , ), forming a complete network topology. .

[0291] Visual verification: Render the pore network using Avizo software to visually observe the pore throat distribution (such as whether there is a "main pore-branch throat" structure) and connectivity paths (such as whether a seepage channel is formed that runs through the model).

[0292] 4. Model optimization (combined with S4 mineral distribution)

[0293] Removing mineral-blocked throats: Mineral distribution based on S4 segmentation If a certain throat passage passes through a mineral voxel ( If any of these are found, they are marked as "invalid larynx" and removed (such larynxes cannot actually participate in seepage).

[0294] Corrected pore throat radius: If the pore edge is close to the mineral (distance < 5 μm), its effective radius needs to be reduced by the part blocked by the mineral (e.g., the original radius is 58 μm, and it is blocked by the mineral by 3 μm, so it is corrected to 55 μm) to ensure that the parameters reflect the actual seepage space.

[0295] S6.3 Output and Interaction with Previous and Next Steps

[0296] Output: Porous network model The database includes parameter tables (including pore throat radius, length, coordination number, etc.), pore throat distribution histograms (e.g., the pore radius of well FB-1 is mainly concentrated in 40~80μm), and three-dimensional pore network visualization diagrams.

[0297] Interaction with S5: S5's It provides accurate pore space boundaries, reducing the pore throat identification error from ±10% in traditional methods to ±3%;

[0298] Support for S7: The parameters (such as throat radius and coordination number) are the core basis for S7 permeability calculation and fracturing design - for example, low coordination number (<2) regions need to improve connectivity through fracturing.

[0299] S7. Model Validation and Application (Based on Optimized Model)

[0300] Core objective: To verify the reliability of the pore network model through laboratory data, and then to formulate targeted coalbed methane reservoir stimulation schemes based on the model parameters (such as fracturing target selection and proppant selection), and predict the development effect, thereby achieving a closed loop of "model-guided engineering".

[0301] S7.1 Model Validation (Multi-dimensional Comparison with Laboratory Data)

[0302] Model validation is a crucial step in ensuring the reliability of the pore network model (PNM). It requires comparison with laboratory data across three dimensions: porosity, permeability, and pore throat distribution, as detailed below:

[0303] 1. Porosity verification

[0304] Model porosity calculation: Based on S6's PNM, the total volume of all effective pores is calculated. Divide by the total volume of the 3D model (Calculated from the dimensions of the S5 model, such as) ):

[0305] ;

[0306] Comparison with laboratory data: The laboratory used the helium method to measure core porosity. (For example, the measured value of coal seam C8 in well FB-1 is 3.098%), requiring .

[0307] Example: PNM calculates the porosity to be 3.05%, compared to... The deviation is 0.048%, which meets the requirements.

[0308] If the deviation is too large (e.g., >5%), backtrack to step S6: check whether micropores (<5 voxels) were missed or mineral-encapsulated pores were misjudged, and re-extract the effective pore space.

[0309] 2. Penetration verification

[0310] Permeability calculation based on PNM: The "pore-throat network seepage model" is adopted, assuming laminar flow of the fluid (coalbed methane) within the pores. Permeability is derived using Darcy's law.

[0311] ;

[0312] In the formula: methane density ( ), The acceleration due to gravity (9.8 m / s²) 2 ); The viscosity of methane ( A is the cross-sectional area of ​​the model (calculated from the dimensions of the S5 3D model).

[0313] Comparison with laboratory data: The laboratory used the pulse decay method to measure core permeability. (For example, the measured value of well FB-1 is 1.92mD), requiring (Relaxed standards for low-permeability reservoirs).

[0314] Example: PNM calculates the permeability to be 1.8 mD with a deviation of 6.25%, which meets the requirements;

[0315] If the deviation is too large (e.g., >20%), backtrack to S6 throat extraction: check whether small throats (<5μm) are missed or the throat length is misjudged, and correct the throat boundary in conjunction with S4 edge constraint.

[0316] 3. Verification of pore throat distribution

[0317] Plot a histogram of pore radius distribution of PNM (e.g., the proportion of the 40~60μm and 60~80μm intervals) and compare it with the pore throat distribution measured by mercury intrusion porosimetry (mercury intrusion porosimetry can reflect the actual pore throat size distribution).

[0318] The deviation of the main range (accounting for >5%) should be <10% to ensure that the throat size of the PNM is consistent with the actual size.

[0319] S7.2 Reservoir Evaluation (Based on PNM Parameters)

[0320] Based on the validated PNM, the reservoir is evaluated from three dimensions: "reservoir capacity", "flow capacity", and "stimulation potential", providing a basis for subsequent development.

[0321] 1. Storage capacity evaluation

[0322] Key metric: Total pore volume Average pore radius ;

[0323] Evaluation criteria:

[0324] High-quality reservoirs: (Corresponding to a core length of 5cm) (Sufficient pore space);

[0325] Example: C8 coal seam in well FB-1 , Its storage capacity is above average.

[0326] 2. Seepage capacity evaluation

[0327] Key metric: Average throat radius Mean coordination number ;

[0328] Evaluation criteria: High-quality seepage channels: (The throat is not easily blocked). (Good connectivity);

[0329] Example: FB-1 well , It has good permeability, but there are local low-permeability areas with throats <20μm (which need to be improved by fracturing).

[0330] 3. Evaluation of Transformation Potential

[0331] Based on the S4 mineral distribution: if the minerals fill along the fractures ( If the overlap between the minerals and the pore edges is greater than 30%, then the area will be difficult to modify (high-pressure fracturing is required).

[0332] Throat heterogeneity: Calculating the coefficient of variation of the larynx radius ( (Standard deviation of throat radius) The area needs to be fracturing in stages (to avoid uneven distribution of proppant).

[0333] S7.3 Reservoir Stimulation Scheme Development (Based on Model Parameters)

[0334] Based on the reservoir evaluation results, a targeted fracturing stimulation plan is developed, which mainly includes three parts: "target selection", "construction parameter design", and "proppane selection".

[0335] 1. Selection of fracturing target points

[0336] Preferred region criteria: Pore volume fraction > 25% (sufficient storage space); average coordination number (Easily forms interconnected networks after fracturing); Mineral volume percentage <15% (reduces fracturing fluid loss);

[0337] Location method: Using the three-dimensional coordinates of PNM, mark the highly connected areas (such as the 760.42-760.48m depth section of the C8 coal seam in well FB-1), and verify the reliability of the target point by combining the logging curves (low gamma value, high resistivity).

[0338] 2. Construction parameter design

[0339] Construction pressure: Calculate the rupture pressure (the minimum pressure required to open the throat) based on the throat radius:

[0340] ;

[0341] In the formula: (Original reservoir pressure). (Displacement) (Cross-sectional area), calculated as follows ;

[0342] Fracturing fluid dosage: designed to be 1.5 times the pore volume ( ), such as the required fluid volume at the FB-1 well target point. .

[0343] 3. Prop selection

[0344] Particle size selection: proppant diameter < 0.8 × minimum throat radius (to avoid blockage), such as PNM with a minimum throat of 20μm, select 12μm ceramic proppant;

[0345] Strength requirements: compressive strength > reservoir closure pressure (calculated from the pore pressure of PNM, such as 8.5MPa), select medium-strength ceramsite (compressive strength 10MPa).

[0346] S7.4 Development Effect Prediction and On-site Feedback

[0347] 1. Gas production prediction

[0348] A seepage model was established based on PNM to simulate gas production after fracturing:

[0349] Input parameters: pore throat parameters of PNM, fracture permeability after fracturing (increased by 5 times);

[0350] Predicted results: FB-1 well target point stable daily gas production 1200~1300 m³ / h 3 The cumulative gas production in the first six months was 210,000 m³. 3 .

[0351] 2. On-site feedback and model correction

[0352] If the measured construction pressure (e.g., 9.6 MPa) deviates from the calculated value (9.4 MPa) by less than 5%, the model is verified as reliable; if the initial gas production pressure is lower than the predicted value (e.g., 800 m³ / min), the model is verified as reliable. 3 / d), check the throat connectivity of PNM (some throats may be missed), and adjust the subsequent fracturing sections in conjunction with production logging (uneven production profile).

[0353] S7.5 Output Results

[0354] Final outputs: 1. Reservoir stimulation plan (including target coordinates, construction pressure, and proppant parameters); 2. Model validation report (comparison results of porosity and permeability); 3. Gas production prediction curve and on-site adjustment suggestions.

[0355] Application value: The optimized scheme reduces the risk of sand blockage during fracturing of FB-1 well by 40% and increases gas production efficiency by 35%, providing quantitative support for the efficient development of coalbed methane in southern Sichuan.

[0356] Overall process summary: From high-precision segmentation in S4 to engineering application in S7, each step is coordinated through a closed-loop data process: the segmentation accuracy of S4 determines the quality of the 3D model in S5, the reliability of the model in S5 affects the extraction of pore throat parameters in S6, and the PNM in S6 directly guides the fracturing design in S7. Through the optimization of the entire process of "segmentation-modeling-verification-application," a precise connection from laboratory core samples to field development has been achieved, providing a replicable technical paradigm for coalbed methane exploration and development.

Claims

1. A method for multi-scale imaging and 3D modeling of digital cores in coalbed methane exploration and development, characterized in that, include: S1. Multi-scale image acquisition; Micrometer-scale CT images are acquired and scan parameters are adapted to obtain micrometer-scale images containing minerals and matrix. Nanoscale CT images are acquired and scan parameters are adapted to obtain nanoscale images containing microcracks and nanopores. S2. Image preprocessing; Denoising and grayscale enhancement are performed on micrometer-level images to obtain micrometer-level preprocessed images; Denoising and grayscale enhancement are performed on nanoscale images to obtain nanoscale preprocessed images; Spatial alignment and fusion are performed on the micron-level preprocessed image and the nanometer-level preprocessed image to obtain the fused preprocessed image; S3. Image Feature Analysis and Optimization; Gray-scale distribution and morphological features of the fused preprocessed image are analyzed to obtain feature analysis results; Based on the feature analysis results, the fused preprocessed image is normalized to obtain an optimized image; S4. Image segmentation; K-means clustering is used to perform initial segmentation on the optimized image to obtain pre-labeled labels. Specifically, based on the physical properties of coal reservoir porosity, matrix, and minerals and the gray-level three-peak distribution characteristics, the pixel gray values ​​of the optimized image are clustered. Pre-labeled labels are generated through initial center setting, pixel distance allocation, iterative update of cluster centers, and convergence judgment. To obtain manually corrected labels, a U-Net model is trained based on pre-labeled labels and manually corrected labels. Semantic segmentation is performed on the optimized image to obtain a category probability map. Specifically, the optimized image is used as input, and a training set is constructed by combining pre-labeled labels and manually corrected labels. An Encoder-Decoder symmetric network structure is adopted, and the training is constrained by a weighted cross-entropy loss function to output a category probability map. To obtain a smooth edge curve, surface fitting is performed on the edge pixels in the category probability map. Specifically, the pore-matrix edge pixel set is extracted from the category probability map. Based on the assumption of natural continuity of the pore and fracture edges in the coal reservoir, a quadratic polynomial is used for surface fitting. The smooth edge curve is generated by solving the coefficient matrix by least squares. By combining the class probability map and the smooth edge curve, a Markov random field is used for global optimization to obtain the final segmentation result. Specifically, the energy function containing data terms, spatial smoothing terms, and edge constraint terms is constructed using the class probability map and the smooth edge curve as inputs. The function is then optimized through iterative conditional mode until the class is stable, and the final segmentation result is output. S5. 3D reconstruction; image sequence alignment is performed on the final segmentation results to obtain an aligned segmentation sequence; the aligned segmentation sequence is stacked along the depth direction and combined with a smooth edge curve for edge optimization to obtain a 3D digital core model.

2. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 1, characterized in that, Also includes: S6. Pore network extraction; Effective pore space is extracted from the three-dimensional digital core model to obtain effective pores; pores and throats are identified and parameters are calculated for the effective pores to obtain pore-throat parameters; a pore network model is constructed based on the pore-throat parameters.

3. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 2, characterized in that, Also includes: S7. Model Validation and Application: The porosity and permeability of the pore network model are compared and validated with laboratory measurement data to obtain a validated pore network model. Based on the validated pore network model, reservoir evaluation and fracturing scheme design are carried out to obtain reservoir stimulation schemes.

4. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 1, characterized in that, In S4, training the U-Net model based on pre-labeled labels and manually corrected labels involves: constructing a training set with 80% pre-labeled labels and 20% manually corrected labels, using a weighted cross-entropy loss function for training, and stopping training when the cross-union ratio of the validation set is not less than 0.

92.

5. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 1, characterized in that, In S4, the global optimization using Markov random fields combines the class probability map and the smooth edge curve. This includes: using the class probability map as a data constraint, the smooth edge curve as an edge constraint, and the class consistency of adjacent pixels as a spatial smoothing term to construct an energy function, and iteratively optimizing until the energy function converges.

6. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 1, characterized in that, In S4, the initial segmentation of the optimized image using K-means clustering includes: determining the number of cluster centers to be 3 based on the gray-level distribution features obtained in S3, corresponding to pores, matrix, and minerals respectively; iterating the clustering until the change in the number of cluster centers is less than 5, and then stopping to obtain the pre-labeled labels of the initial segmentation.

7. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 1, characterized in that, In S4, surface fitting of edge pixels in the category probability map includes: extracting edge pixels of pore-matrix and mineral-matrix in the category probability map, and using quadratic polynomial surface fitting to fit the spatial distribution of edge pixels to obtain a continuous and smooth edge curve, so that the edge direction deviation is less than 2°.

8. The method for digital core multi-scale imaging and three-dimensional modeling in coalbed methane exploration and development according to claim 7, characterized in that, The quadratic polynomial surface fitting method fits the edge pixels using the least squares method and calculates the average distance between the fitted curve and the original edge pixels. When the average distance is less than 3 pixels, it is determined to be a smooth edge curve.

9. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 5, characterized in that, The energy function is iteratively optimized using the gradient descent method. The pixel category label is updated in each iteration until the difference between the energy function values ​​of two adjacent iterations is less than 0.01, thus obtaining the globally optimal final segmentation result.

10. The digital core multi-scale imaging and three-dimensional modeling method for coalbed methane exploration and development according to claim 1, characterized in that, In S5, the alignment segmentation sequence is stacked along the depth direction and combined with the smooth edge curve for edge optimization, including: stacking the alignment segmentation sequence along the depth direction to form an initial three-dimensional model, and correcting the stepped edges of the pores and cracks in the three-dimensional model based on the smooth edge curve obtained in S4, so that the fitting error between the edge voxel and the surface is less than 1 voxel.

Images