3D Deep Learning-Based Segmentation and Morphological Quantification for CT Imaging of Gas-Water Displacement Interfaces

This chapter details the segmentation-reconstruction-quantification integrated method for CT gas-water interfaces. The method description is as follows: 1) Preprocess the CT sequence images and improve the 3D Attention U-net segmentation to obtain binary images of the gas-water interface; 2) Extract a dense point cloud using multi-feature fusion from the segmented interface images, and perform three-dimensional reconstruction using Transformer matching; 3) Calculate multi-dimensional morphological quantification indices based on the reconstructed 3D model.

2.1 CT gas-water interface segmentation based on the improved 3D Attention U-net

2.1.1 Improved 3D Attention U-net network structure

The improved 3D Attention U-net follows the classic U-shaped architecture, with encoding, decoding, and skip connections as the basic framework. The core innovation lies in the introduction of the PAM and the optimization of the skip connection mechanism, which enhances segmentation accuracy of the gas-water interface by feature enhancement and cross-scale fusion. This architecture extracts abstract features through downsampling in the encoding area, restores spatial resolution through upsampling in the decoding area, and compensates for detail loss caused by downsampling via skip connections. PAM and the optimized connection mechanism specifically address the issue of weak features in the gas-water interface, which are prone to interference from the core matrix background.

Figure 1. Improved 3D Attention U-net network structure for CT gas-water interface segmentation

Figure 1 shows the improved 3D Attention U-net network structure suitable for CT gas-water interface segmentation. The encoding area consists of four 3D convolution blocks and three 3D max-pooling layers alternately. Each convolution block contains two 3×3×3 convolution layers, a batch normalization (BN) layer, and a ReLU activation function. In the figure, the BN layers are simplified and omitted. The convolution blocks gradually enhance the texture and morphological features of the gas-water interface by extracting local features. The max-pooling layers perform downsampling with a 2×2×2 kernel size and a stride of 2, which compresses the feature map size while preserving key structural information. The final output feature map size from the encoding area is 1/8 of the input, effectively extracting high-level abstract features. The decoding area consists of four 3D deconvolution blocks, each consisting of a 2×2×2 deconvolution layer, two 3×3×3 convolution layers, a BN layer, and a ReLU activation function. The deconvolution gradually upsamples to restore the spatial resolution of the feature map and provide structural support for precise segmentation.

PAM is the core enhancement module of the network, embedded between the outputs of the encoding area and the attention-guided AG module. It strengthens target feature representation and suppresses background interference through feature re-projection and attention weighting. The specific process is as follows: First, the 3D feature map F_enc∈R^C×H×W×D output by the encoding area, where C is the number of channels, and H, W, and D are the spatial dimensions, is subjected to A/Z/Y/X multi-dimensional 3D average pooling to extract feature information from different spatial dimensions. Then, the pooled features from each dimension are concatenated into a multi-dimensional 3D pooled feature map, which is projected to a lower-dimensional space using a 1×1×1 3D convolution to reduce computational complexity. After BN and ReLU activation function processing, the feature representation ability is enhanced. The processed features are then separated into multiple branches through 3D feature separation, and each branch generates the corresponding attention weights for each dimension through a 1×1×1 3D convolution and a Sigmoid activation function. Finally, the attention weights from each branch are element-wise multiplied with the 3D feature map of the CT gas-water interface output from the original encoding block, resulting in F_att∈R^C×H×W×D. The AG module aligns F_att with the upsampled feature map F_dec from the decoding area through channel alignment, focusing on the target region via attention weight allocation and outputting the fused feature F_fusion. The optimized skip connection no longer directly transmits the original features from the encoding area but fuses them with the features from the decoding area after enhancement via PAM. This effectively reduces interference from background noise during the fusion process and enhances the effectiveness of cross-scale feature fusion. Figure 2 shows the structure of the PAM for CT gas-water interface 3D segmentation.

Figure 2. Structure of the PAM for CT gas-water interface 3D segmentation

2.1.2 Improved hybrid loss function

Gas-water interface segmentation faces two core challenges: 1) The significant pixel ratio difference between the interface region and the core matrix in CT images [17], resulting in severe class imbalance; 2) The small number of boundary pixels and the gradual gray-scale variation [18], making it difficult for traditional loss functions to accurately capture boundary features. To address these two issues simultaneously, a hybrid loss function L_mix with a weight parameter λ is designed, which combines cross-entropy loss and Dice loss through weighted fusion to balance class distribution and boundary segmentation accuracy.

The cross-entropy loss L_CE primarily addresses the class imbalance issue by penalizing misclassified pixels through logarithmic probability. It assigns a higher misclassification cost to the interface pixels, which occupy a very small proportion. The calculation formula is as follows:

$L_{C E}=-\frac{1}{N} \sum_{i=1}^N\left[y_i \log \left(p_i\right)+\left(1-y_i\right) \log \left(1-p_i\right)\right]$                         (1)

where, N is the total number of pixels, y_i is the true label of the i-th pixel, and p_i is the model's predicted probability that the pixel belongs to the interface. The Dice loss L_Dice enhances boundary segmentation accuracy by measuring the overlap between the predicted and true regions, effectively focusing on the subtle features of the interface boundary. The calculation formula is as follows:

$L_{ {Dice }}=1-\frac{2 \sum_{i=1}^N y_i p_i+\epsilon}{\sum_{i=1}^N y_i^2+\sum_{i=1}^N P_i^2+\epsilon}$                        (2)

where, ϵ = 10⁻⁵ is a smoothing term to avoid division by zero in extreme cases.

The final form of the hybrid loss function is: L_mix = λ·L_CE+(1−λ)·L_Dice, where the weight parameter λ is used to adjust the contribution of the two types of losses. The optimal λ value is determined by performing a 5-fold cross-validation over the range λ∈[0,1], with the final optimal value set to λ = 0.3. This value allows the cross-entropy loss to sufficiently suppress the bias caused by class imbalance while enabling the Dice loss to fully optimize boundary features, achieving the best balance between interface pixel recognition and boundary morphology characterization.

2.1.3 Dataset construction and preprocessing

The experimental data are sourced from the Micro-CT scanning results of sandstone core gas-water displacement physical simulation experiments. The spatial resolution of the scanning equipment is 50 μm, with a scan step of 1 mm and a time interval of 30 seconds, resulting in dynamic scan data for 100 time steps. Each time step consists of 200 two-dimensional CT slices, with each slice having a size of 512×512. The data covers the entire dynamic process of gas-water displacement from the initial stage to the stable stage, providing rich dynamic feature samples for the interface segmentation model. The data labeling was independently completed by two experts with more than 5 years of experience in rock mechanics research using the LabelMe3D tool. The labeled object is the gas-water interface region in each slice. The intersection of the labeling results from both experts was taken as the gold standard for segmentation, ensuring the accuracy and authority of the labeling results.

The purpose of data preprocessing is to improve image quality, unify data distribution, and increase sample size, thereby providing high-quality input for model training. First, Gaussian filtering is applied for denoising. A filter kernel with σ = 1.0 is selected to smooth image noise while retaining the subtle gray-scale changes of the interface. This parameter is determined by comparing the gray-scale contrast between the interface and background at different σ values, achieving a balance between denoising and feature preservation. Subsequently, a gray-scale normalization operation is performed, linearly mapping the original image's gray-scale values to the range [0,1] with the mapping formula I_norm = (I−I_min)/(I_max−I_min), where I is the original gray-scale value, and I_min and I_max are the minimum and maximum gray-scale values of a single image. This operation eliminates the gray-scale shift between different scan time steps, unifying the data distribution.

To mitigate the issue of model overfitting, data augmentation strategies are adopted to increase the training samples, including random rotation, horizontal and vertical flipping, and random scaling. The rotation and flipping operations simulate slight pose changes of the core during scanning, while the scaling operation enhances the model's ability to adapt to interface features at different scales. The augmented dataset is split into training, validation, and test sets in a 7:2:1 ratio. The training set is used for model parameter iteration and update, the validation set is used for hyperparameter tuning and overfitting monitoring during the training process, and the test set is used for objective evaluation of the model's final segmentation performance. This splitting ratio follows the conventional setup in medical image segmentation and ensures the reliability of the evaluation results.

2.2 Three-dimensional reconstruction of CT gas-water interface based on multi-feature fusion and transformer matching

2.2.1 Dense feature point extraction using multi-feature fusion

Single feature point extraction methods cannot fully cover the complex characteristics of the gas-water interface. ORB features are fast in real-time but have weak scale adaptability, Harris corner detection is sensitive to edge features but lacks noise resistance, and SURF has scale invariance but high computational complexity. Therefore, this paper combines ORB, Harris, and SURF algorithms to extract dense feature points from the interface, improving feature point integrity, stability, and robustness by complementary advantages, which lays the foundation for subsequent accurate matching.

ORB feature point extraction is based on the FAST algorithm to detect corner points. After selecting an initial threshold T = 20, non-maximum suppression is applied to remove redundant points and ensure a sparse and uniform distribution of feature points. The centroid of the gray-scale region is calculated to determine the main direction, and a 256-dimensional rotation-invariant BRIEF descriptor is generated based on this direction, providing rotational invariance to the feature points. Harris feature point extraction is performed by calculating the second-order moment matrix of the image gray-scale:

$M=\left[\begin{array}{cc}\mathrm{I}_{\mathrm{x}}^2 & \mathrm{I}_{\mathrm{x}} \mathrm{I}_{\mathrm{y}} \\ \mathrm{I}_{\mathrm{x}} \mathrm{I}_{\mathrm{y}} & \mathrm{I}_{\mathrm{y}}^2\end{array}\right]$                        (3)

where, I_x and I_y are the gray-scale gradients in the x and y directions, respectively. The response value is calculated as R = det(M)−k·(trace(M))², where (k = 0.04). Points with a response value greater than R_max×0.01 are selected as corner points. This threshold is experimentally validated and effectively selects key feature points along the interface boundary. SURF feature point extraction detects scale-space extrema using the Hessian matrix, and sub-pixel level feature point localization is achieved through interpolation. The main direction is determined by calculating the Haar wavelet response in the feature point neighborhood, and finally, a 64-dimensional SURF descriptor is generated, ensuring scale invariance.

Multi-feature fusion uses a voting method and weighted fusion strategy to enhance feature quality. The voting method retains the feature point with the highest descriptor matching score at the same pixel location and eliminates redundant information. For non-redundant feature points extracted by the three algorithms, the descriptors are fused using the weighted percentages of ORB (40%), Harris (30%), and SURF (30%). This weight distribution comprehensively considers the advantages of the three features: ORB descriptors have high dimensionality and strong discriminability, Harris corners have accurate positioning, and SURF has good scale adaptability. The weighted fusion descriptors integrate multi-dimensional feature information, significantly improving the robustness of subsequent feature matching.

2.2.2 Feature point accurate matching based on transformer

CT slices of the gas-water interface often have scale differences and gray-scale noise. Traditional feature matching methods are easily affected by scale changes and are difficult to capture the global correlations between feature points, leading to low matching accuracy and high outlier ratios. This paper designs a Transformer-based "from coarse to fine" feature point matching model, which realizes scale-invariant accurate matching through feature encoding, global attention matching, and optimization iterations, effectively solving the above problems.

First, multi-time-step CT gas-water interface slice images are input into a 3D CNN convolutional neural network to extract 3D feature maps. This network strengthens the 3D texture and morphological features of the gas-water interface through multiple convolution operations, providing high-recognition basic features for subsequent matching.

Next, the 3D feature maps are input into a 3D self-attention and cross-attention module. Let the reference time-step feature be S_i and the target time-step feature be S_j. The self-attention mechanism within the module captures the global correlations within the same feature map, and the cross-attention mechanism builds the feature associations between S_i and S_j, calculating the matching relationship S_i/S_j=d_i/d_j. Finally, the CT gas-water interface depth feature map is output, and the coarse matching associations of the feature points are preliminarily determined.

Finally, a feature-point-guided aggregation module optimizes the matching results. This module includes a cross-attention module and a linear cross-attention module. By combining upsampling fusion and downsampling fusion operations, multi-scale feature aggregation is performed on the 3D feature map. The above module structure is repeated four times to fully integrate the interface features at different scales, ultimately generating a CT gas-water interface feature point matching matrix of size H/8×W/8. This process, through multi-scale aggregation and attention correlation, balances local feature consistency and global correlation information, significantly improving the accuracy and stability of the matching results, and effectively reducing the interference of scale differences and noise in the matching process. Figure 3 shows the Transformer modeling method architecture for CT gas-water interface feature point matching.

Figure 3. Transformer modeling method architecture for CT gas-water interface feature point matching

2.2.3 Three-dimensional point cloud reconstruction and optimization

The accurately matched feature point pairs provide reliable two-dimensional correspondences for three-dimensional reconstruction. Based on camera calibration parameters and the triangulation principle, the 3D point cloud of the gas-water interface can be generated. Subsequent denoising and smoothing optimization further enhance the point cloud quality, providing high-precision data support for the subsequent morphological quantification.

The 3D point cloud is generated by calculating the 3D coordinates of the feature points based on the triangulation principle. The camera intrinsic matrix K and extrinsic matrix [R|t] are obtained through preliminary camera calibration, where R is the rotation matrix, and t is the translation vector. For the accurately matched points (u1,v1) and (u2,v2) in the reference and target images, a linear system of equations is constructed:

$\left\{\begin{array}{l}s_1\left[\begin{array}{c}u_1 \\ v_1 \\ 1\end{array}\right]=K[I \mid 0] P \\ s_2\left[\begin{array}{c}u_1 \\ v_1 \\ 1\end{array}\right]=K[R \mid t] P\end{array}\right.$            (4)

where, s₁ and s₂ are scale factors, and P is the 3D point coordinates to be solved. Solving this system of equations gives the 3D coordinates of a single feature point. This operation is repeated for all CT slices across all time steps, and the 3D coordinates of all feature points are integrated to generate the complete 3D point cloud of the gas-water interface, P_cloud∈R^M×3, where M is the total number of points in the point cloud.

Point cloud optimization is performed step-by-step using statistical filtering and moving least squares, balancing denoising and detail preservation. The statistical filtering sets the number of neighboring points (k=20), calculates the mean and standard deviation of the distances between each point and its neighboring points, and removes outlier points with a distance greater than the mean plus 2 times the standard deviation. This parameter setting effectively removes isolated points caused by scan noise and matching errors. The moving least squares method smooths the point cloud while preserving subtle morphological features of the interface. By constructing a locally weighted polynomial surface to fit the point cloud data, local deviations in the point cloud are corrected. The optimized point cloud maintains good smoothness while accurately restoring the true geometric morphology of the gas-water interface, providing high-quality foundational data for subsequent morphological quantification analysis.

2.3 Morphological quantification index system for CT gas-water interface

Based on the high-precision three-dimensional point cloud reconstructed in Section 2.2, a “geometric feature-dynamic feature” dual-dimensional morphological quantification index system is constructed. Geometric features focus on the static spatial morphology of the interface, characterizing the distribution range, occupied volume, and morphological smoothness of the interface. Dynamic features correlate point cloud data from different time steps to reveal the evolutionary patterns of the interface during the displacement process. The two types of indicators complement each other, providing an objective basis for the quantification analysis of the gas-water process, displacement efficiency evaluation, and stability judgment.

2.3.1 Geometric feature indicators

Geometric feature indicators are used to precisely characterize the static spatial morphology of the gas-water interface. Three core indicators are selected: interface surface area, interface enclosing volume, and average curvature. These indicators cover key geometric information from three dimensions: distribution range, gas-phase occupied space, and morphological smoothness. The calculation process is based on the reconstructed three-dimensional point cloud, ensuring the objectivity and accuracy of the indicators.

Interface Surface Area: This is the core indicator for representing the spatial distribution of the interface. The calculation first converts the discrete point cloud into a continuous surface. The Poisson surface reconstruction algorithm is used to fit the three-dimensional point cloud. This algorithm solves the Poisson equation to construct a triangular mesh that fits the topology of the point cloud, effectively preserving the subtle morphological features of the interface. For the generated triangular mesh, the area of each triangle is calculated and summed to obtain the total interface surface area S. This indicator directly reflects the contact range between the gas-water interface and the water phase. The larger the surface area, the broader the gas-water exchange interface, providing basic data for subsequent analysis of mass transfer efficiency. Interface Enclosing Volume: This is used to quantify the spatial size occupied by the gas phase within the core sample. The Axis-Aligned Bounding Box (AABB) method is used to calculate the volume. All three-dimensional coordinates of the point cloud are traversed to extract the maximum and minimum values along the x, y, and z axes, and a rectangular bounding box is constructed based on these extreme values. The volume V of this box is the interface enclosing volume. This method is computationally efficient and reliable in precision, providing a clear reflection of the gas phase's spatial occupancy ability during the displacement process.

Average Curvature: This focuses on the morphological smoothness and convex-concave characteristics of the interface and is a key indicator for revealing the degree of influence of the core pore structure on the interface. The calculation process is based on a triangular mesh: for each vertex in the mesh, the two principal curvatures k₁ and k₂ are calculated based on the geometric relationship with its neighboring vertices. The vertex curvature Hi is the average of the two principal curvatures: H_i=(k₁+k₂)/2. The arithmetic mean of all vertex curvatures gives the average curvature H_avg. The principal curvature is calculated based on quadratic surface fitting to ensure the accuracy of the curvature calculation. The sign of H_avg represents the overall convex or concave shape of the interface: positive values indicate the interface is convex, while negative values indicate the interface is concave. The absolute value reflects the smoothness of the interface: smaller absolute values indicate a smoother interface, while larger values suggest significant fluctuations due to pore throat blockages.

2.3.2 Dynamic feature indicators

Dynamic feature indicators are used to correlate three-dimensional point cloud data from different time steps, characterizing the evolution patterns of the gas-water interface over time. Two indicators are selected: interface migration rate and interface fluctuation amplitude. These indicators quantify the dynamic characteristics of the displacement process from the perspectives of macroscopic advancement efficiency and morphological stability, providing quantitative support for evaluating displacement effects and optimizing displacement parameters.

Interface Migration Rate: This indicator is used to characterize the overall advancing speed of the gas-water interface and is the core dynamic indicator for reflecting displacement efficiency. Before calculation, the overall spatial position of the interface at each time step must be determined. The centroid coordinates are used as the representative position of the interface. For the point cloud at time t, the arithmetic mean of all three-dimensional coordinates of the points is calculated to obtain the centroid G(x,y,z) of the interface at that time. Let the centroids at time t₁ and t₂ be G₁(x₁,y₁,z₁) and G₂(x₂,y₂,z₂), respectively. The interface migration rate v is calculated by the formula:

$v=\frac{\sqrt{(x 2-x 1)^2+(y 2-y 1)^2+(z 2-z 1)^2}}{t 2-t 1}$                      (5)

This indicator, through the ratio of the centroid distance to the time difference, eliminates the interference of local fluctuations on the overall advancement speed and can accurately reflect the macroscopic advancement efficiency of the interface under displacement pressure. The magnitude of this value is directly related to key parameters such as displacement pressure and core permeability.

Interface Fluctuation Amplitude: This indicator is used to evaluate the stability of the interface morphology during the displacement process, focusing on the fluctuation degree of the interface in the displacement direction. The preset displacement direction in the gas-water experiment is the z-axis. The fluctuation of the interface along this direction directly reflects morphological stability. Therefore, the z-direction is selected as the dimension for calculating the fluctuation amplitude. For the point cloud at time t, the coordinates of all points along the z-direction are extracted, and the maximum value z_max and minimum value z_min are determined. The difference between these two values is the interface fluctuation amplitude Δh=z_max-z_min. The physical meaning of this indicator is clear: the smaller Δh is, the smoother the interface in the displacement direction, the more stable the morphology, and the more uniform the displacement process. If Δh increases, it indicates that the interface is influenced by core heterogeneity, pore throat blockages, or other factors, leading to obvious "finger" or local fluctuations and reduced stability. By tracking the changes in Δh at different time steps, the stability evolution of the displacement process can be dynamically monitored.

This paper proposes an integrated method for 3D segmentation, reconstruction, and morphological quantification of CT gas-water interfaces based on an improved 3D Attention U-Net, which has been validated with a real core CT dataset and shown significant effectiveness. The core performance of the proposed method, verified with real core CT gas-water data, is as follows: in segmentation, the Dice coefficient of the improved 3D Attention U-Net reaches 0.92, improving by 5.7% compared to the original model, IoU is 0.86, the average surface distance is 2.1 μm, and the HD is 5.3 μm, which is a reduction of 56.2% and 49.5% respectively compared to the traditional 3D U-Net, achieving pixel-level overlap and sub-micron-level boundary restoration; in morphological restoration, the segmentation results are highly consistent with the measured profiles in the X/Y/Z direction semi-variance curves, with local morphological coincidence in the X direction, inflection point reproduction in the Y direction with base value deviation <10%, and precise matching of displacement features in the Z direction; in 3D reconstruction, the point cloud RMSE is only 0.021mm, with relative deviations of surface area, enclosing volume, and average curvature of 2.3%, 1.8%, and 3.1%, respectively, all meeting experimental precision requirements; in quantitative evolution, the core morphological parameters have a deviation of less than 2% compared to the measured values, which can accurately capture the evolution of the three-stage displacement indicators of "expansion-filling-stabilization," providing high-precision data support for the analysis of the correlation mechanism between "interface morphology" and "displacement efficiency."

This research has certain limitations: the dataset only covers specific core types and displacement conditions, and the generalization to high-noise images and multiphase mixed interface scenarios has not been fully validated. The model processing time for batch data is relatively long, making it difficult to meet real-time analysis requirements. Future research can be deepened in three areas: expanding the dataset to include different lithologies, permeabilities, and displacement pressures to improve method robustness; optimizing processing efficiency through lightweight networks and parallel computing; and integrating pore-scale flow simulation to build a "segmentation-quantification-mechanism" model to further reveal the intrinsic coupling mechanism between interface evolution and displacement efficiency.