0 引言

目前，速度建模的方法首选全波形反演(Full Waveform Inversion，FWI)^[1]。FWI方法使用全部地震波场数据(包括所有频率和位置信息)反演地球模型参数，使地震观测数据与模拟合成数据之间的残差最小以估计地下介质模型，可重建高精度和高分辨率速度模型，并具有揭示复杂地质背景下构造和岩性细节信息的潜力^[2]。然而，FWI存在计算成本高、对初始模型依赖性强、周期跳跃、易收敛于局部极小值等问题^[3]。

近年来，深度神经网络(Deep Neural Network，DNN)展示了在各种数据域之间逼近非线性映射函数的能力^[4]，在解决反问题时受到了广泛关注^[5-7]。深度学习技术在地震勘探领域也展现了良好的应用效果，例如在地震构造解释^[8-9]、地震相识别^[10-11]和地震噪声衰减^[12-13]等方面。在地震波形反演方面，DNN同样表现出很强的能力。Araya-Polo等^[14]针对速度模型建立提出了DNN解决方案，即将原始地震数据转换为与速度相关的特征立方体数据作为网络的输入，但特征提取的步骤可能会引入人为干预错误。随着全卷积神经网络(Fully Convolutional Network，FCN)的发展，数据驱动方法开始着眼于端到端网络的应用。Wang等^[15]将FCN应用于地震反演，使用改进的FCN建立速度模型。随后，Wang等^[16]进一步改进了上述方法，在全卷积层之前应用全连接层对输入数据进行预处理，直接从叠前共炮集资料中预测速度模型。Wu等^[17]设计了一个具有卷积神经网络(Convolutional Neural Network，CNN)特征的编码器—解码器地震反演网络(InversionNet)，并应用于简单断层模型的建立。Yang等^[18]采用并修改了U-Net，从原始地震记录建立速度模型。Li等^[19]提出了基于编码器—解码器架构的地震反演网络(SeisInvNet)，从每个增强的地震道提取局部和全局信息，强制每个地震道对整体速度模型的重建做出贡献。

相较于传统物理方法，上述深度学习数据驱动方法不需要提供初始速度模型，整个过程在数据驱动下进行，但性能受到标记数据集质量及数量的限制。针对该问题，Wang等^[16]使用修改后的自然图像训练网络，为网络的训练提供了大量的数据集。Araya-Polo等^[20]利用基于二维CNN的生成式对抗网络(Generative Adversarial Network，GAN)架构生成速度模型，然后将建模生成的地震数据作为训练数据输入FCN中以重建速度模型。随着上述“速度模型生成—物理建模—速度模型重建”流程的迭代，FCN提高了预测结果的准确性。上述方法都采取了生成更多标记数据集的方法。同时，一些学者尝试引入物理信息约束，缩减深度学习的参数空间。Mosser等^[21]将由低维潜在变量集有效参数化的地质构造模型与物理驱动模型相结合，并在地球物理反演中证明了该方法的有效性。Biswas等^[22]采用近似和线性化的正演模型，在叠前和叠后弹性反演任务中更新CNN，这在一定程度上避免了对标记数据集进行训练的需要。

在改进添加物理信息约束的深度学习地震反演方法时，物理信息的选取是关键。FWI中波动方程正演模拟是至关重要的步骤之一，直接影响FWI的性能。受此启发，本文引入声波方程时域有限差分正演方法，提出了波动方程正演引导的深度学习地震波形反演方法。该方法对网络预测的速度模型正演生成重建的地震数据，通过最小化速度模型误差及地震数据误差以训练网络；同时，通过可微分正演过程传递梯度，根据梯度方向动态调整地震数据误差的权重。模型对比实验和含噪数据适用性实验均证明了本文方法的优越性。

1 方法原理 1.1 物理驱动

物理驱动指通过物理信息求解问题，FWI即为利用地震波形信息获取速度模型的一种物理驱动方法。地震波在地下以类似声波的方式传播，均匀介质下二维声波波动方程表示为

$ \begin{aligned} & \frac{\partial^2 U(x, z, t)}{\partial x^2}+\frac{\partial^2 U(x, z, t)}{\partial z^2} \\ & =\frac{1}{V^2(x, z)} \frac{\partial U(x, z, t)^2}{\partial t^2}+s(x, z, t) \end{aligned} $

(1)

式中：U(x, z, t)表示在空间位置(x, z)和时间t处的声波波场；V(x, z)为(x, z)处的纵波速度；s(x, z, t)是震源函数。FWI中正演模拟通常是由给定的速度模型v模拟地震数据U的过程，即基于式(1)的正演建模问题表示为

$ \boldsymbol{U}=f(\boldsymbol{v}) $

(2)

式中f(·)代表将v映射到U的声波正演算子，是非线性的。

物理驱动的FWI可看作最小化求解问题

$ E(\mathit{v}) = \mathop {\min }\limits_{\bf{v}} \left\{ {\parallel \mathit{\boldsymbol{d}} - f(\mathit{v})\parallel _2^2 + \lambda R(\mathit{v})} \right\} $

(3)

式中：d表示真实地震数据的完整波形信息；λ是正则化参数；R(v)是正则化项，常通过v的L₂或L₁范数表示。

利用式(2)所示的正演建模过程，FWI从初始速度模型开始生成模拟地震数据，并将模拟数据与观测数据持续进行对比，通过梯度迭代更新速度模型，直至模拟地震数据足够接近观测地震数据。然而，这种方法具有一定的局限性，存在计算效率低与资源耗费大的问题。

1.2 数据驱动

不同于物理驱动的FWI，数据驱动方法将式(3)看作回归问题。根据通用近似定理^[23]，DNN可以逼近任意连续函数。因此，数据驱动方法能够直接学习地震数据到速度模型的逆映射

$ \boldsymbol{v}=\operatorname{Net}_{\boldsymbol{\theta}}(\boldsymbol{d}) \approx f^{-1}(\boldsymbol{d}) $

(4)

式中：Net_θ (·)表示由CNN实现的f(·)近似逆算子；θ表示DNN中的可训练权重。在网络训练阶段，需要大量地震数据进行学习。一旦完成训练，便可立即根据输入的地震数据得到预测的速度模型。训练需要花费大量时间，但训练完成后的预测时间几乎可忽略不计。数据驱动方法虽避免了复杂的迭代运算，但对数据集的依赖性强，在数据质量不高、代表性不强时同样会面临计算效率低下的问题。

1.3 波动方程正演引导

本文提出物理信息引导下的数据驱动方案，结合数据驱动与物理驱动，将波动方程正演建模与U-Net相连接。该方案通过数据驱动反演输出速度模型，因此不需要迭代优化每个样本；在数据驱动的基础上引入了物理信息，将网络输出推广到波动方程正演任务中。网络优化过程中地震数据损失的传播限制了参数空间，降低了网络过拟合的概率，提高了泛化性。

如图 1所示，在训练阶段，地震数据d作为标记数据输入反演网络模型，输出相应的预测速度模型$\hat{\boldsymbol{v}} $。然后，使用$\hat{\boldsymbol{v}} $正演重建地震数据$\hat{\boldsymbol{d}}$。有限差分方法^[24]是一种实现简单、模拟精度高的正演模拟方法，也是目前最常用的地震波场正演模拟方法之一。通过有限差分方法将式(1)近似为可微算子，物理正演过程可被看作“DNN层”，将梯度传递至DNN，进行正向传播和反向传播，完成端到端的训练。

如上所述，反演网络可以在正演过程的帮助下进行网络参数的学习，充分利用训练数据拟合速度模型、地震数据与对应标记数据。在预测阶段，直接将地震数据输入网络，便可获得预测速度模型。

1.4 反演网络模型架构

本文采用Ronneberger等^[25]提出的U-Net作为重建速度模型的反演网络模型。相比于普通图像数据，地震数据复杂度更高，但其分布具有较强的规律性。此外，在实际勘探中高质量的地震数据获取较为困难，因此需要对已有数据的底层特征进行充分挖掘以获取更多有用信息。U-Net不仅结合在下采样过程中提取的低分辨率特征与在上采样过程中提取的高分辨率信息，而且在上采样阶段通过跳跃连接充分结合浅层较简单的特征，能够充分挖掘、利用数据信息，对于小样本也不易过拟合，适用于地震反演问题。

U-Net输入和输出都是图像域，而地震反演建立速度模型是一个从时间域地震数据映射到深度域速度模型的过程，两者尺寸并不一致。因此，本文改进了U-Net，在网络训练前，在下采样时，将地震数据设置为与速度模型相近的尺寸输入网络；在上采样时，将最后一次卷积得到的特征图截断为与速度模型相同的尺寸。此外，将网络输入通道数修改为不同震源位置的炮集数，网络输出通道数修改为1。此时，网络从时间域地震数据提取特征并映射到深度域速度模型，实现了时深转换，建立了速度模型。

网络模型架构如图 2所示。其左半部分为收缩路径，用于获取上下文信息，即识别地震数据中某一位置与周围区域的关系，进行地质特征提取；右半部分为扩张路径，即上采样过程，能够将捕捉到的地质特征精确定位。这两条路径完全对称，构成编码器—解码器结构。收缩路径的一次下采样由两个卷积块和一个2×2池化层构成，共进行四次下采样。每个卷积块由3×3卷积层、批量归一化(BN)和ReLU激活函数组成。数据经过卷积操作后，通道数增加但尺寸不变；经过池化操作后通道数不变但尺寸减半。随着层数的加深，网络的感受野不断增大。扩展路径由一个上采样的2×2转置卷积、两个3×3卷积和一次特征拼接重复构成，每次转置卷积后，数据尺寸加倍，同时通道数减半。网络的最后一层通过一个1×1卷积将通道数变为期望的类别数。

1.5 损失函数

未引入正演过程时，反演网络仅由重建速度模型与真实速度模型之间的均方误差(Mean Squared Error，MSE)计算速度模型损失L_v，损失函数为

$ L(\boldsymbol{\theta})=L_v=\frac{1}{N} \sum\limits_{n=1}^N\left\|\boldsymbol{v}_n-\operatorname{Net}_{\boldsymbol{\theta}}\left(\boldsymbol{d}_n\right)\right\|_2^2 $

(5)

式中：N为训练数据的个数；v_n为标记速度模型；d_n为标记地震数据。引入正演过程后，由重建地震数据与标记地震数据之间的均方误差计算地震数据损失

$ L_{\boldsymbol{d}}=\frac{1}{N} \sum\limits_{n=1}^N\left\|\boldsymbol{d}_n-f\left[\operatorname{Net}_{\boldsymbol{\theta}}\left(\boldsymbol{d}_n\right)\right]\right\|_2^2 $

(6)

此时网络同时最小化速度模型的误差及地震数据的误差，损失函数扩展为

$ \begin{aligned} L(\boldsymbol{\theta})= & L_v+\lambda_d L_d=\frac{1}{N} \sum\limits_{n=1}^N\left\|\boldsymbol{v}_n-\operatorname{Net}_{\boldsymbol{\theta}}\left(\boldsymbol{d}_n\right)\right\|_2^2+ \\ & \lambda_d \frac{1}{N} \sum\limits_{n=1}^N\left\|\boldsymbol{d}_n-f\left[\operatorname{Net}_{\boldsymbol{\theta}}\left(\boldsymbol{d}_n\right)\right]\right\|_2^2 \end{aligned} $

(7)

式中λ_d为地震数据损失对反演网络影响的权重参数。一般采用Adam算法更新θ。

总损失函数包括速度损失与数据损失两项。从多任务学习的角度出发，正演任务可以视为反演任务的辅助，它帮助加快速度模型损失的下降，提高反演模型的泛化性。在训练过程中，前期反演网络输出的速度模型与真实速度模型相差较大，且边界模糊，通过正演难以得到有价值的地震数据；后期随着速度模型损失的降低，通过正演得到的地震数据精度越来越高。如果仅对地震数据损失的权重参数λ_d进行人工调整，就会难以适应这一动态的变化过程。针对这一问题，采取Lin等^[26]提出的辅助损失在线学习算法动态调整λ_d。如果正演过程和反演过程对θ产生的梯度方向一致或相近，则认为正演过程有利于反演的学习，给予较大的权重；反之，降低权重，避免影响反演学习。相较于人工调整参数，自适应辅助任务进行加权能够更好地发挥正演过程的引导作用。

2 模型测试

为了评估本文方法的性能，准备了两套不同类型的数据集，其信息如表 1所示。数据集中的地震数据作为反演网络的输入以及正演结果的标签数据，速度模型作为网络输出相比较的标签数据。

2.1 速度模型构建 2.1.1 层状模型

层状速度模型具有3~6层随机速度层，每个速度模型的网格数均为100×100，网格间距为10m，速度范围为1.5~3.5km/s。通过改变地层界面的位置、形状，可形成不同的速度模型(图 3)。

2.1.2 盐体模型

盐体模型采取Yang等^[18]生成的速度模型。其中，每个速度模型具有5~12层背景速度层及1个具有任意形状和位置的盐体，速度范围为2~4km/s，盐体速度值恒定为4.5km/s。为了更好地应用本文方法，将Yang等^[18]生成的速度模型网格数抽取为101×151，网格间距为10m(图 4)。

2.2 地震数据合成

本文采用地表观测方式，在地表均匀放置检波器。正演模拟采取有限差分法求解常密度声波方程，以Ricker子波作为震源，应用完全匹配层(Perfectly Matched Layer，PML)^[27]吸收边界条件。各数据集所采用的正演模拟参数如表 2所示。

2.3 模型对比试验

为了验证本文方法相较于数据驱动方法的优越性，选用U-Net作为基线模型。使用相同数据集分别对基线模型和本文模型进行训练，并比较训练结果。实验在4个NVIDIA Tesla P100 GPU上开展，深度学习框架为PyTorch。在训练过程中，两种模型采用相同的实验设置，具体实验参数如表 3所示。

选取峰值信噪比(Peak Signal to Noise Ratio，PSNR)和结构相似性(Structural Similarity，SSIM)两个指标，定量评估预测速度模型的质量。PSNR由MSE计算而来，是基于误差敏感的图像质量评价指标，其公式为

$ \mathrm{PSNR}=10 \cdot \lg \frac{\mathrm{MAX}^2}{\mathrm{MSE}} $

(8)

式中MAX表示生成图像的像素峰值。

PSNR越大，说明预测数值误差越小。SSIM是衡量图像相似度的指标，它分别从亮度、对比度和结构三个方面进行比较。SSIM值越接近1，说明图像相似度越高，即预测速度模型的质量越好。各个模型在测试集上的定量分析结果如表 4所示。由表可见，对于层状模型两种网络都取得了较好的结果，本文方法在两个指标上的数值略高于U-Net。对于盐体模型两种网络的预测效果有所下降，但本文方法在两个指标上均明显高于U-Net。

2.3.1 层状模型对比

图 5是层状模型测试集上随机选取的预测结果。由图可见，两种方法预测得到的速度模型都与真实速度模型高度相似，各层样貌吻合度高。本文方法对深层区域的样貌特征预测更准确(红框区域)，这说明学习正演重建地震数据有助于预测深层细节。

层状模型地质样貌较为简单且数据集数量较为充足，U-Net和本文方法都能达到较好的预测效果，但本文方法对细节的预测更准确。

2.3.2 盐体模型对比

由图 6可见，U-Net只能预测出盐体的大致形状和位置，整体吻合度较低。相较而言，本文方法预测的地层界面吻合度更高，且在盐体形状、地层样貌等细节方面预测更精确。图 7对比了图 6中样本数据在不同水平位置抽取的垂直方向速度曲线，可见本文方法预测结果更接近于真实数据，进一步证明了本文方法具有更高的预测精度。

不同于层状模型，盐体模型更复杂且数据集量较少，因此速度建模精度有所降低。但实验结果证明，在数据集有限的情况下，本文方法在复杂介质速度建模中更具优势。

2.3.3 效率对比

本文将用于地震反演问题的有限差分领域特定语言程序包Devito^[28]集成至Pytorch中以实现正演建模过程。Devito可以生成高度优化的并行代码，具有计算速度快、占用内存小的优点。

图 8为在本文实验条件下，U-Net及本文方法在两种数据集上计算1000个数据所需的时间对比。在训练过程中，本文方法需要按照表 2所示参数正演生成每一个速度模型的地震数据。由图可见，由于本文方法引入了额外的正演过程，相较于传统U-Net耗时更长。以盐体模型为例，本文方法计算1000个数据的时间比U-Net多约32s，平均每个速度模型的计算时间仅多出约0.03s，计算效率仍维持在较高水平。

2.4 含噪数据适用性实验

在实际勘探中观测得到的地震数据往往含有噪声。为验证本文方法在实际条件下的应用效果，向合成数据集中的地震数据分别加入信噪比(Signal to Noise Ratio，SNR)为20和10dB的高斯噪声以模拟真实含噪数据。其中，20dB信噪比对波形信息的影响较小，用于模拟噪声较小的情况；10dB信噪比用于模拟中等噪声情况，该情况下地震波到达深部地层时产生的较弱反射波受噪声影响较大，难以分辨。

应用本文方法分别以原始数据与含噪数据作为输入进行预测对比实验。图 9是层状模型分别在无噪数据和两种含噪数据上的反演实验结果对比。由图可见，使用20dB含高斯噪声数据作为输入得到的预测结果部分区域结构吻合度及预测数值精度有所降低(图 9c，红框区域)，但仍然能准确地反映绝大部分的地质构造以及速度值；而使用10dB含高斯噪声数据作为输入进行预测时，预测效果相对较差，出现了较为明显的假象(图 9d，红框区域)。

图 10是盐体模型分别在无噪数据和两种含噪数据上的反演实验结果对比。由图可见，20dB含高斯噪声数据作为输入得到的反演结果中除盐体形状略有偏差及局部速度结构假象外，地层位置、速度值与原始数据的预测结果大致相同；10dB含高斯噪声数据作为输入得到的反演结果相对较差，受噪声影响出现更严重的速度结构假象(红框区域)。为进一步评价本文方法，在盐体模型上应用U-Net以10dB含高斯噪声的数据作为输入进行含噪数据对比实验，结果如表 5所示。由表可知，在含噪情况下，本文方法预测精度依旧高于U-Net。同时，对比图 10d与图 10e可知，本文方法虽然出现了一定程度的速度结构假象，但整体的结构吻合度更高。因此，相较于U-Net，本文方法在抗噪能力上仍具有一定优势。

综上所述，本文方法鲁棒性较强，对小至中等强度含噪数据的适应性较强，能获得更准确的预测结果。

3 结束语

本文提出了波动方程正演引导的深度学习地震波形反演方法，结合数据驱动与物理驱动，对神经网络预测的速度模型应用波动方程正演生成地震数据。通过有限差分神经网络近似，声波方程可以传播梯度，神经网络同时学习反演及正演的不确定性参数，最小化速度模型和地震数据与对应标记数据之间的误差。为了充分发挥正演引导作用，采用辅助损失在线学习算法代替人工对地震数据损失的权重进行动态调整。相较于单一数据驱动，本文方法计算时间短，并且可充分利用有限的标记数据，在一定程度上降低了对标记数据量的依赖。在层状及盐体模型上的对比实验证明，相较于U-Net数据驱动方法，本文方法能够在相同条件下得到更准确的预测结果，并在复杂数据上具有一定的应用潜力。含噪声层状及盐体模型数据实验进一步验证了本文方法具有较强的鲁棒性，在数据含噪情况下能够得到较为可靠的预测模型。