一种融合模型和数据的全波形反演方法

A Full Waveform Inversion Method Combining Model and Data

摘要: 全波形反演算法(Full Waveform Inversion, FWI)是一种强大的地球物理成像技术，它通过迭代最小化模拟和观测地震图之间的失配来生成高分辨率地下模型。近年来，随着机器学习方法和深度学习方法的发展，研究者提出了一些基于机器学习和深度学习的全波形反演技术。其中有一种融合模型和数据全波形反演方法(physics-informed training-free frameworks for two-dimensional FWI, FWIGAN)引人瞩目，其使用生成对抗网络的思想与物理学相结合，使用无监督学习的方式来自动地估计出符合物理学规律的模型，但其网络生成图像的质量及稳定性还有待提升。本研究在此方法上引入零中心梯度惩罚(Zero-centered gradient penalty, 0-GP)，它是一种用于生成对抗网络的正则化技术，通过将判别器推向理论最优判别器来提高网络的稳定性和样本质量。结合正则化的新方法称之为(Full Waveform Inversion with Zero-Centered Gradient Penalty, FWILP) FWILP。经过实验，证明了FWILP可以提高模型质量和收敛速度，使得网络更加具有稳定性。

Abstract: Full Waveform Inversion (FWI) is a powerful geophysical imaging technique that generates high-resolution underground models by iteratively minimizing the mismatch between simulated and observed seismic waveforms. In recent years, researchers have proposed some machine learning and deep learning-based FWI techniques with the development of these methods. One such method is the physics-informed training-free frameworks for two-dimensional FWI (FWIGAN), which combines the idea of generative adversarial networks with physics to automatically estimate physically consistent models using unsupervised learning. However, the quality and stability of the generated images still need to be improved. In this study, we introduced the zero-centered gradient penalty (0-GP) to FWIGAN as a regularization technique for generative adversarial networks. This technique improves the stability and sample quality of the network by pushing the discriminator towards the theoretically optimal discriminator. The new method, combining regularization, is called Full Waveform Inversion with Zero-Centered Gradient Penalty (FWILP). Experimental results show that FWILP can improve model quality and convergence speed, making the network more stable.

1. 引言

地震速度模型对于地震勘探，尤其是偏移、成像和解释过程具有重要意义。速度模型的准确估计在地球结构和自然资源的调查中起着至关重要的作用。传统FWI [1] 用模型驱动的方法将基于正演方程数值模拟得到的模拟地震图与实测地震图进行迭代匹配，以获得高分辨率的反演地震速度模型。全波形反演的目标函数最小化问题是通过伴随状态方法 [2] 来解决的，明确计算了代价函数相对于速度模型的梯度。但这种模型驱动在面临正演算子变化时会带来挑战。

近年来深度学习的不断发展推动全波形反演数据驱动方面的研究。研究方向之一是在地震波形和速度模型配对数据训练神经网络，并取得了显著进展 [3] 。这种方法将FWI视为一个数据驱动的机器学习问题，而不是传统的基于波动方程的求解方法。但缺乏偏微分方程的约束可能会影响该方法的准确性和泛化性。另一方面的研究是利用深度学习作为信号处理工具来改进传统FWI的优化过程 [4] ，利用神经网络来推断缺失的低频信息，以减轻循环跳跃问题。

除了上述的数据驱动方法，结合神经网络和偏微分方程，将FWI视为受物理约束的机器学习问题也备受关注。理查等人 [5] 采用生成对抗网络构建了地下地质结构的先验模型，并通过优化低维潜在变量来拟合观测数据。Wu和McMechan [6] 提出了基于卷积神经网络的FWI方法，通过卷积神经网络重新参数化速度模型或梯度场，并通过更新神经网络权值来最小化损失。He和Wang [7] 进一步分析了卷积神经网络的自适应正则化效应。这些方法依赖于预训练的卷积神经网络对初始速度模型进行参数化，并将正则化效应归因于对这些初始速度模型的先验信息的拟合。而FWIGAN [8] 将生成网络用声波方程替代，并且与判别器相结合，加上梯度惩罚机制，形成了一种全新的基于物理方程和生成对抗网络的无监督学习方法。其中(Wasserstein GAN with Gradient Penalty, WGAN-GP) WGAN-GP [9] 的梯度惩罚机制，通过限制鉴别器输出的梯度范数来增加训练的稳定性，但WGAN-GP仍然存在一些问题，在某些情况下，梯度惩罚可能会导致模型的收敛速度变慢或者生成样本的多样性下降等问题。

为了解决这些问题，本研究提出了一种融合模型和数据的全波形反演方法(FWILP)，该方法在FWIGAN的基础上改进了网络，使用了零中心梯度惩罚 [10] ，其是一种有效的正则化技术，可以进一步提高网络的性能，包括更好的生成效果、更快的收敛速度和更好的泛化性能。

2. FWILP

2.1. 整体框架

FWILP运用模型驱动和数据驱动相融合的思想，将基于物理学的波动方程与深度学习中的生成对抗网络相结合。网络模型如图1所示。v代表真实速度模型，vg代表反演速度模型，地震炮集记录通过物理生成器产生，和真实地震炮集记录输入到判别网络当中。判别器网络由六个卷积块和两个全连接层组成。F是判别器的损失，real是真实地震炮集记录通过判别器得到的分数，feak是生成的炮集记录得到的分数，0-GP是零中心梯度惩罚项。首先训练判别器使其分数达到最大，此时判别器的鉴别能力达到最优；然后更新速度模型和雷克子波让判别器的分数达到最小，此时由于速度模型的更新，物理生成器生成质量更高的样本，与判别器形成最大最小的博弈。

Figure 1. FWILP network framework

图1. FWILP网络框架图

2.2. 模型驱动部分

模型驱动是指由观测数据进行正演并且进行多次对比得到实际地下数据的方法。FWILP将生成网络用正演替代，使用正演模拟的方法产生地震炮集数据。在地震波传播模拟中，把二维恒定密度声波方程(AWE)作为控制方程，其中只有纵波速度模型是作为未知参数。AWE方程表示为：

$\frac{1}{v\left(r\right)}\frac{{\partial }^{2}u\left(r,t\right)}{\partial t^2}={\nabla }^{2}u\left(r,t\right)+q\left(r,t\right)\delta \left(r-r_s\right)$ (1)

其中r表示位置，u是波场振幅，v是纵波速度分布，q是震源，t是时间。通过最小化成本函数来反演速度模型，利用伴随状态算法计算梯度，然后利用梯度信息迭代更新速度模型。伴随状态法是用于地震波成像的数值方法，该方法通过计算数据残差和梯度信息之间的关系来反演地下介质的性质，从而实现高分辨率的成像。为了通过AWE模拟波传播约束，将时域二阶精度和空间域四阶精度的规则网格和有限差分方案应用于正演当中，有限差分法的核心思想是将连续的空间或时间域离散化，通过有限差分近似来求解偏微分方程。为了防止来自模拟域边缘多余的反射，使用了完美匹配层(PML)吸收边界条件，其是一种常用的处理边界条件的方法，主要用于模拟波在无限介质中传播时的边界效应，可以有效地吸收入射波并消除反射，从而更准确地模拟波在介质中的传播过程。模型驱动部分由Deepwave [11] 模块实现，其是地震波传播数值模拟软件包，主要用于模拟地震波在地下介质中的传播过程，Deepwave有波传播模块，允许链操作和反向传播而自动重建未知变量。

2.3. 数据驱动部分

数据驱动指的是从观测数据通过学习直接映射到实际地下数据的方法。FWILP通过波动方程正演得到的地震炮集记录输入到判别器中，判别器使用的是一个由六个卷积块和两个全连接层组成的网络，其构成如图1中的判别器部分所示。数据驱动结构如图2所示。

一个块的组成部分是：具有步长1 × 1的3 × 3二维卷积层，步长为2 × 2的2 × 2最大池化层，负斜率为0.1的leaky ReLU函数。最大池化层进行下采样使得特征图的大小减半，最后一个卷积块的特征图被压扁为一个向量后输入完全连接层，最后使用leakReLU函数。最终的输出以标量分数的形式得到，这个分数可以通过测量物理生成器生成的数据和真实数据之间的差异来解释输入图像的真实程度。

Figure 2. Data-driven

图2. 数据驱动图

零中心梯度惩罚是一种用于生成对抗网络的正则化技术，旨在提高(Generative Adversarial Network, GAN) GAN [12] 网络训练的稳定性和样本质量。0-GP的基本思想是，对于鉴别器的输出，计算其关于输入数据的梯度，并要求该梯度的范数除以一个预先设定的常数值(例如1)接近于1。这个约束使得鉴别器输出的梯度不会偏向于太小或太大，把其推向理论最优鉴别器，从而增加了训练稳定性的同时提高鉴别器的质量，且适用于模型和数据融合驱动的方法。

0-GAN网络的基础上增加了以零为中心的梯度惩罚项：

$L_{\text{0-GP}}={{\rm E}}_{x~p_r}\left[\log \left(D\left(x\right)\right)\right]+{{\rm E}}_{z~p_z}\left[\log \left(1-D\left(G\left(z\right)\right)\right)\right]-\lambda {{\rm E}}_{\tilde{x}}\left[{\Vert {\left(\nabla D\right)}_{\tilde{x}}\Vert }^2\right]$ (2)

$p_g$ 是生成器在数据x上的分布， $p_z\left(z\right)$ 是输入噪声变量， $G\left(z;{\theta }_g\right)$ 是到数据空间的映射，其中G是一个可微函数，由具有参数 ${\theta }_g$ 的多层感知机表示。 $D\left(x;{\theta }_d\right)$ 是多层感知器，它输出单个标量。 $D\left(x\right)$ 表示x来自于数据而不是 $p_g$ 的概率。训练D以最大化将正确标签分配给训练示例和来自G的样本的概率，同时训练G以最小化 $\log \left(1-D\left(G\left(z\right)\right)\right)$ 。最后一项是零中心梯度惩罚项，其中 $\tilde{x}=\alpha x+\left(1-\alpha \right)y$ ， $x~p_r,y~p_g$ ， $\lambda$ 是一个超参数，控制惩罚项的大小。

通过使用0-GP，FWILP能够有效地解决GAN训练过程中的不稳定性问题，同时还可以提高模型的样本质量。

3. 实验和结果

3.1. Marmousi速度模型介绍

Marmousi [13] 是勘探地球物理学中最著名和最标准的声速模型之一。其基于地质学和地震学数据构建地下模型，包含了多个地层，包括沉积岩、盐岩、泥岩和砂岩等，被广泛应用于地震勘探和地震成像领域，以帮助研究人员理解地下结构和资源分布。FWILP仅使用纵波速度进行数据模拟和反演，并且对其进行了下采样，大小为100*310，初始模型为线性初始模型。由于FWI方法对于初始模型较为敏感，相较于可以反映宏观背景速度的高斯模型，使用线性初始模型对其反演则更具有挑战性，实验中使用的线性初始模型可以更加反映方法的效果。线性初始模型和重采样真实速度模型如图3所示。

Figure 3. Initial model and resampling model

图3. 初始模型和重采样模型图

3.2. 评价指标

为了评估提出的方法的质量，与FWIGAN使用了相同评估方法，定量评估反演速度模型v_g相对于真实速度模型v的质量。我们计算结构相似性(SSIM)，定义为

$\text{SSIM}\left(v_g,v\right)=\frac{\left(2{\mu }_v{\mu }_{v_g}+c_1\right)\left(2{\sigma }_v{\sigma }_{v_g}+c_2\right)}{\left({\mu }_v^2+{\mu }_{v_g}^2+c_1\right)\left({\sigma }_v^2+{\sigma }_{v_g}^2+c_2\right)}$ (3)

其中 ${\mu }_v$ 、 ${\mu }_{v_g}$ 、 ${\sigma }_v$ 、 ${\sigma }_{v_g}$ 分别是图像v和v_g的局部均值、标准差和交叉协方差。正则化常数 $c_1={10}^{-4}$ ， $c_2=9\times {10}^{-4}$ 。

速度模型v_g与真实速度模型v之间的相对误差，定义为：

$\text{ERROR}\left(v_g,v\right)=\frac{{\Vert v-v_g\Vert }_2}{{\Vert v\Vert }_2}$ (4)

其中 ${\Vert \cdot \Vert }^2$ 表示 ${\mathcal{l}}_2$ 范数。

SNR越大，说明噪声在信号当中的比例越小，图像质量越高，否则相反。SSIM是指两幅图片相似度的指标，其数值越大表示图像和预期图像的差距越小，即图像质量越好。相对误差ERROR是指测量所造成的绝对误差与被测量真值之比，数值越小误差越小。

3.3. 实验结果与结果分析

FWILP、FWIGAN和真实地下的地震炮集数据如图4所示：

由图4可以看出，FWILP的地震炮集数据的构建更加平滑，连接更为顺畅，顶部线条更为平滑，构建相对更加完整。在初始模型更加有难度的情况下，FWILP的地震炮集数据更贴近于真实炮集数据。

FWIGAN与FWILP反演结果如图5所示：

由图5可以看出，FWILP的右部速度较低的岩性恢复更好，纹理清晰，走势明朗；整体下部速度更高的剥蚀地貌脉络清晰，构造完整；右下部速度构造较为复杂的背斜部分受影响更小，岩层分布清晰可见，形态更为相近。由此可见，在初始模型更加具有挑战性的时候，FWILP更忠于真实的地震速度模型。

Figure 4. Seismic gun set data comparison diagram

图4. 地震炮集数据对比图

Figure 5. Comparison of inversion models

图5. 反演模型对比图

本实验三个指标结果如图6所示：

Figure 6. Comparison of metrics

图6. 指标对比图

由图6可以看出，FWILP在三个测量指标上均有一定的进步，SNR更大说明所反演出的速度模型噪声更小，质量更高；SSIM更大说明输出和真实的速度模型之间的差距更小，两个速度模型更为相似；相对误差ERROR更小说明表示输出和真实速度模型的差距更小，错误率更低，反演结果更加真实。

4. 结论

本研究在FWIGAN的基础上提出了一种融合模型和数据的全波形反演方法，称之为FWILP。FWIGAN在地震速度反演任务中显示出了令人鼓舞的结果，但它在捕获精细尺度细节和保留局部特征方面存在局限性。为了解决这些问题，引入了零中心梯度惩罚正则化技术，改善了网络的性能。实验结果表明，FWILP在反演质量和稳定性方面都具有良好的性能。定量评估显示，与FWIGAN相比，误差有所降低，地下结构还原度提高，能够以更高的保真度恢复复杂的地质特征。

尽管FWILP比FWIGAN有进步，但必须承认其局限性。在更加有挑战性的常数初始模型中，并不能很好的构建出速度模型。未来的研究可以探索在不影响反演质量的情况下优化计算效率的策略，降低其对初始模型的敏感性。

综上所述，FWILP算法在地震图像重建方面有了一定程度的改进。结合正则化技术，提高了反演精度，忠实地反映了精细尺度的地质特征。本研究突出了融合模型和数据在各种地球物理成像应用中的潜力，并为地震反演领域的进一步发展铺平了道路。

基金项目

河北省重点研发计划项目(项目名称：基于时空大数据及深度学习的地质灾害风险识别关键技术研究编号：22375415D)；

河北地质大学创新团队项目KJCXTD-2021-11。

NOTES

^*通讯作者。