ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

基于Python的CT岩心裂缝语义分割:从数据预处理到定量分析全流程

基于Python的CT岩心裂缝语义分割:从数据预处理到定量分析全流程 简介这份资源面向计算机视觉入门者、地质图像分析方向的学生以及需要完成期末大作业或课程设计的学习者提供一套基于Python的CT岩芯与岩石裂缝语义分割完整方案用于解决像素级裂隙识别与量化分析问题。压缩包共15个文件约1.15MB包含3个py脚本用于数据增强与均值计算6张jpg图像作为岩石、混凝土及CT扫描的原始图与标注掩码另有zbak备份、md说明文档和gitignore等辅助文件结构紧凑便于直接运行与二次修改。资源覆盖从图像读取、形态学处理到深度学习模型搭建的端到端流程可帮助读者理解语义分割在地质非破坏性检测中的落地方式掌握裂隙分布模式识别与定量分析思路。目前已有72人学习下载适合作为课程实践与算法入门的参考素材。1. CT岩心裂缝语义分割从灰度切片到可训练掩码的完整链路拿到一批工业 CT 扫描的岩心切片想把里面的裂缝自动勾出来这件事在石油地质和岩石力学圈子里需求很实在。人工在几百张切片上逐像素描裂缝一张图少则十几分钟多则半小时标注一致性还差。基于 Python 的 CT 岩心与岩石裂缝语义分割系统本质就是把「切片预处理 → 裂缝像素级标注 → 语义分割模型训练 → 推理出掩码」这条链路用代码串起来让裂缝从灰度图里被自动分离成二值掩码。它适合两类人一类是手里有 CT 数据、想快速搭一套能跑的分割流程的地质或岩土工程师另一类是刚接触语义分割、想找一个真实工业场景练手的 Python 开发者。这套东西不追求 SOTA 指标追求的是数据能进、模型能训、掩码能出、结果能复核。下面按我实际搭过的顺序把每个环节的参数和坑讲清楚。2. 数据准备CT 切片怎么变成语义分割能吃的格式CT 岩心数据通常是 16 位灰度 TIFF 或 DICOM 序列灰度范围能到 0–65535而裂缝在灰度上往往只比基质暗一点点对比度极低。直接丢给模型网络学不到东西。所以第一步不是写模型是把数据整理成「图像 掩码」成对、灰度归一、尺寸统一的格式。语义分割和实例分割的区别在这里很关键裂缝分割只关心「这个像素是不是裂缝」不关心「这是第几条裂缝」所以用语义分割的标签体系一张图对应一张单通道掩码像素值 0 是背景、1 是裂缝不需要给每条裂缝编号。这也是为什么标题里写的是语义分割而不是实例分割——岩心裂缝的工程诉求是统计裂缝面积占比和走向不是数裂缝条数。2.1 从 CT 原始切片到 PNG 的批量转换工业 CT 导出的切片常见是 16 位无符号整型Python 里用tifffile或pydicom读。转成 8 位 PNG 是为了后续标注工具和大部分分割框架好处理但直接线性压缩会丢掉裂缝和基质的微小灰度差所以要先做对比度拉伸。下面这段是我常用的转换脚本。import numpy as np import tifffile from PIL import Image import os def ct_slice_to_png(src_dir, dst_dir, low_pct1, high_pct99): os.makedirs(dst_dir, exist_okTrue) for name in sorted(os.listdir(src_dir)): if not name.lower().endswith((.tif, .tiff)): continue img tifffile.imread(os.path.join(src_dir, name)).astype(np.float32) # 按百分位裁剪避免个别极亮/极暗像素拉垮整体对比度 lo, hi np.percentile(img, [low_pct, high_pct]) img np.clip(img, lo, hi) img (img - lo) / (hi - lo 1e-6) * 255.0 img img.astype(np.uint8) Image.fromarray(img).save(os.path.join(dst_dir, name.rsplit(., 1)[0] .png)) ct_slice_to_png(./ct_raw, ./ct_png)逻辑说明np.percentile取 1% 和 99% 分位作为拉伸上下限比直接用 min/max 稳因为 CT 里常有金属矿物或扫描伪影造成的极端亮斑。参数low_pct和high_pct是可调的裂缝对比度还是不够时可以收到 5/95让拉伸更激进。转换后务必抽查几张确认裂缝没有在拉伸中被压没。2.2 掩码标注规范与目录结构标注工具用 LabelMe 或 CVAT 都行导出成 PNG 掩码。关键是标注规范要统一裂缝边缘怎么算、宽度小于 2 像素的细缝标不标、裂缝和孔洞连在一起时怎么切分。我一般定三条规则宽度小于 2 像素的裂缝不标模型学不稳标注也不一致裂缝与孔洞连通时只标裂缝主体孔洞归背景掩码边缘允许 1 像素误差。目录结构按下面组织训练脚本直接按文件名配对。dataset/ images/ slice_0001.png slice_0002.png masks/ slice_0001.png slice_0002.png掩码必须是单通道、像素值只有 0 和 1或 0 和 255如果标注工具导出的是彩色索引图要转成灰度二值。这一步不做训练时 loss 会直接报类别数不匹配。2.3 数据集划分与增强的边界按切片划分训练/验证/测试比例 7:2:1。注意不能随机打乱所有切片再分因为相邻切片高度相似随机分会导致验证集里出现和训练集几乎一样的图指标虚高。正确做法是按深度区间分块比如前 70% 深度做训练中间 20% 验证最后 10% 测试。增强用水平翻转、垂直翻转、±10° 旋转、亮度抖动就够了。CT 岩心是各向异性的过度旋转会破坏裂缝的真实走向分布所以旋转角度别开太大。3. 语义分割模型选型与训练U-Net 为什么还是裂缝分割的稳妥起点裂缝分割是典型的「细长目标 低对比度」任务目标像素占比往往不到 5%属于强类别不平衡。这个场景下U-Net 系列依然是性价比最高的起点编码器-解码器加跳跃连接的结构对细结构的保留能力好训练数据需求也比 Transformer 类模型低。如果数据量上千张、算力充足可以上 SegFormer 或 DeepLabV3但对大多数岩心项目U-Net 加一个预训练 ResNet 编码器就够用。下面给一套能直接跑的 PyTorch 训练代码。3.1 U-Net 模型定义与损失函数import torch import torch.nn as nn import torchvision class UNet(nn.Module): def __init__(self, pretrainedTrue): super().__init__() resnet torchvision.models.resnet34(weightsIMAGENET1K_V1 if pretrained else None) self.enc0 nn.Sequential(resnet.conv1, resnet.bn1, resnet.relu) # 1/2 self.enc1 nn.Sequential(resnet.maxpool, resnet.layer1) # 1/4 self.enc2 resnet.layer2 # 1/8 self.enc3 resnet.layer3 # 1/16 self.enc4 resnet.layer4 # 1/32 self.up4 nn.ConvTranspose2d(512, 256, 2, stride2) self.dec4 nn.Sequential(nn.Conv2d(512, 256, 3, padding1), nn.BatchNorm2d(256), nn.ReLU()) self.up3 nn.ConvTranspose2d(256, 128, 2, stride2) self.dec3 nn.Sequential(nn.Conv2d(256, 128, 3, padding1), nn.BatchNorm2d(128), nn.ReLU()) self.up2 nn.ConvTranspose2d(128, 64, 2, stride2) self.dec2 nn.Sequential(nn.Conv2d(128, 64, 3, padding1), nn.BatchNorm2d(64), nn.ReLU()) self.up1 nn.ConvTranspose2d(64, 32, 2, stride2) self.dec1 nn.Sequential(nn.Conv2d(64, 32, 3, padding1), nn.BatchNorm2d(32), nn.ReLU()) self.head nn.Conv2d(32, 1, 1) def forward(self, x): e0 self.enc0(x) e1 self.enc1(e0) e2 self.enc2(e1) e3 self.enc3(e2) e4 self.enc4(e3) d4 self.dec4(torch.cat([self.up4(e4), e3], dim1)) d3 self.dec3(torch.cat([self.up3(d4), e2], dim1)) d2 self.dec2(torch.cat([self.up2(d3), e1], dim1)) d1 self.dec1(torch.cat([self.up1(d2), e0], dim1)) return self.head(d1)逻辑说明编码器用 ResNet34 预训练权重跳跃连接把编码器各层特征拼到解码器对应层这是 U-Net 保留细裂缝的关键。head输出单通道 logits配合下面的损失函数。参数上pretrainedTrue在数据少于 500 张时强烈建议开能明显加快收敛。损失函数用 Dice BCE 组合Dice 负责应对类别不平衡BCE 负责像素级稳定梯度。class DiceBCELoss(nn.Module): def __init__(self, dice_weight0.5): super().__init__() self.dice_weight dice_weight self.bce nn.BCEWithLogitsLoss() def forward(self, logits, targets): bce self.bce(logits, targets) probs torch.sigmoid(logits) intersection (probs * targets).sum() dice 1 - (2 * intersection 1e-6) / (probs.sum() targets.sum() 1e-6) return self.dice_weight * dice (1 - self.dice_weight) * bcedice_weight默认 0.5如果裂缝像素占比低于 2%可以提到 0.7让模型更关注裂缝。但别设成 1.0纯 Dice 在训练初期梯度不稳容易震荡。3.2 训练循环与关键超参from torch.utils.data import Dataset, DataLoader from PIL import Image import numpy as np class CrackDataset(Dataset): def __init__(self, img_dir, mask_dir, size512): self.img_dir, self.mask_dir, self.size img_dir, mask_dir, size self.names sorted(os.listdir(img_dir)) def __len__(self): return len(self.names) def __getitem__(self, idx): name self.names[idx] img Image.open(os.path.join(self.img_dir, name)).convert(L).resize((self.size, self.size)) mask Image.open(os.path.join(self.mask_dir, name)).convert(L).resize((self.size, self.size)) img np.array(img, dtypenp.float32) / 255.0 mask (np.array(mask) 127).astype(np.float32) return torch.from_numpy(img)[None], torch.from_numpy(mask)[None] train_ds CrackDataset(./dataset/images, ./dataset/masks) loader DataLoader(train_ds, batch_size4, shuffleTrue, num_workers2) model UNet().cuda() criterion DiceBCELoss(dice_weight0.6) optimizer torch.optim.AdamW(model.parameters(), lr1e-4, weight_decay1e-4) scheduler torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max50) for epoch in range(50): model.train() for img, mask in loader: img, mask img.cuda(), mask.cuda() optimizer.zero_grad() loss criterion(model(img), mask) loss.backward() optimizer.step() scheduler.step() print(fepoch {epoch}, loss {loss.item():.4f})逻辑说明输入统一 resize 到 512×512convert(L)保证单通道。掩码用127二值化兼容 0/255 和 0/1 两种导出。优化器用 AdamW学习率 1e-4配合余弦退火。batch_size 受显存限制4 是 8GB 显存的稳妥值显存够可以上 8。训练轮数 50 是起点看验证集 Dice 曲线决定要不要加。3.3 评估指标别只看准确率裂缝像素占比低准确率accuracy会骗人——全预测背景也能到 95% 以上。必须看 Dice 和 IoU。Dice 对裂缝这种小目标更敏感IoU 更严格。验证时按切片算 Dice 再平均不要把所有像素混在一起算否则大图会主导指标。def dice_score(pred, target, eps1e-6): pred (torch.sigmoid(pred) 0.5).float() inter (pred * target).sum() return (2 * inter eps) / (pred.sum() target.sum() eps)阈值 0.5 是默认实际部署时可以调后面第 5 章会讲怎么调。4. 推理与后处理让掩码从「能出」到「能用」模型训完推理出的原始掩码往往有毛刺、断线、小噪点。直接拿去统计裂缝面积误差会很大。这一章讲推理脚本和三种后处理把掩码修到能进分析流程。4.1 批量推理脚本import torch from PIL import Image import numpy as np import os def predict(model, img_path, size512, threshold0.5): model.eval() img Image.open(img_path).convert(L).resize((size, size)) x torch.from_numpy(np.array(img, dtypenp.float32) / 255.0)[None, None].cuda() with torch.no_grad(): logits model(x) prob torch.sigmoid(logits)[0, 0].cpu().numpy() mask (prob threshold).astype(np.uint8) * 255 return mask, prob model UNet().cuda() model.load_state_dict(torch.load(best_unet.pth)) for name in os.listdir(./test_images): mask, prob predict(model, os.path.join(./test_images, name)) Image.fromarray(mask).save(f./pred_masks/{name})逻辑说明推理时保持和训练一致的 resize 尺寸否则尺度不匹配会掉点。threshold是二值化阈值默认 0.5后面会讲怎么调。保存概率图prob是为了后处理时能重新选阈值不用重跑模型。4.2 三种后处理去噪、断线连接、阈值调优去噪用连通域面积过滤小于 30 像素的连通域直接删掉这些多半是伪影。from scipy import ndimage def remove_small(mask, min_area30): labeled, n ndimage.label(mask 0) for i in range(1, n 1): if (labeled i).sum() min_area: mask[labeled i] 0 return mask断线连接用形态学闭运算3×3 或 5×5 核把裂缝断开的地方接上。核别开太大否则会把两条平行细缝粘成一条。from scipy.ndimage import binary_closing def connect_cracks(mask, kernel_size3): structure np.ones((kernel_size, kernel_size), dtypebool) return binary_closing(mask 0, structurestructure).astype(np.uint8) * 255阈值调优在验证集上扫 0.3 到 0.7每 0.05 一档算 Dice选最高的。裂缝分割里 0.4 左右往往比 0.5 好因为模型对裂缝边缘的置信度偏低降阈值能把边缘捞回来代价是噪点变多配合面积过滤正好。4.3 后处理顺序与参数联动顺序是先阈值二值化 → 闭运算连断线 → 面积过滤去噪。顺序反了会出问题先面积过滤再闭运算可能把刚连上的细缝又当成小连通域删掉。参数联动上闭运算核越大面积过滤的min_area也要相应调大否则连出来的大块伪影删不掉。我一般固定闭运算核 3×3min_area在 20–50 之间试。5. 避坑与排查裂缝分割里最容易翻车的五件事这一章是我踩过的坑按「现象 → 原因 → 解决」写每条都能对上前面章节的操作。现象一训练 loss 一直不降Dice 卡在 0.1 左右。原因多半是掩码没二值化或者掩码和图像文件名没对上模型在学噪声。解决写个检查脚本随机抽 10 对图把掩码叠加到原图上可视化确认裂缝位置对得上再确认掩码像素值只有 0 和 1或 0 和 255。现象二验证集 Dice 很高测试集一塌糊涂。原因是数据集按随机划分相邻切片泄漏。解决改成按深度区间分块划分训练/验证/测试的切片在深度上不重叠。这个坑最隐蔽指标虚高会让人误以为模型能用。现象三推理掩码全是噪点裂缝反而没出来。原因通常是推理时的 resize 尺寸和训练不一致或者归一化方式不同训练除了 255推理没除。解决把预处理封装成一个函数训练和推理共用杜绝两套代码。现象四细裂缝断成一段一段统计面积偏小。原因是模型对细结构召回不足加上阈值 0.5 偏高。解决降阈值到 0.4加闭运算连断线再面积过滤。如果还断考虑在损失里提高 Dice 权重或在训练时加细裂缝的过采样。现象五换一批新岩心的 CT 数据模型直接失效。原因是不同扫描设备的灰度分布、分辨率、伪影特征都不一样模型过拟合了旧数据。解决新数据上做少量标注几十张用低学习率 1e-5 微调别从头训。如果新数据灰度分布差异大重新做百分位拉伸再微调。6. 把分割结果接进裂缝定量分析一个可复现的统计脚本模型出掩码只是中间产物地质上真正要的是裂缝面积占比、裂缝密度、走向分布这些量。这一章给一个从掩码算裂缝面积占比和等效宽度的脚本并讲怎么验证统计结果可信。import numpy as np from PIL import Image from scipy import ndimage def crack_stats(mask_path, pixel_size_mm0.05): mask np.array(Image.open(mask_path).convert(L)) 127 labeled, n ndimage.label(mask) total_pixels mask.size crack_pixels mask.sum() area_ratio crack_pixels / total_pixels # 每条裂缝的等效宽度 面积 / 骨架长度这里用周长近似 widths [] for i in range(1, n 1): comp labeled i area comp.sum() if area 30: continue perimeter np.logical_xor(comp, ndimage.binary_erosion(comp)).sum() if perimeter 0: widths.append(2 * area / perimeter * pixel_size_mm) return { area_ratio: round(area_ratio, 4), crack_count: n, mean_width_mm: round(float(np.mean(widths)), 3) if widths else 0.0 } print(crack_stats(./pred_masks/slice_0001.png, pixel_size_mm0.05))逻辑说明pixel_size_mm是 CT 扫描的体素实际尺寸必须从扫描参数里拿到否则算出来的宽度没有物理意义。等效宽度用2 * 面积 / 周长近似对细长裂缝够用对分叉裂缝会偏大所以统计时最好按连通域分别看。area_ratio是裂缝面积占比直接对应地质上的裂缝孔隙度贡献。验证统计可信度我一般做两件事。一是抽 20 张测试切片把模型掩码和人工标注掩码分别跑这个脚本对比area_ratio的相对误差控制在 10% 以内算可用。二是把掩码叠加回原图肉眼扫一遍重点看有没有把孔洞误判成裂缝、有没有漏掉大裂缝。这两步做完统计结果才敢往报告里写。最后说个习惯这套流程里模型结构可以换、损失可以调但数据预处理和掩码规范一旦定下来就别轻易动。我见过太多项目模型换了三四版指标上上下下最后发现是标注规范中途改了前后数据不可比。把预处理和标注规范固化成脚本和文档比追 SOTA 模型值钱得多。希望帮到你。本文还有配套的精品资源点击获取
返回列表