简介:一种面向地震数据断层解释的深度学习方法研究文献,针对传统人机交互解释效率低、常规算法依赖参数设置等问题,提出结合ResNet残差网络与U-Net编码-解码架构的SeisFault-Net模型,实现端到端断层自动识别。资源为《基于深度卷积神经网络的地震数据断层识别方法》单篇PDF论文,共1个PDF文件,压缩包大小6.49MB,适合地球物理勘探、深度学习应用及地震资料解释方向的研究人员与技术工程师学习参考。文中详细阐述了残差学习块原理、U-Net多尺度特征提取结构及网络训练策略,并通过实际数据实验对比相干算法,验证了该方法在断层垂向连续性和轮廓清晰度方面的优势;同时包含引言中对互相关、相似性、特征结构相干等常规方法的总结,有助于读者系统理解断层识别技术演进。该文档已有272人学习,可作为智能地球物理领域论文研读、方法复现或课堂讨论的补充材料。
1. 断层识别为什么从人工拾取转向深度卷积神经网络
在三维地震解释里,断层拾取一直是个“费力不讨好”的活。解释员盯着相干体切片,一条一条勾断层多边形,一个中型工区几百条断层,几个人干几周是常态;更麻烦的是,不同人拾取的断层位置和延伸范围差异很大,同一个断点,A觉得是、B觉得不是,成果质量全靠个人经验兜底。深度卷积神经网络(CNN)进入这个领域后,把断层识别变成了一个逐体素的图像分割问题:输入地震振幅数据或相干属性,输出一张断层概率体,位置、倾向、延伸范围一次性给全,几分钟跑完一个工区,而且结果可复现。这篇笔记就围绕“基于深度卷积神经网络的地震数据断层识别方法”这个方向,把网络怎么选、训练数据怎么造、参数怎么调、生产落地有哪些坑,一条线讲透。
这套方法适合两类人:一类是被动辄上千平方公里的断层解释压得喘不过气的油藏地球物理工程师,想找个能自动出初稿的工具;另一类是正在做深度学习地震应用研究的在校生或算法工程师,需要一份能快速上手的落地路径,而不是从零搭网络、踩一遍已经被踩烂的坑。
2. 网络怎么选:断层识别不是分类题,是分割题
2.1 为什么普通CNN分类网络做不了断层识别
最早很多人把断层识别当分类任务做:取一个固定大小的窗口,比如32×32×32的体素块,CNN判断窗口中心是否落在断层上。这个思路在小块验证集上效果还不错,一旦放到整口三维工区就废了——断层在三维空间里是一个面,面在体数据里只占极小比例,窗口一滑,绝大多数是背景体素,训练集严重偏斜;而且窗口中心贴着断面但有位移错动时,模型很容易被周围的地层纹理带偏,把层位错断误判成断层。
断层识别本质上是逐体素的像素级分类,也就是语义分割。每个体素都要输出一个“是不是断层”的概率,而不是整个窗口给一个类别。这决定了网络骨架必须带有编码器-解码器结构:编码器逐层下采样提取多尺度特征,解码器把低分辨率特征图上采样回原始分辨率,中间用跳跃连接把浅层细节补回来。U-Net就是这个结构的典型代表,也是我做断层识别默认的首选骨架,原因很简单:它参数量不大、显存占用可控,而且对薄目标的分割效果好——断层在剖面上往往只有几个像素宽,普通分类网络下采样两三次就把这条线抹没了,U-Net的跳跃连接能把这些细节捞回来。
2.2 从U-Net到Attention U-Net,改动点到底在哪
基线U-Net够用,但断层识别有两个特殊问题它处理得不够好。第一是断层体素占比极低,一张256×256的切片上断层像素经常不到5%,U-Net的普通交叉熵损失会倾向于把所有像素预测成背景,因为全猜背景也有95%的正确率;第二是断层附近有大量“类断层”干扰,比如不整合面、河道砂体的边界、相干噪声,这些位置纹理变化剧烈,模型容易给高置信度的假正例。
我的做法是两处改动。第一处是损失函数换成Focal Loss,它给难分类样本更高的权重,公式核心是给交叉熵乘一个调制因子(1-p_t)^γ,γ取2.0,模型就会被迫关注那些被分错的断层体素,而不是躺在大片背景上拿高准确率。第二处是加一个空间注意力模块,具体做法是在解码器的最后一层前插入一个1×1卷积生成空间权重图,让网络学会在断层附近提高注意力、在平坦层位区降低响应。Attention U-Net之所以比普通U-Net在这个任务上稳定,不是因为结构更复杂,而是它把“该看哪里”这件事显式建模了。Transformer结构这两年也有人试,但地震数据标注量通常只有几千个patch,Transformer在这种小数据下的优势发挥不出来,还容易过拟合,我不建议作为首选。
2.3 一个能直接跑的U-Net断层识别模型结构
下面这个模型是3D U-Net的轻量版,输入64×64×64的振幅patch,输出同尺寸的断层概率体。depth设为4,即下采样4次,特征通道从16起步,每次翻倍增加到256。我用PyTorch写,激活函数用LeakyReLU而不是ReLU,因为ReLU在训练初期容易让部分神经元死掉,LeakyReLU斜率0.1就能避免这个问题。
import torch import torch.nn as nn class ConvBlock(nn.Module): def __init__(self, in_ch, out_ch): super().__init__() self.conv1 = nn.Conv3d(in_ch, out_ch, kernel_size=3, padding=1) self.conv2 = nn.Conv3d(out_ch, out_ch, kernel_size=3, padding=1) self.act = nn.LeakyReLU(0.1, inplace=True) def forward(self, x): x = self.act(self.conv1(x)) x = self.act(self.conv2(x)) return x class UNet3D(nn.Module): def __init__(self, in_channels=1, base_filters=16, depth=4): super().__init__() self.encoders = nn.ModuleList() self.decoders = nn.ModuleList() self.pools = nn.ModuleList() ch = in_channels for i in range(depth): out_ch = base_filters * (2 ** i) self.encoders.append(ConvBlock(ch, out_ch)) self.pools.append(nn.MaxPool3d(2)) ch = out_ch self.bottleneck = ConvBlock(ch, ch * 2) ch = ch * 2 for i in range(depth - 1, -1, -1): out_ch = base_filters * (2 ** i) self.decoders.append( nn.Sequential( nn.ConvTranspose3d(ch, out_ch, kernel_size=2, stride=2), ConvBlock(ch, out_ch) ) ) ch = out_ch self.out_conv = nn.Sequential( nn.Conv3d(ch, 1, kernel_size=1), nn.Sigmoid() ) def forward(self, x): skips = [] for enc, pool in zip(self.encoders, self.pools): x = enc(x) skips.append(x) x = pool(x) x = self.bottleneck(x) for i, dec in enumerate(self.decoders): x = dec[0](x) skip = skips[len(skips) - 1 - i] x = torch.cat([x, skip], dim=1) x = dec[1](x) return self.out_conv(x)ConvTranspose3d是转置卷积,作用是上采样。这里用kernel_size=2、stride=2,正好把特征图尺寸翻倍,然后和对应下采样层的输出在通道维拼接。注意拼接前skip和x的空间尺寸必须一致,如果输入patch边长是偶数,MaxPool3d和ConvTranspose3d的配对不会出问题;一旦输入变成奇数尺寸,上采样后尺寸对不上,会直接报错。这就是我固定用64、128这类偶数的原因,训练和推理都省心。
2.4 Focal Loss的实现与参数选择
网络结构定了,训练时真正决定上限的是损失函数和采样策略。Focal Loss的PyTorch实现很短,但参数不能照搬目标检测那套默认值。
import torch.nn.functional as F class FocalLoss(nn.Module): def __init__(self, alpha=0.25, gamma=2.0): super().__init__() self.alpha = alpha self.gamma = gamma def forward(self, pred, target): ce_loss = F.binary_cross_entropy(pred, target, reduction='none') p_t = pred * target + (1 - pred) * (1 - target) modulating_factor = (1 - p_t) ** self.gamma alpha_factor = target * self.alpha + (1 - target) * (1 - self.alpha) loss = alpha_factor * modulating_factor * ce_loss return loss.mean()alpha是正负样本的权重平衡因子,0.25表示负样本权重是0.75、正样本是0.25,只对正样本降权,因为断层体素本来就少,模型容易忽视它们;gamma是难例聚焦参数,gamma越大,对已经分对的样本降权越狠。断层识别任务上,gamma用2.0没问题,alpha建议从0.5开始调,往0.25方向走——如果训练曲线里背景体素的loss降得飞快而断层loss纹丝不动,就得继续降低alpha或者提高gamma。一个容易被忽略的细节是pred在进FocalLoss之前必须经过Sigmoid,这个在模型输出层已经处理了,不要在损失函数里再套一次Sigmoid,否则梯度会翻车。
3. 训练数据从哪来:合成地震数据是断层识别的发动机
3.1 人工标注为什么撑不起深度学习训练集
断层识别的训练数据是个老大难。真实工区的地震数据虽然有断层,但断层标签需要人工解释,解释员在剖面上一层层勾,勾出来的断层位置误差有几个体素宽,而且不同人勾的差异很大,这样的标签喂给网络,学出来的断层边界是模糊的,预测概率体也是模糊一片。更关键的是,真实工区的断层样本有限——一个三维体里可能只有几十条断层,要训练一个几百万参数的3D U-Net,数据量完全不够。
所以这个方向的主流做法是用合成地震数据做训练。通过构造地质模型,设定断层的位置、倾角、走向、断距,然后正演生成地震响应。合成数据的好处有三个:标签是精确到体素的,断层面的位置是数学定义的,没有人工拾取误差;断层的形态、规模可以随意组合,想生成多少条就生成多少条;可以控制噪声、子波频率、地层倾角来制造多样性,让模型在合成数据上见过足够多的变化。这就是这个标题背后最核心的工程逻辑:与其费劲标真实数据,不如批量生成带精确标签的合成数据。
3.2 生成带断层的三维地质模型
断层模型可以简化成这样的思路:先在一个三维网格里生成若干层水平地层,每层给一个随机的厚度和速度扰动,然后把断层面定义为一个平面或曲面,把断面一侧的全部地层沿断层法向平移一个距离——这个位移量就是断距。地层被错开的位置,就形成了天然的断层标签。下面这段Python代码生成一个64×64×64的模型,包含两个不同倾角的断层。
import numpy as np def generate_faulted_model(size=64, num_faults=2): # 初始化层位编号和速度模型 depth = np.arange(size) layers = np.zeros((size, size, size), dtype=np.int32) velocity = np.zeros((size, size, size), dtype=np.float32) # 随机层界面:每层厚度在6~15个体素之间波动 boundaries = [0] while boundaries[-1] < size - 10: boundaries.append(boundaries[-1] + np.random.randint(6, 16)) boundaries.append(size) for i in range(len(boundaries) - 1): layers[boundaries[i]:boundaries[i+1], :, :] = i velocity[boundaries[i]:boundaries[i+1], :, :] = 2500 + 200 * i + np.random.normal(0, 30, (size, size, size))[boundaries[i]:boundaries[i+1], :, :] fault_mask = np.zeros((size, size, size), dtype=np.uint8) for f in range(num_faults): dip = np.random.uniform(30, 60) # 倾角30~60度 strike_offset = np.random.uniform(-10, 10) # 沿走向的位移 x0, z0 = np.random.randint(15, 45, 2) fault_plane = np.zeros((size, size, size), dtype=bool) for z in range(size): x_center = x0 + int((z - z0) / np.tan(np.radians(dip))) for x in range(max(0, x_center - 1), min(size, x_center + 2)): fault_plane[x, :, z] = True fault_mask[fault_plane] = 1 # 断面一侧地层整体平移,制造错断 shift = np.random.randint(4, 10) side = np.random.choice([1, -1]) for x in range(size): if side * (x - x0) > 0: layers[x, :, :] = np.roll(layers[x, :, :], shift, axis=0) return layers, fault_mask这段代码的核心是fault_plane标记断面附近的体素作为标签,然后通过np.roll把断面一侧的地层沿深度方向滚动,模拟断距。np.roll是循环位移,超过边界的数据会从另一头绕回来,这在模型边缘会产生伪断层,所以实际做正演前要把边缘几十个体素裁剪掉或者用padding补上,不能直接拿这个模型做训练样本。倾角范围30到60度、断距4到10个体素这些参数,可以覆盖大多数正断层的形态。
3.3 褶积正演,把地质模型变成地震剖面
有了速度模型和断层标签,下一步要把它变成“看起来像地震数据”的东西。最常用的简化正演是褶积模型:假设地下反射系数由速度差异决定,然后用一个雷克子波去卷积,得到合成地震记录。公式不复杂,但代码里有两个坑要提前避掉:子波主频不能太高,否则合成记录太细碎,断层识别模型学到的特征和真实地震不匹配;反射系数的计算要对速度做对数差分,而不是直接对速度做差分。
from scipy import signal def generate_seismic(layers, velocity, wavelet_freq=35): size = layers.shape[0] reflectivity = np.zeros_like(velocity, dtype=np.float32) # 反射系数 = 相邻层速度差 / 平均速度的2倍(纵波垂直入射近似) v_shifted = np.roll(velocity, 1, axis=0) reflectivity[1:] = (velocity[1:] - v_shifted[1:]) / (velocity[1:] + v_shifted[1:]) * 2 # 雷克子波,主频35Hz,采样间隔1ms,对称长度64 t = np.arange(-32, 32) * 0.001 w = (1 - 2 * (np.pi * wavelet_freq * t) ** 2) * np.exp(-(np.pi * wavelet_freq * t) ** 2) seismic = np.zeros_like(reflectivity) for x in range(size): for z in range(size): seismic[x, :, z] = signal.fftconvolve(reflectivity[x, :, z], w, mode='same') # 加入高斯噪声和轻微空间平滑,模拟真实采集噪声 seismic += np.random.normal(0, 0.02 * np.std(seismic), seismic.shape) return seismic反射系数公式用的是垂直入射近似,没有考虑透射损失和多次波,这是褶积模型的固有局限,但对断层识别来说够用——断层识别的核心特征是地层同相轴的中断和错动,褶积模型已经能把这两点表达清楚。fftconvolve用快速傅里叶变换做卷积,比直接循环卷积快一个数量级;mode='same'保证输出长度和输入一致,不会因为卷积而缩短。噪声强度0.02倍标准差是经验值,加太小模型会过拟合干净合成数据,加太大真实数据的噪声分布又对不上,建议在0.01到0.05之间扫一两个档对比效果。
3.4 训练样本的多样性,决定模型的泛化能力
我见过不少团队把断层识别模型做出来,合成数据上效果惊艳,一到真实工区就垮,多半是训练样本多样性不够。多样性至少包含四个维度:断层倾角要从15度覆盖到75度,不能只生成中角度正断层;断距要从2个体素覆盖到15个体素,小断距是识别难点,大断距容易让模型只学大落差;地层倾角要有平地层也要有倾斜地层,倾斜层位本身就有同相轴斜交,模型得学会区分层位斜交和断层错断;噪声水平要有低噪和高噪,实际地震资料深部信噪比极低,只在干净合成数据上训练等于没训练。常见做法是每个epoch重新生成一批随机模型,而不是用一个固定训练集反复迭代,这样模型每个batch看到的断层形态都是新的,相当于免费的数据增强。这个技巧很便宜,但效果立竿见影——用固定数据集训练的模型在真实数据上的断层连续性通常要差20%以上。
4. 训练与推理:参数设置和生产落地的完整流程
4.1 训练时的关键参数:patch尺寸、batch和优化器
断层识别训练不是把整个三维体直接塞进网络,显存也装不下,必须切patch。patch尺寸我固定用64×64×64,这个尺寸在感受野覆盖和显存消耗之间比较均衡:小于64,断层两侧的地层关系看不全,模型容易把局部扰动当断层;大于128,显存直接爆,而且一个patch里包含多个断层,样本标签统计更难平衡。batch size在单张24G显存卡上设为4比较稳,如果显存不够,优先减batch而不是减patch尺寸。优化器用AdamW,初始学习率1e-3,weight decay设1e-4——AdamW比Adam在训练后期更稳,断层识别这种小目标任务后期loss波动明显,weight decay能起到约束作用。学习率用CosineAnnealing,从1e-3降到1e-5,比StepLR那种阶梯式下降更容易收敛到平滑的局部最优。
4.2 训练时的数据增强:不只是翻转和旋转
合成数据的增强策略和自然图像不一样,需要针对地震数据的物理特性设计。翻转和旋转90度是必须的,地震数据没有上下左右的方向先验,翻转能让模型学到更本质的断层几何特征。除了这些常规操作,还有两个针对性的增强很关键。第一个是振幅增益扰动——真实地震体不同区块的振幅尺度差异很大,有的区域强振幅,有的弱振幅,训练时对patch乘一个0.8到1.2的随机增益,模型就不会把绝对振幅当作断层判据。第二个是波形相位扰动——对patch做随机小角度相位旋转,模拟不同子波相位的影响,这比单纯加噪声更贴近真实资料的多样性。裁剪位置也很有讲究:断层在三维体里是稀疏的,如果完全随机裁剪,很多patch里压根没有断层,浪费计算;常见做法是先让每个patch至少包含一条断层面的一部分,另外再配30%的纯背景patch用于训练负样本。
4.3 推理阶段的滑窗策略:重叠窗口怎么拼不出拼接缝
训练完成后,推理要在整个三维体上逐patch预测。如果每个patch独立预测再拼起来,patch边界处断层概率会出现明显的拼接缝——原因很简单,patch边缘的体素感受野不完整,模型在边缘的预测置信度天然低于中心。解决方法是滑窗重叠加权重融合。常见做法是patch尺寸64,stride设32,每个体素会被多个patch覆盖,预测概率按覆盖次数或权重累加,取平均值。边缘体素权重低、中心体素权重高,加一个三角权重函数在重叠区做平滑过渡,拼接缝基本可以抹掉。
def sliding_window_inference(model, volume, patch_size=64, stride=32): model.eval() cz, cx, cy = volume.shape accum = np.zeros(volume.shape, dtype=np.float32) weight_map = np.zeros(volume.shape, dtype=np.float32) # 三角权重:从patch边缘到中心线性增权 w = np.hanning(patch_size) w3d = np.outer(w, w) w3d = np.stack([w3d * w[i] for i in range(patch_size)], axis=0) w3d = (w3d - w3d.min()) / (w3d.max() - w3d.min()) + 0.01 for z in range(0, cz - patch_size + 1, stride): for x in range(0, cx - patch_size + 1, stride): for y in range(0, cy - patch_size + 1, stride): patch = volume[z:z+patch_size, x:x+patch_size, y:y+patch_size] tensor = torch.from_numpy(patch).unsqueeze(0).unsqueeze(0).float() with torch.no_grad(): pred = model(tensor).squeeze().numpy() accum[z:z+patch_size, x:x+patch_size, y:y+patch_size] += pred * w3d weight_map[z:z+patch_size, x:x+patch_size, y:y+patch_size] += w3d # 最外层stride覆盖不到的位置用最近邻补 idx_z = min((cz - patch_size) // stride * stride, cz - patch_size) idx_x = min((cx - patch_size) // stride * stride, cx - patch_size) idx_y = min((cy - patch_size) // stride * stride, cy - patch_size) for z in range(idx_z, cz - patch_size + 1): for x in range(idx_x, cx - patch_size + 1): for y in range(idx_y, cy - patch_size + 1): if weight_map[z, x, y] == 0: patch = volume[z:z+patch_size, x:x+patch_size, y:y+patch_size] tensor = torch.from_numpy(patch).unsqueeze(0).unsqueeze(0).float() with torch.no_grad(): pred = model(tensor).squeeze().numpy() accum[z:z+patch_size, x:x+patch_size, y:y+patch_size] += pred weight_map[z:z+patch_size, x:x+patch_size, y:y+patch_size] += 1 result = accum / np.maximum(weight_map, 1e-6) return result这段代码的三重循环看着笨重,但对于一个1280×1280×512的真实工区,Patch-based推理是唯一可行的方案。权重用np.hanning生成的是二维窗,再沿深度轴堆叠成三维权重,这样做的好处是patch的六个面都是低权重,避免了拼接缝集中在某一条边上。最后一个补丁处理的是三维体尺寸不能被stride整除的情况——整块数据如果没有补丁覆盖,概率图会出现规则的黑格子,这个坑我踩过,后来加了最后这个兜底循环才彻底解决。
4.4 输出后处理:从概率体到断层骨架
模型输出的是0到1的断层概率体,但解释员需要的是一条条断层,不是一团模糊概率云。后处理一般分三步走。第一步是阈值化:概率大于0.5的体素标记为断层候选点,但这个阈值不是固定的,我更习惯先看概率体的直方图,取双峰之间的谷底作为阈值——断层概率体通常是双峰分布,背景峰在0.1附近,断层峰在0.7附近,谷底一般落在0.35到0.5之间。第二步是形态学处理:用三维连通域分析去掉孤立的噪声体素(小于100个体素的连通域直接删除),再用闭运算把断裂的断层片段连接起来——三维闭运算的腐蚀膨胀核用3×3×3的立方体就够了,核太大会把相邻断层粘连在一起。第三步是骨架化:用形态学细化算法把断层带压缩成单层体素的骨架线,这一步得到的骨架就是断层面的近似表示,可以直接导出成解释软件的断层多边形或者地震属性体。
5. 避坑手册:断层识别最常见的4个翻车现场
5.1 断层和层位混淆:倾角约束为什么不能只在网络里做
现象:模型预测结果里,断层概率体上出现了大量与地层层理平行的连续高亮带,尤其在倾斜地层区,层位错断处和真实断层混在一起,解释员没法用。
原因:单纯靠卷积网络很难区分“地层错断”和“不整合面/河道切割”。这两者在局部纹理特征上高度相似:都是同相轴突然中断,振幅出现不连续。网络在局部patch上看到的信息不足以判断这是断层面还是剥蚀面。
解决:不能只在网络结构里想办法,要在两个层面同时加约束。训练层面,合成数据里必须加入地层褶皱、不整合面、河道充填等非断层的不连续体,明确标注为背景,让模型见过足够多的混淆样本;后处理层面,利用断层的几何先验——真实断层面是近似平面的曲面,倾角变化是渐变的,而层序边界是近乎水平的——把断层概率体做一次倾角滤波,先计算每个体素的局部倾角,把倾角小于10度的响应全部清零,再做骨架化。这个倾角滤波几乎不影响真实断层,但能干掉80%的层位假正例。
5.2 小断距断层被模型直接忽略
现象:断距小于3个体素的断层在预测结果里只剩断断续续的几条短线,甚至完全消失,剖面上去看,断层位置只有微弱的概率响应,达不到阈值。
原因:小断距断层的振幅变化微弱,卷积下采样两次之后,这点微弱差异就被池化操作抹掉了。另外Focal Loss里gamma=2.0对小断距样本也不友好——小断距断层在特征层面和背景很接近,模型对它们给出的概率往往在0.3到0.5之间,调制因子让这些本就不确定的样本权重降得更低,模型就更不愿意学了。
解决:第一,把gamma从2.0降到1.0,让小断距样本的梯度不要被过度压制;第二,在损失函数里显式加入断距大小的权重因子——合成数据生成时记录每个断层体素的断距,训练时把断距小于5个体素的断层体素的loss权重乘1.5;第三,推理时不要只用一个固定阈值,可以输出两个级的概率图,低阈值0.25的结果用于发现弱信号,高阈值0.6的结果用于精确定位,解释员在看弱断层时叠加低阈值图参考。小断距断层是断层识别落地最大的痛点,这块不做专项优化,模型做个demo可以,生产完全不行。
5.3 合成数据到真实数据,预测结果整体漂移
现象:合成数据上验证Dice能达到0.8,同一套权重迁移到真实工区,断层概率整体偏低,原来0.5的阈值切出来的断层稀稀拉拉,阈值调到0.3才勉强能看,但噪声也跟着上来了。
原因:域差异。合成数据用的是理想化褶积模型,子波是零相位雷克子波,真实地震数据经历了吸收衰减、多次波、各向异性传播,波形和振幅关系比合成数据复杂得多。纯粹用合成数据训练的权重,在真实数据上的特征分布发生了偏移。
解决:域适应要分两步走。第一步是训练时的数据增强尽量向真实数据靠拢——相位扰动范围加大,振幅增益范围放宽到0.6到1.5,噪声从高斯噪声换成带通滤波后的有色噪声,这些都能缩小域差距。第二步是微调,用少量真实工区数据对预训练权重做迁移学习——不需要真实标签,用模型在真实数据上的预测置信度做伪标签,挑置信度高的一部分断层体素和确认是背景的体素,继续训练几十个epoch。这个半监督微调套路能把真实工区的预测质量拉回可用水平,比重新标定阈值管用得多。注意微调学习率要降到1e-5以下,不然会把预训练学到的断层几何知识冲掉。
5.4 断层连续性问题:概率体上断层一段亮一段暗
现象:一条延伸几十公里的断层,预测出来的概率体上有若干段概率值明显偏低,导致连通域分析时断层被切成好几段,骨架化无法生成完整断层线。
原因:断层沿走向的位移是变化的,断层中部位移最大,两端逐渐消失,这是断层的自然生长机制。网络对位移大的地方预测置信度高,位移小的尾端几乎和背景没区别,这就造成了概率体沿断层走向的不均匀。
解决:后处理环节引入方向性增强——计算断层概率体的梯度结构张量,在断层延伸方向做平滑,断层法向方向保持尖锐。常见做法是对概率体做各向异性扩散滤波,沿局部主方向扩散,垂直主方向不扩散,这样能把断层从暗处往亮处补连起来。但要注意扩散迭代次数不能太多,我一般控制在10到20次,太多会把两条平行断层扩散成一条。如果做完方向性增强之后断层还是断的,就得回到训练数据里检查是不是断层延伸长度分布不够多样——合成数据里断层长度如果都集中在几十个体素,模型对长断层的连续性就没有概念。
6. 生产验证:从概率体到解释成果的有用性检查
模型训练完、后处理流程也通了,最后一个问题是怎么证明这套结果真的能用。我一般做三件事。第一件事是定量对比:选一个工区里解释员已经解释过的断层多边形,和模型输出的断层骨架做对比,统计模型检出的断层数量、总延伸长度,以及和人工解释的重合率。重合率不是简单算Dice,而是算“人工断层有多少比例落在模型断层的10个体素邻域内”这个召回率,再加一个“模型断层有多少比例落在人工解释邻域内”的精确率。如果召回率能到80%以上、精确率在60%以上,这个模型就有实际生产的价值了。第二件事是交叉验证:同一个模型在相邻工区跑一遍,看看断层空间分布是否和区域构造背景吻合——如果一个工区完全没有断层的地方冒出一堆高概率响应,大概率是模型吃进了振幅异常,要回炉检查训练数据。第三件事最关键,是把模型的断层骨架和人工解释的断层属性交叠显示在三维空间里,让解释员抽查——断层倾向是否合理、断距大小是否和剖面匹配、是否切到了不该切的层位。模型输出的断层骨架是初稿不是终稿,解释员在此基础上修改的效率,比从零开始拾取要高得多。
这套方法做下来,最深的体会是:深度学习断层识别七分在数据、两分在后处理、一分在网络结构。网络结构卷积层加几层减几层,影响远没有合成数据的多样性和后处理参数大。我习惯在每个工区开工前,先用模型跑一个快速测试,看看概率体的双峰分布是否清晰,如果两个峰叠在一起分不开,先别急着调网络,回去看看训练数据里的断层倾角范围、噪声水平是不是和这个工区差太远。合成数据不是选出最优的那一套就完事了,它是一个迭代器,需要根据每个工区的资料特点持续补充。希望这篇文章能帮你绕过我踩过的那些坑,把断层识别模型真正跑进生产流程里。
本文还有配套的精品资源,点击获取