张 岩 李 杰 王 斌 李新月 董宏丽

(①东北石油大学计算机与信息技术学院，黑龙江大庆 163318；②东北石油大学人工智能能源研究院，黑龙江大庆 163318；③黑龙江省网络化与智能控制重点实验室，黑龙江大庆 163318)

0 引言

在理想情况下，地震数据的采样规则且密集。但是，受复杂地表、地下地质条件、采集因素等的影响，易出现空间采样不足、地震数据缺失等不规则的问题，影响后续地震资料的处理、解释结果，因此地震缺失数据的规则化重建非常重要。

从二十世纪八十年代起，学者们就开始研究地震数据规则化重建方法[1]。目前常用的方法大致分为两类：一类是基于传统模型；另一类是基于深度学习。

基于传统模型的地震数据规则化方法主要分为五种，即：①基于相干倾角插值的方法[2-4]，通过扫描时空窗内同相轴的倾角，然后沿着若干个倾角方向通过加权和的形式产生内插的地震道。此类方法处理过程复杂且易受噪声干扰，难以应用到实际场景中。②基于变换域的方法，利用地震数据在某个变换域的稀疏性进行重建[5-10]，如利用Fourier变换[6]、Con-tourlet变换[7]、Ridgelet变换[8]、Curvelet变换[9]、快速离散曲波变换(Fast Discrete Curvelet Transform，FDCT)算法[10]等，可以较好地进行地震数据规则化。③基于波场延拓算子的方法[11-13]，把缺失道作为零道，结合波动方程部分偏移对叠前地震数据进行重建。该类方法可将倾角时差处理(Dip-Moveout Processing，DMO)与反DMO相结合实现地震数据重建，但在地质信息缺失的情况下，存在重建精度低和运算量大的缺点。④基于滤波器的方法[14-16]，通过褶积插值滤波器实现不规则数据重建，通常把某类插值滤波器应用于待规则化的数据，进行褶积操作。此类方法往往具有较高的计算复杂度。⑤基于压缩感知(Compressed Sensing，CS)[17-21]的规则化技术，将不规则地震数据作为完整地震数据的少量信号投影值，在处理端通过稀疏性约束正则化方法实现数据的近似重建，从而突破奈奎斯特采样定理的瓶颈。上述方法共同的问题是所建立的复杂模型通常求解困难，且仅适用于某个特定情况，模型的泛化能力较差。此外，部分模型在缺少地震数据先验知识(如偏移速度、均方根速度、叠加速度等)的情况下，地震数据规则化处理的质量会受到影响。

近年来，深度学习理论与方法逐渐用于地震数据重建。基本原理是通过学习大量的地震数据样本，得到目标区块地震数据分布特征的非线性映射函数，预测相应位置上缺失地震道的实际值，以达到恢复数据中缺失道、规则化地震数据的目的。常见的深度学习重建规则化地震数据的方法包括卷积神经网络(Convolutional Neural Networks, CNN)、生成对抗神经网络(Generative Adversarial Networks，GAN)和自编码器(Auto-Encoder, AE)三类。其中：①基于CNN的规则化方法中，Wang等[22]提出了一种基于ResNets网络的地震资料插值算法，在网络层数较深时效果较好，但需要事先对缺失的数据做预插值处理；Wang等[23]提出利用格林函数的空间互易性重建数据，降低了神经网络在训练过程中对样本多样性的依赖，但同样需要对网络输入的地震数据缺失部分做预插值；高静怀等[24]通过交替迭代求解地震数据重建的最小二乘法问题和预训练的网络模型对地震数据进行重建，同样需要在迭代初始以及输入网络前对缺失数据进行预处理。王钰清等[25]提出了一种基于数据生成和增广的CNN，可用于小样本的网络训练。②基于GAN的规则化方法中，Chang等[26]等提出了一种基于GAN的地震数据插值重建技术；Oliveira等[27]利用GAN对叠后地震数据进行重建。这两种方法均取得了较好的效果，但GAN网络主要存在的问题是训练过程不稳定，结果收敛困难。③基于AE的规则化方法中，郑浩等[28]利用卷积AE，学习完全采样地震数据与缺失重建数据的映射关系，通过残差学习预测缺失数据并进行重建输出，在测试模型上取得了较好的效果；Jia等[29]利用支持向量机重建规则欠采样的地震数据，需要对网络模型进行预插值；宋辉等[30]提出了一种基于卷积降噪AE，可以对地震数据以无监督的方式去噪。

上述深度学习规则化地震数据重建方法利用样本数据时域均方误差(Mean Square Error，MSE)信息损失作为约束，通过网络输出的规则化数据逼近实际完全采样的地震数据，可以达到较好的信噪比评价结果。Zhu等[31]考虑频域特征的提取，利用短时Fourier变换将时域的数据转化到频域，将实部和虚部传入CNN，通过逆变换得到时域的重建地震数据。此方法可以在频域上消除混叠效应，但在能量较弱区域效果不理想。基于深度学习的方法不需要建立复杂的数学模型，相对于传统的基于模型的方法，它能得到数据深层的特征信息，在缺少地下介质先验知识的条件下也可以取得较好的效果。但是，目前基于深度学习的方法还存在一些问题：①通常只关注地震数据单一域特征信息的提取，未挖掘数据联合域的特征信息；②在时域上重建地震数据的方法容易出现细节模糊或过于平滑的现象，丢失纹理信息，影响后续地震数据的解释；反之，若仅关注频域的特征信息，在地震数据能量较弱的区域重建的数据质量较差，无法反映实际的数据特征。

小波具有多尺度分析、多方向性的特点，在图像恢复领域取得了较好的效果。Anbarjafari等[32]提出了利用小波变换将图像从频域分解为平滑子图与细节子图的方法，并利用细节子图辅助恢复高质量的图像；Gao等[33]提出了一种混合小波的卷积网络，将输入数据分解成稀疏特征图谱，然后用另一个卷积网络进行稀疏编码、恢复图像；张岩等[34]考虑其他物理约束信息，提出了一种联合傅里叶域的去噪CNN，取得了较好的纹理保持效果和较高的信噪比；Wang等[35]引入阻抗域的损失，提出了一种基于闭环CNN测井约束地震反演方法，应用于真实地震数据，获得了较好的效果。

据此，本文提出基于联合小波域深度学习的地震数据规则化方法。首先，建立联合深度CNN学习地震数据在时域与小波域的分布特征，得到规则化数据的预测模型，将不规则地震数据的重建转化为CNN小波系数的预测；其次，设置结合时域与小波域的损失函数，通过地震数据的整体特征和局部细节信息约束网络模型，调节联合损失函数的平衡系数以调整网络模型学习的注意力。模型测试和实际数据应用均表明，该方法具有很好的细节保持效果，更具鲁棒性。

1 方法原理

1.1 不规则地震数据重建模型

假设完整的地震数据为x。实际上，在缺道或稀疏采样等条件下，采集到的不规则地震数据y可以表示为

y=Rx

(1)

式中R为采样矩阵，表示从M道地震数据中采样得到M1道(M>M1)。基于深度学习的地震数据规则化重建过程，就是通过学习样本特征从y中重建得到完全采样地震数据的近似估计x′(x′≈x)，利用多次迭代使训练结果趋于稳定，并逐渐逼近x。

1.2 基于小波变换的地震数据规则化

以样本x作为时域的标签。通过滤波器组对x进行二维小波分解，即

(2)

式中：hψ为低通滤波器；hφ为高通滤波器；下标中的“-”表示卷积操作中的翻转操作；⊕表示卷积运算；j=0,1,2,…,J-1，J为最大尺度因子；m=n=0,1,2,…,2j-1；A、V、H和D分别表示低频分量、水平高频分量、垂直高频分量和对角高频分量，对应小波域的标签C=(C1,C2,C3,C4)。

以y作为网络输入，利用网络G训练得到各频率分量的小波系数

(3)

式中：GA、GV、GH、GD分别为低频分量、水平高频分量、垂直高频分量、对角高频分量对应的映射；A′、V′、H′、D′分别表示为网络训练所得小波系数不同方向低频分量、水平高频分量、垂直高频分量、对角高频分量，对应网络训练的预测小波系数C′=(C′1,C′2,C′3,C′4)。

小波反变换后得到重建时域数据x′。联合计算时域与小波域的误差，设置全局损失函数ltotal作为约束，并通过正向传递和反向传播调整网络参数。经过多次迭代，网络训练结果趋于稳定，得到最终网络模型，将地震数据规则化问题转化为小波系数预测问题。

损失函数采用在数据采样点处x与x′的均方误差，即

(4)

1.3 小波域特征提取

地震数据的波前信息在时域上表现为复杂的纹理状曲线，尽管多层CNN具有较强的特征提取能力，但仅利用时域信息提取特征时具有很大的局限性。小波变换通过缩、放母小波的宽度获取信号的频率特征，平移小波基获取信号的时间信息[36]。信号的小波变换相当于利用母小波的缩、放和小波基的平移，并与原始信号卷积，得到小波系数(图1)。

图1 小波变换示意图

选择haar小波作为小波基，根据二维快速小波变换(Fast Wavelet Transform，FWT)[37]计算haar小波系数，图2为地震数据进行小波分解的实例。由小波变换的系数分布特性可知，低频小波系数的预测可以保留地震数据全局的主要特征信息，高频小波系数的预测有利于地震数据细节的恢复。

图2 时域(a)与小波域(b)地震数据对比

为了充分利用小波域的特征，本文设计了小波预测损失和纹理损失两种损失函数。前者是小波域均方误差的加权形式，定义为

(5)

式中：λk是平衡不同小波子带重要性的权重系数；n为小波系数分量数。赋予高频系数较大的权重，可以将训练注意力集中在局部高频细节重建上，生成细节效果较好的规则化地震数据。为了防止过拟合导致高频小波系数收敛到0，定义纹理损失为

(6)

式中：q为约束高频小波系数的初始子带序号，本文采用一级小波分解，生成的低频分量子带保留数据主要信息分量，其余三个子带为不同方向的高频分量；γk、α和ε均为平衡系数，其中γk为各高频分量的权重，α略大于1，ε略大于0。ε保证了ltexture不为0，从而使高频小波系数非零，防止小波高频系数训练过拟合。

1.4 联合小波域深度学习网络设计

本文设计的联合学习CNN模型G如图3 所示，由嵌入层网络、小波系数预测网络和联合损失计算3个子网组成。

图3 网络模型结构

1.4.1 嵌入层网络

不规则的地震数据输入嵌入层网络后，经过多层卷积特征提取，将得到的特征图谱传递给后续网络。输入的不规则地震数据尺寸为128×128，嵌入层网络的所有卷积核尺寸均为3×3，步长为1。通过补“0”操作使特征图谱尺寸与输入数据相同。卷积后的特征图谱经过归一化和激活函数，再进入下一层操作。卷积层、归一化层和激活函数构成一个残差块，前、后残差块之间设置残差连接，可以加速收敛和防止梯度消失。每层卷积核的数量沿前向递增，分别是128、256、512、1024，为小波系数预测网络提供足够的特征信息。

1.4.2 小波系数预测网络

小波系数预测网络由4个独立的并行子网络组成，它利用嵌入层网络提取的特征图谱学习、预测不同频率分量的小波系数。地震数据具有较强的曲线纹理特征。以一级小波分解为例，将小波预测网络分成4个独立的小波预测子网络；再通过CNN单独处理对应方向子带的小波系数，重建规则化的地震数据。每层网络的卷积核的尺寸为3×3，步长为1，网络预测的小波系数尺寸与不规则地震数据输入的尺寸相同。由于每个子网预测的小波系数相对独立，因此网络更具鲁棒性。

1.4.3 联合损失计算

联合损失计算是对地震数据小波域和时域的损失进行评价。网络将各个独立小波预测子网的小波系数反变换为时域的规则化地震数据，使其与完全采样的地震数据进行对比、计算损失，再反向传递到网络中，从而更新网络参数权重。

1.5 联合损失函数

网络的联合损失函数由三部分组成，即全局信息预测损失、小波系数预测损失和纹理细节预测损失。其中，全局信息损失是计算重建的时域规则化地震数据与完全采样地震数据之间的损失，以均方误差作为约束条件；小波系数预测损失是利用完全采样地震数据的小波系数对预测子网得到的小波系数施加约束；纹理细节预测损失是计算输入数据与标签之间小波变换高频分量的误差，加强此部分的约束可改善地震数据规则化的纹理细节效果。

设空间域的均方误差为全局损失lfull，联合误差函数为

ltotal=lwavelet+μltexture+νlfull

(7)

式中μ和ν是平衡参数，可根据网络训练目标的侧重点进行调整，以改变网络学习的注意力。

由于小波系数与不规则的地震数据输入尺寸相同，因此通过CNN可使每个特征图谱的大小保持一致，可以降低训练难度、兼顾地震数据的全局拓扑信息和纹理细节信息。

2 模型测试

2.1 参数设置

选择Marmousi模型测试本文方法。震源和检波器置于地表，通过中间放炮、两端接收的方式进行正演获取地震道集数据，检波器采样间隔为4ms，道间距为10m。完全采样实验样本数据裁剪尺寸为128×128的切片数据x，作为训练样本时域标签。从完整地震数据抽取比例为r的地震道作为不规则的地震数据样本，其他道为空道，分别以随机抽取和均匀抽取的方式仿真稀疏采样两种不规则情况。

地震数据重建效果的衡量指标采用信噪比(Signal-to-Noise Ratio，SNR)

(8)

式中N表示样本数量。

将10000组Marmousi模型数据按照8∶1∶1比例且不交叠的方式分别划分为训练集、验证集和测试集。在训练阶段，使用训练集数据进行训练并用验证集评估网络训练效果；当验证集结果趋近收敛时保存网络模型和参数，并用测试集数据对预训练模型进行评估。

2.2 网络模型测试

针对均匀采样和随机采样两类情况，制作从10%～90%共9组不同采样率下的样本(每组采样率间隔为10%，两类采样样本共18组)，分批次输入网络进行训练。不同采样率下训练的不同网络模型用来测试相应采样率下的地震数据。

由图4可见，在不同采样率条件下，本文方法(联合小波域的深度学习方法)均较好地保留了实际地震数据的特征；随着采样率不断提高，规则化效果也随之提升。

图4 不同采样率条件下均匀采样数据(上)及其规则化重建结果(下)

图5为本文方法在验证集测试的评估结果。由图5a可见，随着迭代次数增多，SNR逐步提高，重建效果不断增强；当达到1000次时，算法取得收敛，SNR基本稳定在最大值附近，说明本网络模型具有较好的收敛性和稳定性。为了使网络更好地收敛至最优解，防止训练后期学习率过大，本文采用指数衰减算法，设置初始学习率为2×10-4，每次迭代衰减指数为0.995(图5b)。迭代初期用较大的学习率使结果快速收敛，所以前期联合误差数值变化相对较为剧烈(图5c)，有助于加速重建地震数据；后期减小学习率，使目标函数收敛至局部最小值，联合误差变化相对缓慢(图5c)，有利于重建规则数据细节信息。另外，在训练过程中网络调整参数导致结果振荡，但是网络会根据损失函数进行调整，从而使最终结果趋于稳定。

图5 SNR(a)、学习率(b)和ltotal(c)随迭代次数变化趋势

2.3 纹理细节保持效果

为了验证本文方法对细节重建的有效性，测试样本在50%均匀采样条件下，对比联合小波域学习网络(SNR=19.8550dB)和仅使用时域学习网络(SNR=17.2667dB)的结果(图6)。由图可见，联合小波域学习的网络局部细节特征更准确，更接近于真实地震数据。这证明了联合小波变换的CNN具有更好的纹理保持性能。

图6 50%均匀采样条件下全采样(a)与未采用联合小波(b)、采用联合小波(c)的重建波形(上)及其局部(红框)放大(下)

2.4 算法对比测试

将本文方法应用于地震数据随机缺失的情况，并与当前较先进的重建算法进行对比(图7)。由图可见，本文方法重建的地震数据波形连续性好、无突变，更逼近真实地震数据。

图7 全采样数据与50%随机采样条件下不同算法重建结果对比

本文方法是从全局拓扑信息中预测小波系数，因此在低采样率情况下具有一定的优势。由图8可见，在20%的低采样条件下，空道占据大部分，很难重建地震数据，因此基于物理建模重建的地震数据均出现大幅度缺失或失真；数据驱动方法重建地震数据整体效果得到大幅提升，但与真实数据相比仍存在细节缺失或重建精度不够的问题；在较低采样率(≤20%)且未考虑采样数据位置的情况下，本文方法重建的地震数据较好地保留了原始地震数据特征和波形信息，在细节处更逼近实际数据。这证明了本文方法对地震数据缺失位置不敏感，在地震数据规则化的过程具有鲁棒性。

图8 全采样数据与20%均匀采样条件下不同算法重建结果对比

表1和表2分别对比了均匀采样与随机采样两种方式不同采样率条件下本文方法与其他方法测试集样本数据重建时SNR的均值对比。由表可见，本文方法显著优于基于模型的规则化方法，也优于同类基于深度学习的规则化方法。对于本文方法而言，在高采样率(>80%)时，无论是均匀采样还是随机采样，数据整体结构的特征均保持较完整，重建的数据效果比较接近；在中等采样率(20%～80%)时，均匀采样比随机采样能更好地保留数据的结构特性，均匀采样重建的效果要优于随机采样；在低采样率(≤20%)时，均匀采样和随机采样严重破坏了实际地震数据的结构特性，两种方法重建效果均较差。

表1 均匀采样条件下不同方法重建地震数据的SNR均值对比 dB

表2 随机采样条件下不同方法重建地震数据的SNR均值对比 dB

在50%随机采样条件下，比较不同深度学习方法训练效率。不同方法训练迭代1000次的训练完成时间如表3所示。由表可见，联合小波学习方法相比于时域学习方法增加了计算时间，但增加幅度相对较小，重建地震数据SNR却得到了较大程度地提高。

表3 不同深度学习方法训练所用时长及重建地震数据的SNR均值对比

3 实际地震数据测试

选取A油田实际地震数据测试本文方法规则化处理效果。震源和检波器置于地表，检波器采样间隔为2ms，道间距为12.5m。将实际样本数据共5000个按照8∶1∶1分别划分为训练集、验证集和测试集，使用训练集数据训练网络，再用测试集测试网络的有效性。任选1个测试集样本，实际地震数据不同方法重建波形对比如图9所示。由图可见，本文方法重建的地震数据同相轴光滑、连续，可较好地重建缺失道和细节，具有较好的鲁棒性。

图9 实际地震数据及不同方法重建波形对比

为了验证本文网络的泛化能力，将B采油厂叠前和叠后地震数据各抽取5000组样本，制成一个样本增广的地震数据集，训练集、验证集和测试集按不交叠的方式划分比例为8∶1∶1，通过结合两种样本数据集进行网络训练。当训练收敛时，分别使用测试集中叠前和叠后数据分别进行测试。当测试数据为叠前数据时，规则化地震数据SNR=16.0539dB；当测试数据为叠后数据时，规则化地震数据SNR=16.331dB。这表明本文网络具有一定的泛化能力和适用性。

4 结束语

本文提出的联合小波域深度学习的地震数据规则化方法有效利用了频域和时域的特征。对比基于模型的方法以及仅在时域学习条件下的重建效果，本文方法具有细节保持效果好、对地震数据缺失位置不敏感、在采样率较低的情况下具较好的重建效果的特点。实际地震数据重建结果验证了本文方法的准确性和有效性。需要指出的是，基于卷积神经网络的方法往往需要大量数据学习相应地震数据的特征，所以如何仅在有限量样本的情况下，提高模型的泛化能力，取得相对理想的效果，是进一步研究的方向。