0 引言

地震资料解释一般可分为资料准备、剖面解释、空间解释和综合解释四个部分。其中地震反射层位追踪是剖面解释中的重要一环。传统的层位追踪方法主要依靠人工，并以三维可视化技术作为辅助手段^[1]。传统方法的优点是层位追踪精准度高，符合地质规律，但极其耗费人力和时间。于是，发展了基于种子点的层位自动追踪技术。郑公营等^[2]在主测线和联络线地震剖面上插入种子点，再根据一定的特征计算、寻找地震道之间的相似种子点，完成了层位自动追踪。常用的特征主要有波形特征^[3]和相关性^[4]。这种技术可以在一定程度上降低人力成本^[5]，但因受算法和数据的限制，在复杂地质结构且地震资料受噪声的影响时，难以达到预期效果。

随着深度学习技术的发展，越来越多的深度学习方法应用于地震资料解释领域。杨梦琼等^[6]利用知识蒸馏的回归方法追踪层位，利用教师模型引导学生模型训练，在训练参数不冗杂的情况下实现了单一层位的追踪。常德宽等^[7]提出了基于主动学习的地震多层位解释方法，利用半监督思想和大量的未标注数据实现了多层位的追踪。李海山等^[8]用深度学习方法和相对年代地质体追踪层位，精度更高。但在复杂的地质结构中，当断层多、地震反射较弱时，这些方法预测精度会较低，并且训练时间更长。

为了解决上述问题，需要深度学习网络具备像素级别处理的能力^[9]。由此，可采用语义分割的思想，首先将地震剖面切割成多个含有语义的区块，并区分每个区块的语义类别，得到像素级别的语义标注图^[10]；再利用二维卷积神经网络(CNN)能很好地学习图像特征的特点进行训练。如Wu等^[11]采用了编码―解码卷积网络进行层位追踪；Yang等^[12]采用了CNN的层位自动追踪；朱振宇等^[13]采用了U-Net解释了地震反射层位。然而这些二维卷积网络的感受野有限^[14]，需要加入多尺度特征融合模块或空洞卷积才可以更好地利用上、下文特征去提高语义分割模型的准确性^[15-17]，但这仍会损失一部分细节特征，导致在用于非全区分布层位或者过断层层位时预测效果不佳。

为了选择合适的网络模型，需要详细了解语义分割和地震反射层位追踪的特点。经典的语义分割算法是：首先勾勒出图像中各个物体的准确轮廓，对图像中每个像素点进行判断，并划分一个对象类别^[18]；再将同一对象类型规定为同一个类标签。与常规图像处理方法不同，地震反射层位追踪难点是：①需要区分的类别较多。一般1幅图片仅包含2~4个需要区分的主体，而1幅地震剖面追踪的层位可达10个以上。②训练数据较少。由于层位标注困难，部分地震数据体中只有少数剖面可用于训练。③在图像中层位位于相邻主体区域的边缘处时，要求算法对边缘处的计算精度高。④分割的区域较大，判断单个像素点的分类时，需要参考更多区域数据的上、下文特征，也就是需要更大的感受野。

U²-Net是一种U型结构(双层U型嵌套)网络^[19]，设计了多个ReSidual U-blocks(RSU)模块，每个模块中的感受野不同，使U²-Net网络能够融合不同大小感受野，充分利用上、下文信息提升预测精度，能够找回在下采样提取特征时所损失的一部分边缘特征，数据量较少时也不容易过拟合。并且每个RSU都是U型结构，大部分计算都在下采样操作中，在不加大计算成本的情况下加深了网络深度。

综上所述，本文选用U²-Net，利用多尺度特征和尽量少的标签量提高预测精度。

深度学习的训练中标签至关重要。在Wu等^[11]的基于编码―解码卷积网络进行层位追踪和Yang等^[12]的基于CNN的层位自动追踪中，都采用了非充填标签，因背景占比过大，在预测时极易出现跳点和连续性不好的情况。朱振宇等^[13]在基于U-Net的地震反射层位追踪中使用了充填标签，但是充填标签的制作方法没有涵盖非全区分布的层位或过断层的层位。

为此，本文采用全新的充填标签制作方法，即遍历地震数据每个像素点，判断当前像素点所在位置并为其划分一个层位区域。对于穿过断层的层位，则自动搜寻相邻层位，实现非全区分布层位、断层等复杂条件下的地震反射层位及不整合面的充填标签的制作。最后，采用U²-Net网络模型对F3数据体和M工区地震资料进行训练，以测试本文方法的效果和适用性。

1 方法原理 1.1 网络模型

与U-Net不同，U²-Net采用了RSU模块。

1.1.1 RSU模块

RSU模块的作用主要是捕获数据的多尺度特征，根据深度的不同可以进一步细分为RSU-7、RSU-6、RSU-5和RSU-4模块。以RSU-7 (图 1)为例，RSU模块由三部分组成：①卷积输入层。包含1个Conv+BN+ReLU层。②U型的编码―解码结构。左侧为编码结构，由L-2个(L为当前RSU模块的深度)上采样×2+Conv+BN+ReLU层和1个Conv+BN+ReLU层组成；右侧为解码结构，由L-2个下采样×1/2+Conv+BN+ReLU层和1个Conv+BN+ReLU层组成。③特征拼接层。底部卷积采用膨胀卷积，膨胀率为2。

在RSU模块中，首先，在输入层将输入数据$ \boldsymbol{X}(\mathrm{尺}\mathrm{寸}\mathrm{为}H\times W\times {C}_{\mathrm{i}\mathrm{n}} $)转化为$ \mathrm{尺}\mathrm{寸} $$ W\times M\times {C}_{\mathrm{o}\mathrm{u}\mathrm{t}} $的中间特征图，提取局部特征$ F1\left(\boldsymbol{X}\right) $。然后，在编码阶段以$ W\times M\times {C}_{\mathrm{o}\mathrm{u}\mathrm{t}} $中间特征图作为输入，逐步下采样将数据尺寸减半、通道数翻倍以增大感受野。在解码阶段，先与左侧同一深度的输出数据进行级联合并，再逐步上采样。据此可以减少下采样过程中造成的细节损失，得到L深度下的尺度特征$ U\left[F1\left(\boldsymbol{X}\right)\right] $。RSU-L模块深度越大，其内部包含的池化操作越多，经由RSU-L模块训练的数据的感受野也就越大，局部特征也就越丰富。同样，越浅层的RSU模块，全局特征越丰富。RSU模块可以根据L的不同从输入特征图中提取多尺度特征。最后，在特征拼接层，将局部特征$ F1\left(\boldsymbol{X}\right) $和深度L下的尺度特征$ U\left[F1\left(\boldsymbol{X}\right)\right]\mathrm{进}\mathrm{行}\mathrm{拼}\mathrm{接} $，可得到此模块下输出的特征数据$ F1\left(\boldsymbol{X}\right)+U\left[F1\right(\boldsymbol{X}\left)\right] $。

除上述模块外，还有一种特殊的RSU模块，即RSU-4F(图 2)。相较于RSU-4，RSU-4F使用了空洞卷积替换了原本的上、下池化操作。这是因为在编码四层之后特征图的分辨率已经很低，如果进一步降低可能会导致重要的上、下文信息丢失。

1.1.2 外部结构

U²-Net结构如图 3所示，从左往右依次为编码器、解码器和显著特征融合模块。

在编码阶段，En-1~En-4分别采用RSU7~RSU-4模块，而En-5和En-6均采用RSU4F模块。每个RSU模块之间插入一个核尺寸为2×2的下采样，特征图尺寸逐渐变为原来的1/32。

在解码阶段，De-1~De-4分别采用RSU7~RSU4模块，De-5采用RSU4F模块。每个RSU模块之间插入一个核尺寸为2 ×2的上采样，并且同时将上一层上采样的输出与来自对称编码器阶段的输出的级联作为输入，特征图尺寸逐渐恢复为原来的大小。

在显著特征融合阶段，在En-6、De-5~De-1的右侧均插入一个Conv+Sigmoid的卷积+激活函数层，得到通道数为1的特征图，并将其分别上采样至输入图，得到$ {S}_{\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}}^{\left(1\right)} $~$ {S}_{\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}}^{\left(6\right)} $；然后，经过1×1的卷积和Sigmiod激活函数，得到所需的预测概率图Sup1~Sup6；再将它们进行拼接融合后得到$ {S}_{\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{e}} $；最后，经过卷积和Sigmoid函数得到最终概率图Sup0。

1.2 标签制作

将人工标注的层位坐标作为标准去标注多套层位，并将每一层赋予不同的标签数值1，2，…，得到如图 4所示的语义分割训练数据标签。

如果直接用图 4中的标注数据作为标签进行训练，将会因为前景与背景数据分布特别不均衡而极易出现过拟合的现象，导致模型训练失败。为了解决这一问题，本文采用了充填标签生成策略，即将两个相邻层位之间的部分标注为1个层位区域，使其下边界同属于该层位区域。据此，可使原5个层位和顶底平层(图 4)变为8个层位区域(图 5)。此时，8个区域(层位标签)占比相对均衡，两层之间的分界也很明显。

充填标签制作方法如下。

(1) 遍历地震数据。如图 6a所示，首先将三维地震数据体V(x, y, z)沿着线号(x)拆分为众多二维垂直地震剖面(y, z)，然后在每一条剖面中根据道号(y)垂直拆分为众多一维地震道(z)。遍历每一个地震道，即可遍历每条地震剖面，进而遍历整个地震数据体。

(2) 充填上边界为全区分布层位的层位区域。如图 6b所示，对于每一个地震道而言，根据如下规则确定地震道上数据点属于的层位区域，即

$ \begin{array}{l}\mathrm{I}\mathrm{f}\mathrm{ }\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\left(\mathrm{c}\mathrm{m}\mathrm{p}1;h1\right) < Z < \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\left(\mathrm{c}\mathrm{m}\mathrm{p}1;h2\right);\\ \mathrm{T}\mathrm{h}\mathrm{e}\mathrm{n}S\left(x, y, z\right)==h2\end{array} $

(1)

式中：$ z $为当前地震数据点的垂向坐标值；$ S\left(x, y, z\right) $表示当前数据点S的线号、道号和时间。

然后，赋予数据点和该区域下边界层位相同的标签数值1，2，…。同时，该规则明确了充填层位区域含完整的下边界，这适用于层位非全区分布或断层发育的情况。

(3) 充填上边界为非全区分布层位或有断层发育的层位区域。当上边界为非全区分布层位或断层时，上边界的$ \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}(\mathrm{c}\mathrm{m}\mathrm{p}1;h1) $采样时间在层位断开或不存在处为nan值，无法用规则进行判断，此时下边界不完整。对于还没有充填的下边界部分，上边界$ \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}(\mathrm{c}\mathrm{m}\mathrm{p}1;h1) $选用距原上边界最近、采样时间最少的层位，再利用上述规则进行判断。

对于穿过断层的层位，因其缺少下边界约束，利用判断规则划分层位区域时，断层线区域不能填充为层位类而被划分为背景类，如图 7右侧存在两条断层。

1.3 U²-Net损失函数

本文采用有辅助分类器的深度监督损失函数^[20]

$ \mathrm{L}\mathrm{o}\mathrm{s}\mathrm{s}=\sum\limits_{m=1}^{M}{w}_\text{side}^{\left(m\right)}{l}_\text{side}^{\left(m\right)}+{w}_\text{fuse}{l}_\text{fuse} $

(2)

式中：M是编码或解码阶段数，当M = 6时，符合本文网络模型需求；$ {l}_{\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}}^{\left(m\right)} $为概率图Supm与真实值之间的损失值；$ {l}_{\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{e}} $为最终概率图Sup0与真实值之间的损失值；$ {w}_{\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}}^{\left(m\right)} $和$ {w}_{\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{e}} $是每个损失项的权重，本文网络模型中全部设置为1。对于每一个损失值$ l $，使用标准二进制交叉熵计算损失，具体定义为

$ l=-\sum\limits_{\left(r, c\right)}^{\left(H, W\right)}\left\{{P}_{G\left(r, c\right)}\mathrm{l}\mathrm{g}{P}_{S\left(r, c\right)}+\left(1-{P}_{G\left(r, c\right)}\right)\times \right.\\ \left.\mathrm{l}\mathrm{g}(1-{P}_{S\left(r, c\right)})\right\} $

(3)

式中：$ (r, c) $为概率图的像素点坐标，其中$ r、c\mathrm{分}\mathrm{别}\mathrm{为}\mathrm{其}\mathrm{横}\mathrm{坐}\mathrm{标}\mathrm{、}\mathrm{纵}\mathrm{坐}\mathrm{标}\mathrm{值} $；$ {P}_{G(r, c)} $为标签图中$ \left(r, c\right) $点的像素值；$ {P}_{S(r, c)} $为预测的概率图中$ \left(r, c\right) $点像素值。

1.4 评价指标

本文采用准确率(Pixel Accuracy，PA)、类准确率(Class Accuracy，CA)、交并比((Intersection over Union，IoU)、平均交并比(Mean Intersection over Union，mIoU)评价本文方法效果。

(1) PA计算概率图中分类正确的像素点占所有像素点的比例，用于评估模型整体效果。

$ \mathrm{P}\mathrm{A}=\frac{\sum\limits_{i=1}^{k}\left|{F}_{i}\bigcap {G}_{i}\right|}{\sum\limits_{i=1}^{k}\left|{G}_{i}\right|} $

(4)

式中：$ {G}_{i} $表示在地震剖面上属于第i类层位区域的像素点集合；$ {F}_{i} $表示被预测为第i类层位区域的像素点集合；$ k $为总类数。

(2) CA计算第i类层位区域在概率图中的准确率，用于评估模型对单个层位的效果。

$ \mathrm{C}\mathrm{A}=\frac{\left|{F}_{i}\bigcap {G}_{i}\right|}{\left|{G}_{i}\right|} $

(5)

(3) IoU计算第i类层位区域的概率图结果和标签的交集并集的比率，用于评估全局预测效果。

$ \mathrm{I}\mathrm{o}\mathrm{U}=\frac{\left|{F}_{i}\bigcap {G}_{i}\right|}{\left|{F}_{i}\bigcup {G}_{i}\right|} $

(6)

(4) mIoU计算整个概率图结果和标签的交集并集之比。

$ \mathrm{m}\mathrm{I}\mathrm{o}\mathrm{U}=\frac{1}{k}\underset{}{\sum\limits_{i=1}^{k}}\frac{\left|{F}_{i}\bigcap {G}_{i}\right|}{\left|{F}_{i}\bigcup {G}_{i}\right|} $

(7)

PA和CA都是基于概率图中单个层位区域分割效果的评价。当层位区域类别分布不均匀、小类误差严重时，它们无法反映全局的分割效果。网络模型除了将i层位区域本身预测为i层位区域外，还有可能将大量其他层位区域预测为i层位区域，预测结果就会出现不连续、跳点等情况，这种预测错误无法仅根据PA和CA评价。因此，结合IoU指标可以从全局角度评估复杂的层位预测效果。

除上述量化指标之外，仍要结合地质认识评价层位预测结果，即沿着预测的层位提取地震属性，再通过地震属性与已知的地质认识进行对比分析。如果二者一致，就可以认为预测的层位结果合理。

2 实际资料应用

将本文方法应用于F3地震数据体以测试层位预测效果，并利用M工区地震数据进行验证。F3地震数据体包括200~679线共480条剖面、400~1199道共800道、597个采样点。从上到下层位依次为horFS8、horTOP、horFS6，horMSF4、horFS4。其中，horFS6在每条测线上局部分布，horTOP在部分测线上局部分布(图 8)，这可测试网络预测层位的全局能力和细节能力。另外，剖面两侧地层较为平缓，右侧1/3处发育一条断层，其附近地质结构复杂，这可用于测试网络的准确性、稳定性。

2.1 数据预处理

在测试之前，需要对数据预处理以便于模型提取数据特征。

2.1.1 归一化

采用MinMaxScaler方法对F3地震数据进行归一化处理，消除数据特征大小的差异性，提高模型收敛速度，提高模型预测精度。

2.1.2 数据增强

采用重叠步长窗口的方式进行数据增强，目的是扩充样本多样性，避免过拟合，但不改变地下地层的结构。同时，窗口的重叠避免了部分数据处于边界状态而降低预测精度的缺陷。

对于F3地震数据体，选取尺寸为512×512的窗口，设置重叠步长为256，从上到下、从左往右依次选取。当窗口到数据体边界时，则将窗口内缺少的数据用0补齐。分块后标签结果如图 8所示。

2.1.3 独热编码

每个层位类标签都是独立数值。为了增强标签数值的稀疏性，减少深度学习预测结果的错误，提升泛化能力，本文将标签数据值转化成独热编码(One-Hot)。即用N位状态寄存器对N个标签值进行编码，使每个独立的标签值都有独立的寄存器位。实现方法是：对于地震剖面上第n个像素点，其标签值c用1个0与1组成的向量$ {\boldsymbol{T}}_{n} $表示，即$ {\boldsymbol{T}}_{n}=[{T}_{1n}, {T}_{2n, }\dots , {T}_{cn}, \dots , {T}_{kn}] $。若层位区域标签为c，则相应的$ {T}_{cn} $为1，具体对应关系如表 1所示。

2.2 测试过程

每间隔50线取1线作为训练集，其余线作为验证集。采用Adam优化器，学习率设置为0.001，训练轮数为200，使用GPU、3060显卡进行训练。为了评价U²-Net网络效果，将U²-Net模型与语义分割模型中认可度较高的U-Net+PPM模型进行对比分析。训练完成后记录训练时间，其中U²-Net模型每轮训练时间为27 s，训练200轮所需时间为90 min；而U-Net+PPM模型每轮训练时间为53 s，训练200轮所需时间为176.4 min，这说明了U²-Net训练效率更高。

2.3 测试结果

为进一步分析多套层位追踪效果，选取F3地震数据体中具有代表性的3条线(204线、563线、475线)和M实际工区的1条线(7450线)进行测试。其中，F3地震数据体的3条线用于测试方法预测结果的精度和稳定性、模型的泛化性能，M工区的1条线用于测试本文方法的适用性。

2.3.1 204线层位预测结果分析

204线靠近标签线，U²-Net模型(图 9)与U-Net+PPM模型(图 10)的预测精度均较高。其中，horFS6的层位区域数据量较少，且为非全区分布层位，同相轴尖灭处预测难度比较大，可用于分析网络模型的精度。

由表 2可见，U²-Net模型的PA和mIoU均较高，而U-Net+PPM模型的PA较高，但mIoU较低。观察U-Net+PPM模型各层位IoU可以发现，horFS6的IoU非常低，这说明U-Net+PPM模型在预测难度大时精度较差。

由图 9、图 10可见，U²-Net模型的预测误差比U-Net+PPM模型的小，预测层位结果更接近人工标注，且不易出现跳点、连续性差等问题。

对于所处同相轴振幅变化速度不快的层位(如horFS8、horMSF4、horFS4)，U²-Net模型和U-Net+PPM模型层位预测效果均很好。但在经过断裂等地质结构复杂区时，U-Net+PPM模型的追踪结果存在跳点(图 10c)，而U²-Net模型的预测层位连续性更好(图 9c)。

对于所处同相轴振幅变化速度较快的层位(如horTOP)，U-Net+PPM模型所预测的层位振幅变化较大处与标签层位(图 10c)不吻合，而U²-Net模型预测结果中二者吻合较好(图 9c)。这是因为U²-Net模型拥有RSU模块，RSU-7~RSU-4利用不同的模块深度，使U²-Net模型在编码和解码过程中均能多尺度利用数据的上、下文信息。

对于非全区分布层位(如horFS6)，U-Net+PPM模型预测结果存在部分CA比较高(97.58%)，但IoU仅为59.49%(表 2)；在非全区分布层位尖灭处出现跳跃现象(图 10c)。而U²-Net模型的预测结果(图 9c)不仅CA为98.69%，而且IoU为93.3%，非全区分布层位尖灭处预测效果好于U-Net+PPM。

综上所述，U²-Net模型的预测精度相较于U-Net+PPM模型有了很大的提升。

2.3.2 363线结果分析

363线(图 11~图 12)中horTOP和horFS8均为非全区发育层位，且两个非全区发育层位相邻，同相轴追踪难度较大，可用于分析模型的稳定性。

由表 3可见，U²-Net模型与U-Net+PPM模型的PA和mIoU均较高，但U-Net+PPM模型horFS6和horMFS4的IoU相对偏低，这说明这两层的预测结果存在明显波动。

从图 11、图 12可以看出，U²-Net模型的稳定性比U-Net+PPM模型有较大的提升。U-Net+PPM模型预测的horTOP层位和标签层位吻合性差；预测的horFS6层位有较多的跳点(图 12c)，挤占了horMFS4的层位区域，导致horMFS4的IoU偏低。而U²-Net模型预测的horTOP层位和标签层位吻合更好，预测的horFS6仅在同相轴尖灭处有少量上升的跳点。

2.3.3 475线结果分析

F3地震数据体起始线号为200，每间隔50线选取了1线作为训练集，因此从225线开始的每50线(如475线)都距离标签剖面最远。选择475线可测试网络模型预测精度和泛化性。

由图 13可见，U-Net+PPM模型的预测层位的跳点更多，右侧边缘处不够准确(图 13a)；而U²-Net模型的预测层位的跳点相对较少，预测效果更好。由此可证实同一数据体U²-Net模型的泛化性更好，并且在标签占比仅为1.875%的情况下，也能取得较好的预测结果，大大减少了人工标注层位的难度。

值得注意的是，在不同地质背景的资料中已经证实网络迁移性较差，例如本文F3地震数据体训练的模型不能预测其他工区的层位。而在地质背景类似的工区，选取一定差异性数据测试，训练出的网络模型具有一定的迁移性。

2.4 有效性验证

为了进一步验证本方法的有效性，选用M工区地震数据体进行测试。该数据体左侧资料信噪比较低，右侧资料信噪比较高(图 14a)。解释的目标层位处于1230~2000 ms，从上至下共有5个目标层位，其中horT5和horT1可见2个小断层。数据体共1001条剖面，选取14条剖面作为训练集，其余剖面作为验证集，其他参数设置与前文一致。

取验证集中7440线结果分析。由图 14可见，无论是在低信噪比资料处还是高信噪比资料处，U²-Net模型预测结果都比较合理，对于切割horT5和horT1层断层，预测层位均能与地质规律一致。

3 结束语

本文利用一种全新的充填标签生成策略，弥补了传统标签在处理不整合面时的缺陷，能满足更加多样的层位解释。结果表明，本文方法预测结果在精度、稳定性、泛化性等方面均显著提高，这为多套层位追踪提供了一种新的高精度层位解释方法，并且可以适应不同信噪比地震资料的反射层位预测。但是，限于标签样本的数量，深度学习技术需要结合地质规律对预测不合理区域进行人工修正，经过校正后的层位才能应用于实际生产。