简介这份资源面向计算机视觉入门者、地质图像分析方向的学生以及需要完成期末大作业或课程设计的学习者提供了一套基于Python的CT岩芯与岩石裂缝语义分割完整方案。包内共15个文件以py脚本、jpg示例图像、zbak备份文件及md说明为主压缩包约1.15MB涵盖数据增强、均值计算等处理脚本与岩石、混凝土、CT三类样本及对应标注图便于直接运行与对照实验。资源围绕像素级分类任务展开涉及Pillow、OpenCV等基础图像处理库以及TensorFlow、PyTorch配合Keras搭建分割网络的常见技术路线可用于断层扫描图像的裂隙分布量化分析。目前已有71人学习下载适合希望快速获取可复现代码、理解端到端分割流程并积累地质工程场景实践经验的读者参考。1. 从一张 CT 岩心切片说起为什么通用语义分割模型在岩石裂缝上会翻车拿一块直径 5 厘米的碳酸盐岩岩心去做 CT 扫描重建出来的体数据里裂缝往往只有 1 到 3 个体素宽灰度值和周围基质、方解石充填带的差异可能不到 8%。你把这张切片丢给在 Cityscapes 或 VOC 上训过的 DeepLabV3它大概率会把整条裂缝判成背景或者把高密度的矿物颗粒误判成裂缝。这不是模型不行是数据域差得太远——自然图像里物体有纹理、有边界、有语义上下文而 CT 岩心切片里裂缝就是一条灰度略低的细线上下文几乎为零。基于 Python 的 CT 岩心与岩石裂缝语义分割系统要解决的就是这件事把岩心 CT 序列切成二维切片用像素级分类把裂缝、基质、高密度矿物、孔隙这几类分开输出可直接量化的裂缝面积占比、开度分布和连通性指标。它适合做岩石力学、数字岩心、油气储层表征的工程师和研究生也适合已经会写 Python、想找一个真实工业数据集练语义分割的开发者。整套东西的核心不是模型多深而是数据怎么标、类别怎么定、损失函数怎么压住类别不平衡——这三点决定了你是拿到 0.6 的 mIoU 还是 0.85。2. 数据准备与标注从 CT 序列到可训练的裂缝掩码2.1 岩心 CT 数据的读取与切片策略CT 岩心数据通常是 16 位无符号整数的体数据格式可能是 DICOM 序列、TIFF 堆栈或 RAW 二进制。第一步不是急着写 Dataset而是先确认体数据的维度顺序和灰度范围。常见做法是用SimpleITK或pydicom读入统一转成numpy的(D, H, W)数组再沿轴向切片。import SimpleITK as sitk import numpy as np import os def load_ct_volume(series_dir): 读取 DICOM 序列返回 (D, H, W) 的 float32 体数据 series_dir: DICOM 文件所在目录 reader sitk.ImageSeriesReader() dicom_names reader.GetGDCMSeriesFileNames(series_dir) reader.SetFileNames(dicom_names) volume reader.Execute() arr sitk.GetArrayFromImage(volume).astype(np.float32) # (D, H, W) # CT 值截断到岩石有效范围去掉空气和金属伪影 arr np.clip(arr, -1000, 3000) # 归一化到 [0, 1] arr (arr - arr.min()) / (arr.max() - arr.min() 1e-8) return arr def slice_along_axis(volume, axis0, stride1): 沿指定轴切片stride 控制采样密度 slices [] for i in range(0, volume.shape[axis], stride): if axis 0: slices.append(volume[i, :, :]) elif axis 1: slices.append(volume[:, i, :]) else: slices.append(volume[:, :, i]) return np.stack(slices, axis0)逻辑说明load_ct_volume先做灰度截断再归一化这一步很关键。CT 值里空气约 -1000 HU金属伪影可能到 3000 HU 以上不截断的话归一化会被极端值拉偏裂缝的微弱对比度直接消失。slice_along_axis的stride参数控制采样密度岩心轴向分辨率通常 0.1 到 0.5 毫米如果原始层间距是 0.3 毫米stride1就够如果层间距远小于像素尺寸可以stride2或3降冗余。参数上np.clip的上限我一般设 3000下限 -1000这是岩石 CT 的常规窗口。如果你的岩心含黄铁矿等高密度矿物上限可以放到 4000但要注意归一化后裂缝和矿物的对比度是否还够。2.2 裂缝标注的类别定义与掩码生成语义分割的类别定义直接决定项目能不能落地。岩心 CT 里我一般分四类0 背景基质、1 裂缝、2 高密度矿物、3 孔隙。裂缝的标注最麻烦因为细裂缝在切片上只有几个像素宽手工勾画效率极低。常见做法是先用阈值分割加形态学操作生成粗掩码再人工修正。阈值不是拍脑袋定的要看灰度直方图的双峰分布。import cv2 import numpy as np def generate_crack_mask(slice_2d, low_thresh0.35, high_thresh0.55): 基于灰度阈值的裂缝粗掩码生成 slice_2d: 归一化后的二维切片 low_thresh: 裂缝灰度上限 high_thresh: 基质灰度下限 # 双阈值分割中间区域用形态学闭运算连接 crack_candidate (slice_2d low_thresh).astype(np.uint8) matrix_region (slice_2d high_thresh).astype(np.uint8) # 去除小噪点 kernel cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3, 3)) crack_candidate cv2.morphologyEx(crack_candidate, cv2.MORPH_OPEN, kernel) crack_candidate cv2.morphologyEx(crack_candidate, cv2.MORPH_CLOSE, kernel) # 只保留细长区域去掉圆形孔洞 contours, _ cv2.findContours(crack_candidate, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) mask np.zeros_like(slice_2d, dtypenp.uint8) for cnt in contours: area cv2.contourArea(cnt) if area 5: # 去掉太小噪点 continue x, y, w, h cv2.boundingRect(cnt) aspect_ratio max(w, h) / (min(w, h) 1e-6) if aspect_ratio 3: # 细长区域才可能是裂缝 cv2.drawContours(mask, [cnt], -1, 1, -1) return mask逻辑说明low_thresh和high_thresh之间是过渡带不直接参与分割避免把裂缝边缘的模糊像素误判。aspect_ratio 3这个条件是我踩过坑之后加的——孔隙在切片上也是低灰度圆形区域不加形状过滤会把孔隙全标成裂缝。area 5去掉的是 CT 重建噪声不是真实裂缝。参数上low_thresh和high_thresh需要根据你的岩心灰度直方图调。我一般先用matplotlib画直方图找裂缝峰和基质峰之间的谷底谷底位置就是low_thresh谷底往右 0.1 到 0.15 是high_thresh。如果你的岩心裂缝充填了方解石灰度可能比基质还高这时候阈值逻辑要反过来先做矿物识别再做裂缝。提示粗掩码只是给标注人员打底不能直接当标签用。细裂缝的连通性和开度必须人工修正否则训练出来的模型会把断续的裂缝段连成一条量化结果完全不可信。2.3 数据集划分与增强策略岩心 CT 切片有个特点相邻切片高度相似如果随机划分训练集和验证集验证集里的切片可能在训练集里有几乎一样的副本mIoU 虚高。正确做法是按深度区间划分比如前 70% 深度做训练中间 15% 做验证最后 15% 做测试。def split_by_depth(num_slices, train_ratio0.7, val_ratio0.15): 按深度顺序划分避免相邻切片泄漏 train_end int(num_slices * train_ratio) val_end int(num_slices * (train_ratio val_ratio)) train_idx list(range(0, train_end)) val_idx list(range(train_end, val_end)) test_idx list(range(val_end, num_slices)) return train_idx, val_idx, test_idx增强策略上岩心切片不能随便用旋转和翻转。轴向切片做水平翻转是合理的因为岩石结构没有方向性但垂直翻转会改变层理方向对沉积岩来说不物理。我一般只用水平翻转、小角度旋转±10 度和灰度抖动±5%不做弹性形变——弹性形变会把裂缝拉成不真实的形状。3. 模型选型与训练U-Net、DeepLabV3 和 SegFormer 在裂缝上的真实表现3.1 为什么裂缝分割更偏向 U-Net 系裂缝是细长结构语义分割模型里对细结构最友好的是 U-Net 的跳跃连接。编码器下采样会丢失细裂缝的位置信息跳跃连接把浅层高分辨率特征直接送到解码器这是 U-Net 在裂缝上比 DeepLabV3 稳的原因。DeepLabV3 的空洞卷积虽然扩大了感受野但对 1 到 2 像素宽的裂缝空洞卷积的采样点可能直接跳过裂缝像素。SegFormer 是 Transformer 系全局注意力对裂缝的连通性建模有优势但需要更多数据。岩心 CT 标注成本高通常只有几百到几千张标注切片SegFormer 容易过拟合。我的经验是标注数据少于 2000 张U-Net 加残差编码器如 ResNet34最稳超过 5000 张可以试 SegFormer 的 MiT-B2 backbone。import torch import torch.nn as nn from torchvision.models import resnet34 class DecoderBlock(nn.Module): def __init__(self, in_ch, skip_ch, out_ch): super().__init__() self.conv1 nn.Conv2d(in_ch skip_ch, out_ch, 3, padding1) self.bn1 nn.BatchNorm2d(out_ch) self.conv2 nn.Conv2d(out_ch, out_ch, 3, padding1) self.bn2 nn.BatchNorm2d(out_ch) self.relu nn.ReLU(inplaceTrue) def forward(self, x, skip): x torch.cat([x, skip], dim1) x self.relu(self.bn1(self.conv1(x))) x self.relu(self.bn2(self.conv2(x))) return x class ResUNet(nn.Module): def __init__(self, num_classes4): super().__init__() backbone resnet34(pretrainedTrue) self.encoder0 nn.Sequential(backbone.conv1, backbone.bn1, backbone.relu) # 1/2 self.encoder1 nn.Sequential(backbone.maxpool, backbone.layer1) # 1/4 self.encoder2 backbone.layer2 # 1/8 self.encoder3 backbone.layer3 # 1/16 self.encoder4 backbone.layer4 # 1/32 self.decoder4 DecoderBlock(512, 256, 256) self.decoder3 DecoderBlock(256, 128, 128) self.decoder2 DecoderBlock(128, 64, 64) self.decoder1 DecoderBlock(64, 64, 64) self.final nn.Conv2d(64, num_classes, 1) def forward(self, x): e0 self.encoder0(x) e1 self.encoder1(e0) e2 self.encoder2(e1) e3 self.encoder3(e2) e4 self.encoder4(e3) d4 self.decoder4(e4, e3) d3 self.decoder3(d4, e2) d2 self.decoder2(d3, e1) d1 self.decoder1(d2, e0) return self.final(d1)逻辑说明编码器用 ResNet34 的预训练权重岩心 CT 数据量小从头训收敛慢且容易过拟合。解码器每层先拼接跳跃连接再卷积DecoderBlock里两个 3x3 卷积加 BN 和 ReLU这是 U-Net 的标准配置。最后 1x1 卷积输出num_classes通道不做 softmax因为训练时用CrossEntropyLoss内部会做。参数上num_classes4对应背景、裂缝、矿物、孔隙。如果你的岩心只有裂缝和基质两类改成 2。预训练权重用resnet34(pretrainedTrue)如果下载不了可以换成weightsNone从头训但学习率要降到 1e-4 以下。3.2 损失函数Dice Focal 怎么压住裂缝的类别不平衡裂缝像素在整张切片里占比通常不到 3%用纯交叉熵训练模型会倾向于全预测背景准确率看着有 97%但裂缝的 IoU 接近 0。常见做法是 Dice Loss 和 Focal Loss 加权组合。class DiceFocalLoss(nn.Module): def __init__(self, alpha0.5, gamma2.0, dice_weight0.7): super().__init__() self.alpha alpha self.gamma gamma self.dice_weight dice_weight def forward(self, logits, targets): # Focal Loss 部分 ce nn.functional.cross_entropy(logits, targets, reductionnone) pt torch.exp(-ce) focal self.alpha * (1 - pt) ** self.gamma * ce focal focal.mean() # Dice Loss 部分只对裂缝类计算 probs torch.softmax(logits, dim1) crack_prob probs[:, 1, :, :] # 裂缝类 crack_target (targets 1).float() intersection (crack_prob * crack_target).sum() dice 1 - (2 * intersection 1e-6) / (crack_prob.sum() crack_target.sum() 1e-6) return self.dice_weight * dice (1 - self.dice_weight) * focal逻辑说明Focal Loss 的gamma2.0让模型更关注难分类的裂缝像素alpha0.5平衡正负样本。Dice Loss 只对裂缝类算因为背景和矿物的 Dice 本来就高加进去会稀释裂缝的梯度。dice_weight0.7是我在岩心数据上试出来的Dice 占主导Focal 做辅助。参数上如果你的裂缝更细、占比更低dice_weight可以提到 0.8gamma提到 3.0。但gamma太高会让训练不稳定损失值震荡这时候把学习率降到 1e-4 再试。3.3 训练配置与显存优化岩心 CT 切片尺寸通常是 512x512 或 1024x10241024 的切片在 8GB 显存上跑 ResUNet 会 OOM。常见做法是随机裁剪到 512x512 训练推理时用滑窗。from torch.utils.data import Dataset, DataLoader import random class CoreDataset(Dataset): def __init__(self, slices, masks, crop_size512, augmentTrue): self.slices slices self.masks masks self.crop_size crop_size self.augment augment def __len__(self): return len(self.slices) def __getitem__(self, idx): img self.slices[idx] mask self.masks[idx] h, w img.shape if self.augment: # 随机裁剪 top random.randint(0, h - self.crop_size) left random.randint(0, w - self.crop_size) img img[top:topself.crop_size, left:leftself.crop_size] mask mask[top:topself.crop_size, left:leftself.crop_size] # 水平翻转 if random.random() 0.5: img np.fliplr(img).copy() mask np.fliplr(mask).copy() # 灰度抖动 img img * random.uniform(0.95, 1.05) return torch.from_numpy(img).float().unsqueeze(0), torch.from_numpy(mask).long() # 训练循环关键参数 optimizer torch.optim.AdamW(model.parameters(), lr1e-4, weight_decay1e-5) scheduler torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max50)逻辑说明crop_size512是显存和感受野的折中。np.fliplr后加.copy()是因为 numpy 翻转返回的是视图不 copy 的话后续灰度抖动会改到原数据。AdamW的weight_decay1e-5比 SGD 的动量更适合小数据集。CosineAnnealingLR的T_max50对应 50 个 epoch学习率从 1e-4 余弦降到 0。参数上batch size 在 8GB 显存上设 4 到 6用accumulation_steps2模拟更大 batch。如果验证集 mIoU 在 20 个 epoch 后还在涨把T_max加到 80。4. 推理、后处理与裂缝量化从掩码到开度分布4.1 滑窗推理与大图拼接训练时裁剪到 512推理时整张 1024 切片要滑窗。滑窗重叠 128 像素重叠区域取概率平均避免拼接缝。def sliding_window_inference(model, image, window_size512, overlap128): 对单张大图做滑窗推理返回类别概率图 model.eval() h, w image.shape stride window_size - overlap prob_map np.zeros((4, h, w), dtypenp.float32) count_map np.zeros((h, w), dtypenp.float32) with torch.no_grad(): for y in range(0, h - window_size 1, stride): for x in range(0, w - window_size 1, stride): patch image[y:ywindow_size, x:xwindow_size] tensor torch.from_numpy(patch).float().unsqueeze(0).unsqueeze(0) logits model(tensor) probs torch.softmax(logits, dim1).squeeze(0).numpy() prob_map[:, y:ywindow_size, x:xwindow_size] probs count_map[y:ywindow_size, x:xwindow_size] 1 count_map np.maximum(count_map, 1e-6) prob_map / count_map return prob_map逻辑说明stride window_size - overlap重叠区域概率累加后除以计数得到平均概率。count_map防止除零。最后argmax得到类别掩码。参数上overlap128是 512 窗口的 25%再小会有拼接缝再大推理时间线性增加。如果显存够window_size可以设 768重叠 192。4.2 裂缝开度与面积占比的量化拿到裂缝掩码后量化指标才是岩心分析真正要的东西。开度用距离变换算面积占比直接统计像素。from scipy.ndimage import distance_transform_edt def quantify_crack(mask, pixel_size_mm0.05): mask: 二值裂缝掩码 pixel_size_mm: 单个像素的物理尺寸毫米 crack_pixels (mask 1) area_ratio crack_pixels.sum() / mask.size # 开度距离变换的局部最大值乘 2 dist distance_transform_edt(crack_pixels) # 骨架上的距离值代表到最近边界的距离乘 2 是开度 from skimage.morphology import skeletonize skeleton skeletonize(crack_pixels) openings dist[skeleton] * 2 * pixel_size_mm return { area_ratio: area_ratio, mean_opening_mm: openings.mean() if len(openings) 0 else 0, max_opening_mm: openings.max() if len(openings) 0 else 0, crack_length_mm: skeleton.sum() * pixel_size_mm }逻辑说明distance_transform_edt算每个裂缝像素到最近背景的距离骨架上的距离值就是裂缝半开度乘 2 得全开度。pixel_size_mm必须从 CT 扫描参数里拿不能瞎设否则开度数值没有物理意义。参数上pixel_size_mm常见 0.02 到 0.1 毫米取决于扫描分辨率。如果你的 CT 是 1024x1024 覆盖 50 毫米直径像素尺寸约 0.05 毫米。5. 避坑与排查岩心裂缝分割里最容易翻车的 5 个地方5.1 现象验证集 mIoU 0.85测试集掉到 0.5原因按切片随机划分数据集相邻切片泄漏。岩心 CT 相邻切片差异极小训练集里见过几乎一样的图验证集自然高。解决改成按深度区间划分训练、验证、测试在深度上完全不重叠。如果数据量够按岩心样本划分不同岩心的切片不混。5.2 现象裂缝预测断断续续连通性差原因损失函数只用了交叉熵模型对细裂缝的梯度不够。或者后处理直接 argmax没有做连通性修复。解决损失函数加 Dice Loss权重 0.7 以上。后处理用形态学闭运算连接断点但 kernel 不要超过 3x3否则会把两条独立裂缝连成一条。5.3 现象开度量化结果偏大 2 到 3 倍原因pixel_size_mm设错了或者距离变换后没有乘 2。距离变换得到的是半开度直接当开度用会小一半如果像素尺寸用了 CT 重建的体素尺寸而不是切片像素尺寸又会偏大。解决确认 CT 扫描的PixelSpacing标签单位是毫米。距离变换结果乘 2 再乘像素尺寸。5.4 现象训练 loss 震荡mIoU 不收敛原因学习率太高或者gamma设太大导致 Focal Loss 梯度爆炸。解决学习率从 1e-4 降到 5e-5gamma从 3.0 降到 2.0。加梯度裁剪torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm1.0)。5.5 现象推理时显存够但速度极慢原因滑窗推理没有用torch.no_grad()或者 batch 维度没加导致逐像素推理。解决推理包在with torch.no_grad():里patch 加 batch 和 channel 维度。如果还慢把模型转 ONNX 或 TensorRT512 窗口的 ResUNet 在 RTX 3060 上单张 1024 切片约 1.5 秒。6. 进阶技巧用测试时增强和伪标签把 mIoU 再推 3 个点标注数据不够的时候测试时增强TTA是最便宜的涨点手段。对同一张切片做水平翻转、小角度旋转分别推理后把概率图平均再 argmax。裂缝分割里水平翻转和小角度旋转是物理合理的垂直翻转不行。def tta_inference(model, image): 测试时增强原图 水平翻转 旋转 ±10 度 probs [] # 原图 probs.append(sliding_window_inference(model, image)) # 水平翻转 flipped np.fliplr(image).copy() p_flip sliding_window_inference(model, flipped) probs.append(np.fliplr(p_flip).copy()) # 旋转 10 度 from scipy.ndimage import rotate rotated rotate(image, 10, reshapeFalse, order1) p_rot sliding_window_inference(model, rotated) probs.append(rotate(p_rot, -10, reshapeFalse, order1)) # 平均 return np.mean(probs, axis0)逻辑说明每个增强版本推理后要逆变换回原图坐标系再平均。rotate的order1是双线性插值order0是最近邻概率图用双线性更平滑。TTA 的代价是推理时间乘 3但 mIoU 通常能涨 2 到 3 个点。伪标签是另一条路先用训练好的模型对未标注切片推理置信度高于 0.9 的像素当伪标签加入训练集重新训。岩心 CT 里裂缝的伪标签要人工抽检因为模型容易把矿物边缘误判成裂缝。我一般只对裂缝类做伪标签背景和矿物类不碰。还有一个技巧是深监督在 U-Net 解码器的每个尺度加辅助损失让浅层也直接学裂缝特征。辅助损失权重从 0.1 到 0.4 逐层增加深层权重最大。这个在裂缝细、下采样丢失严重的时候特别有用mIoU 能再涨 1 到 2 个点。我自己的习惯是每训完一个模型先拿 20 张测试切片做 TTA 和不开 TTA 的对比如果 TTA 涨点不到 1 个点说明模型本身已经够稳不用加这个推理开销。伪标签只在标注数据少于 500 张的时候用多了反而引入噪声。裂缝分割这件事数据质量比模型结构重要得多标注的时候多花一小时修正细裂缝比调一天超参管用。希望帮到你。本文还有配套的精品资源点击获取