0 引言

在地震勘探中，初至拾取要求尽可能准确地确定地震道中初至出现的时刻^[1]，是静校正、层析成像等处理的前提，初至拾取精度在很大程度上影响后续处理精度^[2]。近年来，随着“两宽一高”地震采集技术的普及，炮、检密度越来越大，采集数据达到PB数量级，并且地表向复杂山地、巨厚黄土塬等延伸，导致采集资料信噪比较低，初至不明显。如果采用人工拾取方式，工作量大，拾取效率低，容易引入系统的人为误差^[3-4]。如何快速、精确地拾取海量低信噪比数据的初至，是亟需解决的一个关键问题。

过去几十年，业界提出了许多初至拾取方法，目前使用较多的是基于地震信号瞬时特征的方法^[5]，有能量比法、相关法、图像处理法、分形维法、神经网络法等。能量比法^[6-7]就是计算初至上、下时窗的能量比，该方法简单、计算效率高，但对于低信噪比数据拾取精度较低。由于能量比法选取了初至附近的两个时窗，很难确定真正的起跳时间^[8-9]。相关法^{[4, 10]}虽然初至拾取精度较高，但相邻道互相关计算的复杂度较高，对于邻近道为干扰道或空道的情况无能为力。图像处理法^[11-12]将地震记录转化为灰度图，使用边缘检测法检测初至位置，并基于边缘追踪技术追踪初至波。该方法需要将地震记录转化为图像，拾取初至后又要将初至位置映射到地震记录，不易大规模应用。分形维法^[13-14]通过分析分形维数随时间的变化自动拾取初至，只能适用于信噪比较高、初至起跳明显且初至振幅与初至前、后的振幅有较大差异的地震记录，计算的稳定性和可靠性不高。传统的神经网络法^{[2, 15-16]}由于网络层数有限，也不易提取相邻道的空间特征信息，拾取精度不高。

近年来，以深度学习为代表的人工智能技术在计算机视觉、自然语言处理、语音识别等领域广泛使用。深度学习通过构建具有多个隐藏层的网络模型和海量数据挖掘数据深层特征^[17]，从而更好地发现数据的有效特征表示。人们也将深度学习技术应用于初至拾取中，取得了一些研究进展。Yuan等^[18]利用卷积连接方式的优势，首次将卷积神经网络(一种卷积连接和全连接方式的组合网络结构)应用于地震初至拾取。使用滑动窗方式提取时空振幅制作标签，规避了属性提取、属性选择，也极大减小了数据预处理工作量。由于仍采用全连接网络进行分类，且标签只有初至和非初至，因此网络结构不宜太深，在面对复杂地表或复杂问题时可能精度有限。刘佳楠等^[19]提出了基于全卷积神经网络拾取地震初至的方法，但该方法的训练样本尺寸较大，训练时占用的内存空间和计算资源较多，使用手工标注训练数据的标签，消耗了大量人力。丁建群等^[20]提出基于U-Net网络的地震初至拾取方法，但该方法的神经网络模型基于图像训练，无法大规模应用。王君等^[21]、帅威等^[22]和浦义涛等^[23]基于深度残差网络拾取地震初至及其软件开发，由于采用152层的深度残差网络，且对原始地震数据做线性动校正和反线性动校正处理，训练和应用效率较低。

为了克服上述方法的缺点，本文首先提出了一种快速制作初至拾取训练数据集的方法，并研发了相应的软件模块，然后提出一种基于U-Net和SegNet混合网络的地震初至自动拾取方法，最后用实际数据进行了测试。

1 基本原理

初至拾取的本质是要确定地震道上纯噪声和噪声与地震叠加信号之间的分界时刻^[24]，非常适合用图像分割算法确定该时刻，即初至位置。在深度学习中，常用的图像分割网络结构基于全卷积网络(FCN)的改进型(U-Net和SegNet)，将它们应用于地震初至拾取就是实现从地震数据到初至起跳分界时刻端到端的映射。

U-Net是一个基于卷积神经网络的图像分割网络，因其在遥感图像及医学领域的良好应用效果，获得了广泛的研究和应用^[20]。该网络结构呈U型对称(图 1)^[25]，在编码器和解码器对称的层提供跳跃连接，可以使解码过程拼接编码过程对称层的特征向量，避免了解码器直接在高级特征图中进行监督和计算损失函数值，而是结合低层特征图中的特征，使最终特征图中既包含高层特征，也包含很多低层特征，以实现不同尺度下的特征融合，从而提高网络模型图像分割精度。然而，由于上采样过程是反卷积操作，U-Net需要学习更多的参数，训练相对较慢^[26]。

SegNet被设计成一个用于像素级语义分割的核心分割引擎^[27]，解决端到端的二分类和多分类问题(图 2)。SegNet的核心在于其上采样的方法，编码器网络池化时记录最大值的位置，解码器网络反池化时将输入值直接赋给记录位置，其他位置值置零，这样不会损失图元的轮廓细节，更利于图像分割。由于SegNet将池化索引传递给上采样层，需要更少的参数，训练更快^[28]。

本文提出了一种U-Net和SegNet混合网络结构，命名为U-SegNet。U-SegNet以SegNet结构为基础，通过在解码器网络反卷积层之前融合跳跃连接信息，提供编码器网络的多尺度信息，以获得更好的性能，并且其上采样操作将U-Net中的反卷积改为反池化，池化索引被传递到上采样层，网络模型收敛更快。因此，U-SegNet网络结构更利于分割背景噪声区域和含噪信号区域，从而提高初至拾取精度。基于U-SegNet网络结构的初至拾取，可用输入炮集样本S与最终初至上、下范围分割概率矩阵M之间的一个非线性映射表示为

$ \mathit{\boldsymbol{M}} = {\rm{USNs}}\left( {\mathit{\boldsymbol{S}}, m} \right) $

(1)

式中USNs(·)表示U-SegNet网络结构，m为USNs(·)的权重或偏置参数。

2 方法流程

基于U-SegNet的初至自动拾取流程包括制作训练数据集、设计网络模型、训练网络模型、测试网络模型和实际资料应用。

2.1 制作训练数据集

为了克服以往神经网络初至拾取方法中人工标注标签耗时、费力的缺点，本文提出并编程实现了一系列在单炮记录和标签数据显示中自动联动或单独拾取训练样本及标签的方法，包括手动拾取、定位拾取、自动拾取。自动拾取方式又分为沿直线固定步长拾取(图 3)、沿曲线固定步长拾取(图 4)和矩形固定步长采样拾取(图 5)，并形成了训练数据集制作软件SeismicNetV1.0。为了方便后续网络模型的测试和应用，该软件将每个样本和标签以全局唯一标识符命名，在指定的磁盘目录下保存为占用存储空间较小的二进制文件以及JPG格式图片文件，并且将它们的文件名、文件路径、道方向长度、时间方向采样点数、中心点的道序号和时间等信息存储到数据库。

训练及测试数据来自四川三维工区A，该区地表条件复杂，低降速带岩性变化较大，很多排列的远炮检距道初至不明显。制作训练数据集过程如下。

(1) 使用某商业软件拾取炮集数据的初至，人工修正未能拾取或拾取精度不高的初至，导出文本格式的初至文件。

(2) 根据初至文件中的初至标注SEG-Y标签文件，标签文件道数和采样点数与炮集数据文件相同，其中初至时间以上样点值写为0，初至及以下样点值写为1。

(3) 利用原始炮集数据和加噪炮集数据分别制作训练数据集，并将其合并为最终的训练数据集，具体为：

1) 将原始炮集数据文件、标签文件导入SeismicNet软件中，设置样本大小为80×160(80道、160个采样点)，以矩形固定步长采样方式在不同单炮最小至最大初至范围内拾取10000个样本及对应标签，保存为10000对二进制文件。

2)首先对原始炮集数据使用

$ {A_{\rm{n}}} = {A_{\rm{r}}} + \frac{{{A_{\max }} - {A_{\min }}}}{5} $

(2)

添加噪声。式中：A_r为原始炮集数据采样点振幅；A_n为加噪炮集数据采样点振幅；A_max为以当前道为中心道的99道数据的最大振幅；A_min为以当前道为中心道的99道数据的最小振幅。然后利用拾取原始炮集训练数据的方式，从加噪炮集数据文件与标签数据文件中以矩形固定步长采样方式拾取5000对加噪训练样本及其标签，保存为5000对二进制文件(图 6)。

3) 合并两种数据源制作的训练数据集作为最终的训练数据集，训练数据集中初至样本占该区炮集数据初至的29%。

采用矩形固定步长采样方式拾取的样本千差万别，包含面波和各种随机噪声；沿理论初至曲线拾取的样本只体现初至附近的特征。因此前者拾取的训练数据集泛化能力更强，并且制作的训练数据集具有以下优点：①拾取训练样本大小为80×160，较前人采用的256×4096训练样本^[19]的网络模型训练和应用效率更高；②直接在原始炮集数据中拾取训练样本，无需对原始炮集进行各种预处理，如线性动校正、自动增益控制、高斯滤波等，制作训练数据集和测试应用较前人方法^[20-21]更高效。

2.2 设计网络模型

本文结合U-Net和SegNet网络结构，搭建了U-SegNet_5layer(图 7)和U-SegNet_4layer(图 8)网络结构。U-SegNet_5layer网络结构的编码器部分(左边)和解码器部分(右边)都包含5层。U-SegNet_4layer的网络结构与U-SegNet_5layer相似，但比后者少了最底部一层。

最终网络输出一张大小为80×160的、以0~1概率表示的分割结果图。Sigmoid激活函数为

$ {p_{v, n}} = \frac{1}{{1 + {{\rm{e}}^{ - {x_{v, n}}}}}} $

(3)

式中：p_{v, n}为第v个样本、第n个采样点的概率值；x_{v, n}为第v个样本、第n个采样点通过U-SegNet解码器的输出值。

2.3 训练网络模型

训练和测试应用网络模型的软、硬件环境为：Windows 10操作系统PC机(12核Intel酷睿i7-8700 CPU，主频为3.2GHz，内存为8GB)，深度学习平台为Tensorflow 1.14.0。

将制作的训练数据集以98:1:1的比例划分为训练集、验证集和测试集，对应的样本数分别为14700、150和150。训练过程使用小批量梯度下降法，在保证不发生内存溢出的情况下，尽可能大地设置批次大小，可以有效地覆盖训练集和测试集损失函数的尖锐极小值，从而避免损失函数极小值导致的较差泛化^[30]。在训练中批次大小设置为60，每轮迭代245次，并用验证集验证1次。

由于网络模型分割的背景噪声区域和初至及以下区域，在标签数据中分别用0和1表示，所以损失函数采用二元交叉熵函数(Binary Cross-Entropy)并将其最小化，该损失函数具有信息论解释基础、良好的数理表现和优越的性质^[31]。某个训练样本的二元交叉熵损失函数为^[26]

$ \begin{array}{l} {\rm{Cross\;Entropy}}\left( {\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Y}}, \mathit{\boldsymbol{\hat Y}}} \right){\rm{ }} = - \sum\limits_{{p_x} \in X} {\left[ {{y_{{p_x}}}\lg {{\hat y}_{{p_x}}} + } \right.} \\ \;\;\;\;\;\;\;\;\left. {\left( {1 - {y_{{p_x}}}} \right)\lg \left( {1 - {{\hat y}_{{p_x}}}} \right)} \right] \end{array} $

(4)

式中：X为训练样本；Y为X的标签；$\mathit{\boldsymbol{\hat Y}}$为网络模型预测的概率矩阵；p_x表示X中的一个采样点，y_px为p_x的标签值；${{\hat y}_{{p_x}}}$为网络模型预测的p_x的概率。

对于一个包含n个训练样本的批次，损失函数为^[26]

$ J = \frac{1}{n}\sum\limits_{i = 1}^n {{\rm{Cross\;Entropy}}({\mathit{\boldsymbol{X}}_i}, {\mathit{\boldsymbol{Y}}_i}, {{\mathit{\boldsymbol{\hat Y}}}_i})} $

(5)

式中：X_i为训练样本批次中的第i个样本；Y_i为X_i的标签；${{\mathit{\boldsymbol{\hat Y}}}_i}$为网络模型预测X_i的概率矩阵。

网络模型训练的梯度下降优化算法采用适应性矩估计(Adaptive Moment Estimation，Adam)，该算法是扩展后的随机梯度下降算法，结合了自适应梯度算法(Adaptive Gradient Algorithm，AdaGrad)和均方根传播算法(Root Mean Square Prop，RMSProp)的优点^[28]。首先计算梯度的一阶矩估计和二阶矩估计，然后为各个参数设置不同的不相关的自适应学习率^[32]，即

$ {m_t} = {\beta _1}{m_{t - 1}} + (1 - {\beta _1}){g_t} $

(6)

$ {v_t} = {\beta _2}{v_{t - 1}} + (1 - {\beta _2})g_t^2 $

(7)

$ \left\{ \begin{array}{l} {{\hat m}_t} = \frac{{{m_t}}}{{1 - \beta _1^t}}\\ {{\hat v}_t} = \frac{{{v_t}}}{{1 - \beta _2^t}} \end{array} \right. $

(8)

$ {W_{t + 1}} = {W_t} - \frac{\eta }{{\sqrt {{{\hat v}_t}} + \varepsilon }}{{\hat m}_t} $

(9)

式中：m_t和v_t分别为梯度的一阶矩估计和二阶矩估计；β₁和β₂为动力值，通常分别取0.9和0.999；${{\hat m}_t}$和${{\hat v}_t}$分别为一阶矩和二阶矩的修正值；g_t=ΔJ(W_t)为第t次迭代损失函数关于W_t的梯度，W_t为第t次迭代模型的参数；ε是为了避免分母为0而取的一个很小值，一般为10^-8；η为学习率，默认值为0.001。Adam优化算法易于实现，计算效率高，对内存的要求小，对梯度的对角重新缩放是不变的(将梯度乘以仅带正因子的对角矩阵是不变的)，并且非常适合于海量数据或多参数问题，对超参数有直观的解释，通常只需很少的调优^[33]。

训练过程中可以对某个样本在网络结构中某单元的输出进行可视化显示，可以看到不同卷积核提取输入的不同特征。图 9为某个样本及其在U-SegNet_5layer编码器网络第二层的第二个卷积层输出的64张特征图。

提前停止(Early Stopping)法可以在网络模型达到最佳状态时自动停止训练，是一种能够避免网络发生欠拟合或过拟合的方法^[32]。训练过程中设置了提前终止条件，当训练集的损失函数值在三轮迭代内出现上升即停止训练，并且训练过程中只保存当前最好的模型参数。U-SegNet_5layer迭代32轮即停止训练过程，保存了第29轮的模型参数，最终训练集准确率为99.52%，验证集准确率为99.44%；U-SegNet_4layer迭代了48轮停止训练过程，保存了第45轮的模型参数，最终训练集准确率为99.37%，验证集准确率为99.32%。图 10为U-SegNet_5layer和U-SegNet_4layer训练过程的损失函数值曲线和准确率曲线。由图可见，两个网络模型在很少的迭代次数内完成训练，说明训练效率都很高，并且U-SegNet_5layer训练和验证的准确率略高于U-SegNet_4layer。表 1为U-SegNet、U-Net、SegNet训练准确率，其中UNet_5layer使用的标签数据初至位置样点值写1，初至以上和以下样点值写0。由表可见，U-SegNet的训练准确率高于U-Net和SegNet。

2.4 测试网络模型

由于U-SegNet_5layer训练和验证的准确率略高于U-SegNet_4layer，故选择U-SegNet_5layer进行测试和应用。

首先，在制作训练数据集的炮集中随机拾取四个包含初至的80×160大小的测试样本，调用U-SegNet_5layer保存的最佳权重参数分割初至上、下范围(图 11)。可见，对于不同位置的初至，利用U-SegNet_5layer分割的初至上、下范围都很准确。

测试数据输入到网络模型中的输出只是初至上、下范围分割结果图，还需要从这些分割结果图中确定具体的初至位置。分割结果图实际上是预测为初至以上范围或初至及以下范围的概率矩阵。在分割结果图中如果某个采样点确定为初至以上范围，则概率值近于0，如果某个采样点确定为初至以下范围，则概率值近于1，初至附近上、下几个采样点的概率值为0~1。本文从上往下遍历分割结果概率矩阵的每一道，将第一个概率值大于0.5的采样点确定为该道的初至。但某些样本存在很多不包含初至的道，如图 11中第3个和第4个样本，这些道的初至位置不能由这些样本确定，而由它们上方或下方相邻的样本确定。图 12为由10张分割结果图确定的初至位置。

确定每张分割结果图的初至位置之后，需要将分割结果图的每一个道序号映射为炮集中的道序号，每一道初至位置映射为炮集中的具体初至时间，映射过程分为三步：

(1) 根据文件名从数据库中查询出某个测试样本的道方向长度、时间方向采样点数、中心点的道序号和时间。

(2) 根据以上信息和分割结果图中的初至位置，由

$ I = {I_C} - \frac{C}{2} + i $

(10)

$ T = {T_C} - \frac{{\left( {H - 1} \right)R}}{2} + sR $

(11)

计算出每一道在炮集中对应的道序号和初至时间。式中：I为炮集中的道序号；i为分割结果图中的道序号；I_C为样本中心点在炮集中的道序号；C为样本道方向长度；T为某一道在炮集中的初至时间；s为某道在分割结果图中的采样点序号；R为采样率；T_C为样本中心点在炮集中的初至时间；H为样本时间方向采样点数。

(3) 剔除可能存在的重复道及其初至。

图 13为某商业软件拾取的初至在炮集中的显示。由图可见，排列右边远炮检距道信噪比很低，初至不明显，该软件未能拾取这些道的初至，并且第三个排列右边拾取的初至存在向下跳变，拾取精度不高。图 14为在炮集数据中拾取的测试样本，图 15为U-SegNet_5layer拾取的测试样本初至在炮集中的显示。由图 15可见，U-SegNet_5layer精确拾取了图 13的远炮检距道的初至，且不存在初至跳变问题。

本文将U-SegNet_5layer用于青海M区戈壁滩实际地震资料，测试其泛化性能。由于U-SegNet_5layer是由复杂山地地表的低信噪比数据训练而成，将其应用到信噪比相对较高的M区资料，初至拾取效果很好(图 16)。

2.5 实际资料应用

在制作训练数据集过程中，选取连续250炮的887500道实际数据测试U-SegNet_5layer在大数据量炮集中的初至拾取效果。

目前几乎所有的深度学习初至拾取方法在测试应用网络模型时输入的样本包含炮集数据所有采样点，但在实际初至拾取过程中不需要输入初至以上背景噪声和初至以下深层的大部分样本，可节省大量的计算资源。为了减少输入到网络模型中的数据量，提高大数据量炮集初至拾取效率，本文采用SeismicNet软件的沿曲线自动拾取方式，在250炮数据中沿理论初至曲线以30道为步长拾取80×160大小的样本(图 17)，这样输入到网络中的数据量约为输入炮集所有采样点的1/10，应用效率提高了10倍。

在与训练网络模型相同的软、硬件环境下，U-SegNet_5layer拾取250炮实际数据用时285s，每秒拾取3114道数据，拾取效率约为某商业软件的2.2倍。

本文将网络模型的初至拾取率定义为单炮中网络模型成功拾取初至的道数占单炮总道数的比例。图 18和图 19分别为某商业软件和U-SegNet_5layer拾取250炮数据初至的拾取率曲线。由图可见：该商业软件单炮初至拾取率最高为99.7%，最低为80.8%，平均值为92.5%(图 18)；U-SegNet_5layer单炮初至拾取率最高为99.9%，最低为96.1%，平均值为98.4%(图 19)。将U-SegNet_5layer拾取250炮初至的拾取率投影到卫星遥感影像中(图 20)可见：炮点分布的西北方向靠近公路和村庄，干扰较大，初至拾取率较低；炮点分布的东南方向靠近村庄、机场和工厂，初至拾取率更低；炮点分布的中心位置距干扰源较远，初至拾取率最高。

3 结论

本文提出的基于混合网络U-SegNet的地震初至自动拾取方法克服了传统方法和近年提出的基于深度学习的方法在初至拾取效率、精度、稳定性和实用性等方面的缺点，表现为：

(1) 实现了矩形固定步长采样方式和沿理论初至曲线自动拾取训练数据集的功能，不仅省去了繁琐的训练数据预处理过程，极大地提高了制作训练数据集的效率，而且缩小了样本尺寸，减少了输入网络模型中的数据量，提高了网络模型训练、测试和应用效率。

(2) 设计的混合网络结构U-SegNet结合了U-Net可以融合不同尺度特征和SegNet可以很好地保留图元轮廓信息的优点，更精确地分割初至以上纯背景噪声部分和初至及以下部分，使初至拾取更准确。

(3) 将训练的网络模型用于实际数据测试时，提出了一种将炮集数据快速转换到初至时间的方法。应用结果表明，该方法初至拾取率和效率均高于某商业软件，易于工业化生产，具有良好的发展前景。