简介：基于DenseUnet的岩石盐体图像分割实战项目，为深度学习入门者与地质遥感研究人员提供了一条完整可复现的路径。资源包含Python训练、评估、预测三个核心脚本，代码注释详尽，配合约20MB轻量压缩包，可快速完成从数据准备到模型推理的全流程。包内共2000个文件，以1984张PNG图像为主，涵盖岩石盐体原始图与标注掩膜；另有8个py脚本、5个xml配置文件、2个txt说明和1个README文档，便于理解项目结构与自定义扩展。训练脚本自动生成loss、IoU、学习率衰减曲线及训练日志，帮助用户直观监控模型收敛情况；评估脚本计算测试集的IoU、召回率、精确率等指标；预测脚本输出GT与叠加掩膜图像，方便结果可视化。目前已有56人学习使用，适合需要动手实践分割算法、或希望替换自有数据集进行训练的用户，参照README即可一键运行。
1. DenseUnet盐体分割：从地震剖面到像素级掩膜的实战起点
盐体是油气勘探中的重点构造，但它在二维地震剖面上往往表现为一团高振幅、弱层理的杂乱反射区，边界被断层和低速带磨得模棱两可。基于DenseUnet的盐体图像分割，本质上是把问题转成逐像素二分类：每个点要么属于盐体，要么属于围岩。相比普通U-Net，DenseUnet在编码器里用稠密连接替代单路卷积，让深层网络仍能拿回浅层的细纹理，解码器再用跳跃连接把高层语义与低层空间细节拼在一起。下面会从网络结构、数据预处理、损失函数到推理后处理，把这套方案完整走一遍。适用对象是已经跑通U-Net、想在低信噪比图像上提升分割边界质量，同时习惯看参数和具体踩坑点的工程师。
2. DenseUnet网络结构拆解：Dense Block、Transition Down 与跳跃连接
2.1 Dense Block 的稠密连接：每一层都“见过”前面所有层
DenseNet 的核心思想是：一个Block内，每一层卷积的输出都会和之前所有层的输出在通道维上拼接，而不是像ResNet那样做加法。设Block初始输入通道数为 K0，growth rate 为 g，每层产出 g 个新通道，则第 L 层输入通道数为 K0 + (L-1)*g。这样的好处很直接：反向传播时，梯度可以从网络深度任意位置直接流向早期层，缓解了深层语义网络常见的梯度消失。对盐体分割来说，盐体内部纹理变化不大，真正难的是边界处微弱振幅差异，而稠密连接相当于把不同感受野、不同抽象程度的特征全部暴露在卷积层面前，特征重用率高，同等效果下所需参数更少。
落地到 Keras，一个简洁的Dense Block函数如下：
def dense_block(x, growth_rate=16, n_layers=4): """Dense Block：每一层输出 growth_rate 个特征图并拼到主干上。""" for _ in range(n_layers): y = tf.keras.layers.BatchNormalization()(x) y = tf.keras.layers.ReLU()(y) y = tf.keras.layers.Conv2D(4 * growth_rate, 1, padding='same')(y) # 降维 y = tf.keras.layers.BatchNormalization()(y) y = tf.keras.layers.ReLU()(y) y = tf.keras.layers.Conv2D(growth_rate, 3, padding='same')(y) # 只新增 g 个通道 x = tf.keras.layers.Concatenate()([x, y]) # 通道维拼接 return x这里每个子层先经过BN和ReLU，再用1x1卷积把通道压到 4*g，目的是降低3x3卷积的计算量；随后3x3卷积真正产出 g 个特征图并与输入 x 拼接。参数上，growth_rate 建议在12~24之间取，n_layers 即Block内卷积层数，取4~6层。注意 g 不要贪大，否则第一个Block之后通道数快速膨胀，后层的显存占用会明显增加；在101x101分辨率的地震剖面上，g=16、深度4层的配置通常比 g=32 更稳，后面我会给出完整参数表。
2.2 Transition Down 与 Transition Up：编码器如何降采样、解码器如何恢复分辨率
U形结构的关键在于两件事：编码器逐步把空间尺寸降下来，换取更大的感受野；解码器再逐步把特征图还原到原图尺寸。DenseUnet里，编码器每个Dense Block之后接Transition Down，常见做法是 BN + ReLU + 1x1卷积 + 2x2平均池化或最大池化；解码器则用转置卷积或双线性上采样做Transition Up，并把对应编码器输出拼回来。我一般用平均池化做降采样，因为盐体边界是连续的灰度变化，平均池化能保留背景趋势，而最大池化会更突出局部极值，容易放大噪声。
这个阶段的Keras实现片段：
def transition_down(x, out_channels): x = tf.keras.layers.BatchNormalization()(x) x = tf.keras.layers.ReLU()(x) x = tf.keras.layers.Conv2D(out_channels, 1, padding='same')(x) x = tf.keras.layers.AveragePooling2D(2, strides=2)(x) return x def transition_up(x, out_channels): x = tf.keras.layers.Conv2DTranspose(out_channels, 3, strides=2, padding='same')(x) return xTransition Down的1x1卷积不只是降维，它会在降采样前对Dense Block产出的冗余通道做一次压缩，例如把通道压到上一层的 2/3，控制解码器拼接后的计算量。Transition Up的输出通道一般设成下一层Dense Block的输入，或者编码器对应层的通道数，这样跳跃连接才能对齐通道数直接concat。实际工程里，转置卷积的kernel可以取2或3：取2步长2不会重叠，但缺少跨像素信息；取3步长2会有1个像素重叠，对边界更友好，我推荐3。
2.3 为什么盐体图像分割选DenseUnet而不是普通U-Net
普通U-Net每个层只有一条前向路径，信息要逐层“压缩-恢复”，深层特征图会丢掉很多高频细节；DenseUnet把每个Block内部的特征通道全部保留下来，解码器上采样后拿到的跳跃连接里既有高层语义也有局部纹理。对盐体分割来说，最典型的问题是盐体边界并不总是落在灰度突变处，有时要依靠周围沉积层位的形状推断；DenseNet式连接让网络更容易捕捉到这些跨尺度的上下文。
另一个差别在参数效率：要达到相近精度，DenseUnet往往比同深度U-Net少20%到30%的参数，但训练时显存占用反而更高，因为中间特征图全都被保留下来用于拼接。所以在显存有限时，我通常会把growth_rate调小，而不是减少Block层数，这样精度的损失相对可控。训练时如果发现GPU占用告警，优先检查是否在保存特征图时显存峰值翻倍，而不是急着换更大的显卡。
2.4 DenseUnet关键参数表
下表是盐体分割任务上比较常用的一组起点参数：
| 参数 | 推荐值 | 说明 |
|---|---|---|
| growth_rate | 16 | 每层新增通道数，影响显存和容量 |
| 每Block层数 | 4 | 深度，6层以上要配合dropout |
| Block数量 | 4 | 编码器4个Dense Block，对应4次下采样 |
| 初始通道 | 32 | 第一层卷积通道，数据集小就降到24 |
| Transition Down压缩率 | 0.5 | 把通道压到原来的一半 |
| drop_rate | 0.1~0.2 | 仅在深层Block使用 |
| 输入尺寸 | 128x128 | 训练时统一尺寸，推理可滑窗 |
工程上常见的做法是先把这组参数跑通，再调growth_rate和Block数量。如果测试IOU在0.5以下徘徊，优先看数据预处理而不是网络层数。后面第4章还会给出训练相关参数。
3. 盐体数据集准备与预处理：尺寸归一化、数据增强与生成器写法
3.1 数据目录结构、尺寸与mask格式
盐体分割数据集一般包含两类文件：地震剖面灰度图和同名mask图。地震剖面通常以16位或32位浮点保存，灰度范围不稳定，不能直接除以255作为归一化；常见做法是先把每张剖面按2%和98%分位数做截断，再映射到[0,1]，这样可以顶住噪声尖峰。图像尺寸上，公开的盐体挑战赛数据集的剖面是101x101，但对于DenseUnet，我建议用随机裁剪或双线性插值统一到128x128，方便4次下采样后特征图尺寸为8x8。mask是单通道二值图，0为围岩，1为盐体。
目录结构上保持简单：
salt_dataset/ train/ images/ # 地震剖面，命名 001.png masks/ # 同名mask，命名 001.png val/ images/ masks/这类数据的特殊点在于：同一个数据集里可能有的剖面完全没有盐体，有的盐体占了大半张图。如果直接按文件顺序切训练集和验证集，很容易把某一类样本全分到验证集里。常见做法是按“是否含盐”分层抽样，确保训练和验证都同时包含两类样本，比例尽量接近。
3.2 用 Keras Sequence 写数据生成器
数据量不大，但建议写生成器避免一次性载入内存。用tf.keras.utils.Sequence的子类实现：
import tensorflow as tf import numpy as np import cv2 class SaltGenerator(tf.keras.utils.Sequence): def __init__(self, image_paths, mask_paths, batch_size=8, size=(128, 128), augment=False): self.images = image_paths self.masks = mask_paths self.batch_size = batch_size self.size = size self.augment = augment def __len__(self): return int(np.ceil(len(self.images) / self.batch_size)) def __getitem__(self, idx): batch_img = [] batch_mask = [] for i in range(idx * self.batch_size, min((idx + 1) * self.batch_size, len(self.images))): img = cv2.imread(self.images[i], cv2.IMREAD_UNCHANGED) msk = cv2.imread(self.masks[i], cv2.IMREAD_GRAYSCALE) img, msk = self._preprocess(img, msk) batch_img.append(img) batch_mask.append(msk) x = np.stack(batch_img, axis=0) y = np.stack(batch_mask, axis=0)[..., np.newaxis] # 变成 HxWx1 return x, y def _preprocess(self, img, msk): if img.ndim == 2: img = img[..., np.newaxis] # 统一尺寸 img = cv2.resize(img, self.size, interpolation=cv2.INTER_LINEAR) msk = cv2.resize(msk, self.size, interpolation=cv2.INTER_NEAREST) # 分位数截断归一化 lo, hi = np.percentile(img, [2, 98]) img = np.clip((img - lo) / (hi - lo + 1e-6), 0, 1) if self.augment: img, msk = self._augment(img, msk) return img.astype(np.float32), msk.astype(np.float32)这段代码的要点是：mask使用INTER_NEAREST最近邻插值，避免缩放时盐体边界被插出中间值；图像用双线性插值保留灰度渐变。__getitem__返回的 y 会补一个通道维，因为二分类输出头是1个通道的sigmoid。归一化放在每个样本内做而不是全数据集统一计算，原因在于地震剖面之间测井增益不一致，逐样本截断更符合实际资料特点，但代价是不同剖面的相对灰度强弱被抹掉了；如果地质背景差异不大，问题不大。
3.3 盐体分割专用的数据增强组合
盐体剖面与其他遥感图像不一样，翻转有明确物理意义：地质剖面可以左右平移但不能随便转90度，因为地震数据有垂直深度方向。推荐的增强组合是：随机水平翻转（概率0.5）、随机平移/缩放8%以内、随机0~10度旋转，以及轻度对比度扰动。不建议用弹性形变，因为会破坏层位的平行关系。代码里可以在_augment方法中使用cv2.warpAffine配合随机变换矩阵。
def _augment(self, img, msk): rows, cols = img.shape[:2] # 水平翻转 if np.random.rand() < 0.5: img = cv2.flip(img, 1) msk = cv2.flip(msk, 1) # 小角度旋转和平移 angle = np.random.uniform(-10, 10) dx = np.random.uniform(-0.08, 0.08) * cols dy = np.random.uniform(-0.08, 0.08) * rows mat = cv2.getRotationMatrix2D((cols/2, rows/2), angle, 1.0) mat[0, 2] += dx mat[1, 2] += dy img = cv2.warpAffine(img, mat, (cols, rows), flags=cv2.INTER_LINEAR) msk = cv2.warpAffine(msk, mat, (cols, rows), flags=cv2.INTER_NEAREST) return img, mskangle取±10度而不是90度，是因为地震剖面的垂向是深度，旋转90度后的样本物理上不成立。若你的数据源包含三维地震体的inline切片，才可以把90度旋转当作数据扩充。增强概率和幅度不宜过大，否则模型会花容量去学习无关形变；判别增强是否合适的标准是增强后的mask和图像边界是否仍然大致吻合。如果发现盐体边缘被旋转“切”出锯齿，可以加一点mask腐蚀后处理。
| 增强方式 | 参数 | 说明 |
|---|---|---|
| 水平翻转 | p=0.5 | 沿测线方向对称，保留层理 |
| 旋转 | ±10度 | 超出范围会破坏深度方向 |
| 平移缩放 | ±8% | 模拟剖面位置轻微偏移 |
| 对比度扰动 | 增益0.8~1.2 | 补偿不同测井增益差异 |
4. 训练DenseUnet的关键配置：损失函数、学习率与类别不平衡
4.1 BCE + Dice组合损失函数怎么写
盐体在整幅图里通常只占10%到30%，很多样本甚至完全没有盐体，普通二分类交叉熵（BCE）会让网络快速学会“全部预测为0”，因为空图的loss最低。Dice Loss对前景占比不敏感，更适合小目标，但单独用Dice会导致训练早期梯度不平滑。常见做法是把两者按比例相加：
def dice_loss(y_true, y_pred, smooth=1e-6): y_true = tf.cast(y_true, tf.float32) y_pred = tf.sigmoid(y_pred) intersection = tf.reduce_sum(y_true * y_pred) return 1 - (2.0 * intersection + smooth) / (tf.reduce_sum(y_true) + tf.reduce_sum(y_pred) + smooth) def combined_loss(y_true, y_pred): bce = tf.keras.losses.BinaryCrossentropy()(y_true, y_pred) return bce + 0.5 * dice_loss(y_true, y_pred)这里的 y_pred 是未过sigmoid的logits，dice_loss 内部计算sigmoid，避免与BCE的数值稳定性冲突。权重0.5表示Dice部分只起辅助作用；如果空盐体样本占比超过50%，可以把权重提到1.0或改成 BCE + 2*Dice，让网络更重视形态匹配。观察训练日志时注意，loss数值的绝对值没有直接意义，要看训练集与验证集的Dice或IOU趋势。
4.2 模型编译与训练调参表
我一般用Adam，初始学习率1e-3，配合ReduceLROnPlateau在验证损失三个epoch不降时把学习率减半。batch size在512MB显存上取8，128x128输入下相对稳妥。
model.compile(optimizer=tf.keras.optimizers.Adam(1e-3), loss=combined_loss, metrics=[tf.keras.metrics.MeanIoU(num_classes=2, name='miou'), tf.keras.metrics.BinaryAccuracy(name='acc')]) callbacks = [ tf.keras.callbacks.ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=3, min_lr=1e-6), tf.keras.callbacks.EarlyStopping(monitor='val_miou', patience=10, restore_best_weights=True), tf.keras.callbacks.ModelCheckpoint('denseunet_salt_{epoch:02d}_{val_miou:.3f}.h5', monitor='val_miou', save_best_only=True) ] model.fit(train_gen, validation_data=val_gen, epochs=60, callbacks=callbacks)MeanIoU在这里的 num_classes=2，但输出通道只有1个，所以IOU计算会自动把预测概率以0.5为阈值转成类别再统计。实测下来，训练集miou超过0.7、验证集还在0.5附近徘徊，多半是过拟合，解决办法不是加更多数据，而是把 drop_rate 从0.1提高到0.3，或在Transition Down里加SpatialDropout2D。如果验证集一直小于0.4，先检查mask是否错位，或者归一化是否把盐体区域压没了。
| 训练参数 | 推荐值 | 适配场景 |
|---|---|---|
| optimizer | Adam | 默认即可，不需要Nadam |
| 初始学习率 | 1e-3 | 配合衰减 |
| batch_size | 8 | 128x128输入，显存紧张就降到4 |
| 最大epochs | 60 | 配合EarlyStopping使用 |
| 学习率衰减 | 验证损失3轮不降减半 | 盐体分割用loss比miou更稳 |
4.3 空盐体样本如何处理
空盐体样本不能直接丢弃，因为它教会模型“什么不是盐”。但大批量空样本会让Dice Loss变成常数，干扰优化。常见的做法有两个：一是按“含盐样本与空样本比约为2:1”重新采样训练集；二是在损失里加一个类别权重，比如给正样本更高的权重。我用的是重新采样，在构造数据集时先统计每个mask的面积，把含盐样本复制两份并入列表。另一种更简单的方式是在生成器里按空样本概率跳过，但要注意不能跳过太多，否则每个epoch的有效样本数不等于列表长度。
| 空盐体占比 | 处理方案 | 适用阶段 |
|---|---|---|
| 小于20% | 不必特殊处理 | 小规模试跑 |
| 20% ~ 50% | 损失加权或重采样 | 正常训练 |
| 大于50% | 重采样确保每个batch至少有1张含盐图 | 工程落地 |
4.4 训练曲线怎么看、什么时候停
重点看 val_miou 和 train_miou 的gap。gap小于0.05说明模型刚从欠拟合走出来，可以继续训练；gap在0.1到0.2之间说明开始过拟合，EarlyStopping生效后会自动恢复最佳权重；gap超过0.25基本可以断定数据分布出了问题。
盐体分割里最常见的诡异曲线是：前几个epoch miou是0，因为Dice Loss在初始阶段梯度很小，logits接近0.5附近，预测出的mask全是半透明灰色，IOU按0.5阈值一算就是0。这时不要停，等10个epoch左右再看。如果BCE Loss在降但Dice不降，把学习率调低一半再续训。
5. 推理与后处理：把预测掩膜变成干净盐体区域的三个技巧
5.1 滑窗推理避免显存瓶颈
训练用的128x128输入到了推理阶段不一定够。完整地震剖面通常有几百到几千列，直接整图预测很容易超出显存或丢失盐体的整体连续性。常见做法是128x128滑窗、重叠32像素，把每个窗口的预测概率按位置平均后再拼回原图。重叠区域取平均可以有效抑制窗口边缘的条纹伪影，这在地震剖面上特别明显。
stride = 96 probs = np.zeros((H, W), dtype=np.float32) weight = np.zeros_like(probs) for y in range(0, H - size + 1, stride): for x in range(0, W - size + 1, stride): patch = cv2.resize(img[y:y+size, x:x+size], (128, 128)) p = model.predict(patch[..., np.newaxis][np.newaxis, ...], verbose=0) probs[y:y+size, x:x+size] += p[0, ..., 0] weight[y:y+size, x:x+size] += 1.0 probs /= np.maximum(weight, 1e-6)stride=96表示步长小于窗口尺寸，即窗口间有32像素重叠，相当于每张剖面最终输出概率是多个窗口预测的加权平均。注意边缘不满128的窗口需要pad，推理后再裁掉，否则剖面上下边界会出现异常高亮。显存充足时可以把步长减小到64，重叠更多，结果更平滑，但推理时间会翻倍。
5.2 后处理三件套：阈值、连通域过滤、形态学闭合
模型输出的连续概率要先转二值。阈值不是固定的0.5，我常用0.55或0.6，因为盐体识别偏保守更实用，宁可漏掉边缘也不能把非盐体圈进来；建议在验证集上扫一遍0.4到0.8选出最优阈值。第二步是基于连通域过滤：统计每个盐体连通域的面积，删除小于面积阈值的孤立点，比如小于100像素的“盐体”多半是噪声。第三步是形态学闭合，用5x5或7x7的结构元素把盐体内部断裂的小孔隙填上，但不要做腐蚀，因为盐体边界已经很细，腐蚀会把真实边缘削掉。
ret, binary = cv2.threshold(prob, 0.55, 1, cv2.THRESH_BINARY) n, labels, stats, _ = cv2.connectedComponentsWithStats((binary * 255).astype(np.uint8), 8) for i in range(1, n): if stats[i, cv2.CC_STAT_AREA] < 100: # 小连通域当噪声删除 binary[labels == i] = 0 kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5, 5)) binary = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel)连通域按8邻域统计，面积阈值要按图像分辨率调整：如果把剖面插值到512宽，面积阈值可以提到500。处理完后再对比原始地震剖面，重点看盐体边界是否把断层切出来的非盐区圈进来了。如果出现又大又散的异常连通域，多半是阈值偏低，先用0.6再跑一遍。
5.3 验证：用IOU复算掩膜质量
后处理改变了预测结果，所以不能再用训练时的miou来说明最终效果。常见做法是单独写一个函数，对验证集逐样本计算标准IOU，同时统计“空盐体预测为空”的准确率。容易被忽略的是IOU会把大盐体和小盐体等同看待：如果验证集里有一张盐体面积达到全图30%的大目标，它的IOU会主导平均值，掩盖小盐体效果差的问题。所以我会按mask面积分层统计，把验证集分成1000像素以下、1000到5000像素、5000像素以上三档，分别报IOU。DenseUnet通常在5000像素以上的大盐体上表现最好，小盐体更多靠后处理保留。把验证脚本保留下来，跑一遍完整验证集，记录每一档的IOU分布，再决定要不要回去调growth_rate或后处理阈值。
本文还有配套的精品资源，点击获取
