论文阅读记录之《VelocityGPT 》

一、基本介绍

  标题：《Propagating the Prior From Shallow to Deep With a Pre-Trained Velocity‐Model Generative Transformer Network》——通过预先训练的速度模型生成变压器网络将先验信息从浅层传播到深层

  作者：Randy Harsuko, Shijun Cheng, and Tariq Alkhalifah

  期刊：JGR: Machine Learning and Computation（CFF-A类，3区）

  发表时间： Received 26 AUG 2024；Accepted 29 JAN 2025

  摘要：构建地下速度模型，对于利用地震数据进行地球发现、勘探以及监测至关重要。随着机器学习的兴起，这些速度模型（更准确地说是它们所服从的分布）可以被精准、高效地保存在一个生成式模型之中。一旦完成存储，这些速度模型的分布便可在未来的诸多应用中发挥作用，例如在求解反问题时作为正则化项或量化不确定性的工具。然而，现有的多数生成器（如归一化流或扩散模型）将图像（即速度模型）视为整体，并未显式地利用空间依赖关系及分辨率随深度的变化。为弥补这一缺陷，我们提出 VelocityGPT——一种全新的实现方案：它采用自回归训练的 Transformer 解码器，自上而下（由浅至深）地逐层生成速度模型。由于地震数据通常在地表采集，这种“自上而下”的生成器可把从浅层反演得到的信息作为先验，用以引导深层建模。具体而言，浅层速度分布会进一步约束并指导深部速度的预测。为便于实现，我们引入了一个附加网络对速度模型进行压缩。此外，我们还将测井数据或构造信息（以偏移成像结果表示）等先验注入生成过程。在合成数据上的实验表明，VelocityGPT 在地震速度模型构建的生成应用中展现出巨大潜力。

二、介绍

  理解底层数据分布在统计学、机器学习及数据科学等诸多领域都至关重要。其优势包括但不限于：量化不确定性、合理权衡决策权重以及生成新的观测样本（Tomczak, 2021）。传统上，学习数据分布主要依赖统计方法，可分为参数化与非参数化两类。参数化方法假设数据服从某种特定分布，该分布由有限个参数刻画，例如由均值和方差定义的正态分布；典型代表包括线性回归和高斯混合模型（Reynolds, 2009）。相比之下，非参数化方法不对数据分布预设固定形式，因而能适应更广泛的数据分布形态；典型代表包括核密度估计（Parzen, 1962）与近邻方法。近年来，深度学习方法，尤其是生成式模型，已被证明是刻画复杂数据分布的高效工具。

  生成式模型因其能够捕捉复杂的数据分布，在众多科学领域备受关注。在地震数据分析中，它们可从已学习到的分布中进行采样，从而被广泛应用于多种任务，包括扩充观测数据。典型应用包括：地震反演（如 Meng 等，2021）、作为反演约束（如 Wang 等，2023），乃至不确定性量化（如 Siahkoohi 等，2020）。此外，Wang 等（2024）指出，在训练阶段可将我们对地下的先验知识（例如测井数据和地震图像）纳入网络，使其在推断阶段能够进行条件采样。Taufik 等（2024）进一步展示了生成式模型在地震反问题中的高级应用：扩散模型被证明可有效存储多参数地下模型的分布，并用于正则化弹性全波形反演（FWI）算法。然而，上述方法大多基于 CNN 架构，而 CNN 是为图像任务量身定制的。

  在地震反演问题中，核心任务之一是重建地下速度模型。理论上，最佳做法是在目标区域周围同时布设震源与接收器；然而这种布置在经济与技术层面几乎不可行。因此，工业界普遍将地震数据采集局限在地表。结果，我们对深部（资源主要富集区）的分辨率，完全取决于对浅部地层的刻画精度。当近地表存在显著横向非均质性时——这在陆上地震勘探中屡见不鲜（Bakulin & Silvestrov, 2021）——这一需求变得尤为迫切。在速度建模领域，最著名且广泛应用的技术之一便是“层剥离”（layer stripping）（Yilmaz, 2001）：自浅至深，一次反演一层或一个垂直网格点，迭代推进。近期，一些地质建模方法也利用观测到的近地表资料（如测井数据）来约束地下重建（Guo et al., 2021；Jia et al., 2021；Wu et al., 2005），进一步彰显了数据驱动的“自上而下”策略的优势。回到生成式模型，当前主流框架——如基于流的模型（如归一化流 Rezende & Mohamed, 2015）与潜变量模型（如生成对抗网络 Goodfellow et al., 2014；扩散模型 Ho et al., 2020）——并未显式具备自上而下逐层生成速度模型的能力。相比之下，自回归模型天然支持序列预测，能够以先前生成的部分作为先验。过去，我们曾使用 CNN（如 PixelCNN，Van Den Oord et al., 2016）和 RNN（如 PixelRNN，Van Den Oord et al., 2016）作为自回归骨干网络，但前者感受野受限，后者计算效率低下。ImageGPT（Chen et al., 2020）首次将 Transformer 解码器（Vaswani et al., 2017）引入自回归图像合成，凭借 Transformer 的可扩展性与长程依赖建模能力，为这一问题提供了全新的解决思路。

  本文提出了一套面向地震速度建模的全新生成式框架——VelocityGPT。该框架倡导“自上而下”逐层生成速度模型的理念，并以自回归方式训练 Transformer 解码器来实现。VelocityGPT 通过训练捕获速度模型的内在分布，从而为多种潜在应用奠定基础，例如扩充训练样本、正则化速度反演以及量化不确定性。不过，本文仅聚焦于评估框架对速度分布的学习与表征能力，核心依据是生成样本的质量。具体贡献总结如下：

1. 提出一种新颖框架，可自然地以“自上而下”顺序合成地震速度模型。
2. 展示多种将先验知识（地质特征、井速度、叠后地震图像）融入生成过程的方法。
3. 在 OpenFWI 基准数据集上验证所提框架，证实 VelocityGPT 具备稳健的速度模型生成潜力。
4. 证明在小尺度速度模型上训练的 VelocityGPT，无需再训练即可直接生成更符合真实尺度的速度模型。

  本文结构如下：第 3 节首先概述整体框架，包括对速度模型进行编码并随后存储其分布。接着，第 4 节详细阐述每一步的实现细节及训练策略，并展示如何引入地质类别、测井信息和构造特征（来自图像）等多种条件，使得在推断阶段能够进行条件采样。最后，第 5 节和第 6 节分别对结果进行讨论并给出结论。

三、理论

3.1 框架

  我们提出的多条件自回归速度生成器框架如图1a所示。在该框架的第一阶段，我们训练了一个向量量化变分自动编码器（Vector-Quantized Variational Auto Encoder，VQ-VAE，Van Den Oord等人，2017），该编码器可将速度模型的重叠区块转换为离散形式（类似类别）。随后，通过训练生成式预训练Transformer（Generative Pre-trained Transformer，GPT）模型（Radford，2018）来存储速度模型的分布。该GPT模型由Transformer解码器堆叠而成，包含一个嵌入层和一个输出层，训练过程采用自回归目标函数。速度模型在作为GPT模型输入时的表示方式，自然决定了模型会采用自上而下的预测方式。如图2所示，速度模型（或推理阶段中的不完整速度模型）经VQ-VAE压缩（编码）后，会进入离散潜在空间，其中每个层都由一组整数表示。之后，这些整数按光栅化顺序（从左上到右下）展平为一个长整数序列，该序列即作为GPT模型的一个输入样本。GPT模型会预测对应速度模型更深层的下一组（或多组）整数。最后，通过对整数序列进行重塑，并利用VQ-VAE对离散潜在空间进行解码这两个反向操作，得到完整的速度模型。各阶段的详细内容将在后续章节中展开讨论。

图 1.（a）VelocityGPT 的整体框架与网络结构。该框架分为两个阶段：首先训练向量量化变分自编码器（VQ-VAE），将速度模型转换为离散表示；其次训练自回归的生成式预训练 Transformer（GPT），以学习速度模型的分布。（b）单个 Transformer 解码器模块的示意图。

图 2. 速度模型编码/解码过程示意图。离散潜空间（中间）中的纯色区域对应速度域（左侧）中相应阴影区域；虚线框表示即将生成的区域。

3.2 向量量化变分自编码器（VQ-VAE）

  框架的第二阶段涉及GPT模型训练，该阶段通过自回归分类来生成新的速度样本。由于此方法具有分类性质，我们首先需要将速度模型从连续空间转换到离散空间。获取输入的离散表示有多种方法，其中Van Den Oord等人（2017）提出了一项开创性工作，他们使用一种中间包含嵌入向量的改进型自编码器，将多种模态（音频、图像、视频）转换为其离散表示。近期研究在此概念基础上，对训练过程或网络组件进行了改进（例如Esser等人，2021；Yu等人，2021）。为简单起见，本研究采用向量量化变分自编码器（VQ-VAE），该模型包含一个卷积自编码器，并涉及嵌入实现（Van Den Oord等人，2017）。

  向量量化变分自编码器（VQ-VAE）的结构如图1a所示，该模型为双层编解码器结构，以嵌入矩阵作为瓶颈层。每一层均采用核大小为4、步长为2的二维卷积（2D convolution）。在该配置下，VQ-VAE的压缩比为 $2^l$，其中l代表层数。因此，本文所采用的VQ-VAE可将速度模型从64×64的连续值压缩为16×16的离散（整数）值。为将速度模型的连续空间量化为离散空间，需选定码本大小（即唯一离散值的数量）$D$。给定数据集 $D={(x_i,e_i)}_{i=1}^N$（i从1到N），其中$x_{ᵢ}$表示输入速度模型，$e_{ᵢ}$表示模型训练过程中得到的离散潜在表示（初始值随机初始化），则采用以下目标函数对网络进行优化（Van Den Oord等人，2017）：

$${\mathcal{L}}_{\text{VQ-VAE}}=\log p(x\mid z_d(z_e(x)\left)\right)+\Vert \text{sg}[z_e(x)]-e{\Vert }_2^2+\beta \Vert z_e(x)-\text{sg}[e]{\Vert }_2^2,\text{\hspace{1em}(1)}$$

  其中，$z_e$ 是编码器，$z_d$ 是解码器，$s_g$ 是停止梯度算子，其定义为在前向传播中为恒等映射，在反向传播中为零映射。目标函数中各项的作用如下：第一项是重建损失，用于优化编码器和解码器；第二项是嵌入损失，用于更新字典；最后，第三项是承诺损失，用于约束编码器和嵌入层。在公式1中，我们采用重建损失的一般形式，即给定网络输出时输入$x$的对数概率分布，因为可以使用任何相似性函数。在我们的方法中，我们测试了常用的L2和L1范数作为重建损失，发现L2范数在整体像素重建方面效果更好。总目标函数表示为所有训练样本的损失之和：

$${\mathcal{L}}_{\text{total}}=\sum _{i=1}^N{\mathcal{L}}_{\text{VQ-VAE}}(x_i,e_i)，\text{\hspace{1em}(2)}$$

3.3 生成式预训练Transformer（Generative Pre-trained Transformer）

  在自回归训练中，我们的思路借鉴了ImageGPT（Chen等人，2020）的方法。在ImageGPT中，每个像素被视为一个标记（token），其训练目标是根据前面的像素（或其量化版本）预测下一个像素（或其量化版本）。ImageGPT基于Transformer解码器构建（Vaswani等人，2017）。我们采用了ImageGPT模型，并对其进行改进，使VelocityGPT能更好地适配输入——即速度模型的离散潜在编码（而非像素）及其对应条件。该模型包含三个主要组件：嵌入块（embedding block）、解码器块（decoder block，可多个）和预测头（prediction head）。输入为经过训练的VQ-VAE输出的编码速度模型的独热向量（one-hot vector）e∈{0,1}T×De \in \{0,1\}^{T \times D}e∈{0,1}T×D（满足${\sum }_{d=1}^D{e}_d=1$）。其中，$S$表示速度模型经展平后的独热编码离散潜在表示的长度，$D$表示速度码本大小。在嵌入块中，我们使用嵌入层$E_1$，将每个独热编码的离散潜在编码投影到隐藏维度$H$，具体过程如下：

$$E_1(e)=e{W}_{E_1}+b_{E_1},\text{\hspace{1em}(3)}$$

其中，$W_{E1}\in {\mathbb{R}}^{D\times H}$ 和 $b_{E1}\in {\mathbb{R}}^H$ 分别是嵌入层的权重和偏置。然后，将一个可学习的位置嵌入 $E_2$ 加到经线性投影的数据上：

$$E_2(r)=r{W}_{E_2}+b_{E_2},\text{\hspace{1em}(4)}$$

其中，$r\in {\mathbb{Z}}^S$ 是离散潜在编码的相对位置，$W_{E2}\in {\mathbb{R}}^{S\times H}$ 和 $b_{E2}\in {\mathbb{R}}^H$ 分别是位置嵌入层的权重和偏置。因此，嵌入块的输出为：

$$E=E_i(e)+E_2(r).\text{\hspace{1em}(5)}$$

  接下来，嵌入后的速度模型$\hat{E}$被输入到由$L$个解码器块组成的堆叠结构中。单个解码器块的示意图如图1b所示，该解码器块由一个掩码多头自注意力（Masked Multi-Headed Self-Attention，MMSA）层、一个逐位置前馈（Position-Wise Feed-Forward，FFN）层，以及层归一化（Layer Normalization）和残差连接（Residual Connections）共同构成。在掩码多头自注意力（MMSA）层内部，首先将嵌入后的速度模型$\hat{E}$通过线性投影转换为三个不同的矩阵，即查询矩阵（Query，$Q$）、键矩阵（Key，$K$）和值矩阵（Value，$V$）：

$$Q=\hat{E}{W}_Q+b_Q,K=\hat{E}{W}_K+b_K,V=\hat{E}{W}_V+b_V,\text{\hspace{1em}(6)}$$

其中，$W\in {\mathbb{R}}^{H\times H}$ 和 $b\in {\mathbb{R}}^H$ 分别是查询、键和值矩阵线性投影的可学习权重和偏置（通过各自的下标来区分）。随后，对 $Q$、$K$ 和 $V$ 执行掩码注意力（MA）操作，其定义为：

$$X=MA(Q,K,V)=\text{softmax}\left(\frac{Q{K}^{\mathrm{T}}+M}{\sqrt{H}}\right)V,\text{\hspace{1em}(7)}$$

其中，$M\in {\mathbb{R}}^{S\times S}$ 是一个掩码矩阵，其对角线及下三角区域的元素为0，其余区域（上三角区域）的元素为负无穷（$-\infty$）。该掩码可强制注意力图（attention map）成为下三角矩阵，从而避免模型在输入序列中利用未来信息（即后续的速度层）。速度模型离散潜在编码之间的关联（权重）存储在查询-键点积（query-key dot product）中。这些操作会在注意力头（attention heads）数量 $A$ 上分布式执行，具体方式为：先对输入进行拆分，执行操作后再将结果拼接回去。

  随后，矩阵 $X$ 会与残差连接（residual connection）相加，并通过层归一化（layer normalization）进行归一化处理。接着，该矩阵被输入到前馈网络（FFN）中。此FFN是一个简单的双层全连接神经网络，中间采用高斯误差线性单元（Gaussian Error Linear Unit，GELU）作为激活函数（Devlin等人，2018；Vaswani等人，2017）。FFN层的输出表示为：

$$F={\psi }_1(X{W}_1+b_1)W_2+b_2.\text{\hspace{1em}(8)}$$

其中，${\psi }_1$ 为高斯误差线性单元（GELU）函数，$W$ 为权重，$b$ 为偏置。同样地，需对 $F$ 施加残差连接。最后，最终解码器块的输出会传入预测头 $P$ 中。该预测头的作用是将隐藏维度 $H$ 投影到速度嵌入维度 $D$，因此采用了一个线性层和双曲正切（Tanh）激活函数 ${\psi }_2$，具体形式如下：

$$P(F)={\psi }_2(F{W}_P+b_P),\text{\hspace{1em}(9)}$$

  原始ImageGPT通过逐像素预测来生成图像。然而，我们并非预测像素，而是基于离散向量的前序序列，预测一个代表速度模型区块（经VQ-VAE编码）的离散向量。在此过程中，输入概率分布可表述为：

$$p(e)=\prod _{t=1}^{T}p(e_{{\pi }_t}\mid c,e_{{\pi }_1},\dots ,e_{{\pi }_{t-1}},\theta ),\text{\hspace{1em}(10)}$$

其中，${\pi }_i$ 是离散向量 $e$（从VQ-VAE编码器输出获得）按光栅化顺序（从左上到右下）展平后的位置，$T$ 是序列长度，$c$ 是应用于生成器的速度模型条件。因此，我们采用交叉熵损失（CE）来优化模型：

$${\mathcal{L}}_{\text{CE}}(p,q)=-\sum _{t=1}^T\sum _{d=1}^{D}p(e_{t,d})\log q(e_{t,d})\text{\hspace{1em}(11)}$$

其中，$D$ 是码本大小，$p(e_{t,d})$ 是第 $t$ 个标记的真实概率分布，$q(e_{t,d})$ 是预测分布（即VelocityGPT模型的输出）。公式9为条件的直接实现提供了可能，具体可通过在速度模型的离散向量序列前附加相应条件来实现。正如Zhan等人（2022）所展示的，多种图像条件（如草图、文本、音频）被附加在图像的离散潜在表示输入之前，形成一个长向量。在本研究中，为了引入类别信息，我们构建了一个新的嵌入层，将速度模型的10个类别转换到潜在空间，并将其附加在嵌入后的离散向量之前。为了注入井信息，在训练过程中，我们从每个速度模型中随机选取一个垂直剖面，并使用其由16个离散元素组成的离散表示。然后，通过另一个嵌入层，将这16个离散元素连同井位信息一起投影到潜在空间，并将其附加在输入序列之前。

  类似的方法可用于整合地震图像中的构造信息。为了获取每个速度模型的图像，我们采用一维叠后线性建模算子：

$$I(z)=\frac{1}{2}s(z)⊛\frac{\mathrm{d}\ln Z(z)}{\mathrm{d}z},\text{\hspace{1em}(12)}$$

其中，$I$、$s$ 和 $Z$ 分别代表叠后数据（post-stack data）、子波（wavelet）和作为深度函数的声波阻抗（acoustic impedance）。声波阻抗的正式定义为：

$$Z(z)=V(z)p(z),\text{\hspace{1em}(13)}$$

其中，$V$ 和 $\rho$ 分别代表深度方向上的速度模型和密度模型。我们通过加德纳方程（Gardner’s equation，Gardner等人，1974）由速度估算密度。随后，为压缩图像中丰富的构造信息，我们采用了另一个向量量化变分自编码器（VQ-VAE），其架构与用于处理速度模型的VQ-VAE相似。

  然而，由于该图像与速度模型具有相同的维度，若将图像直接附加在序列前，意味着要向输入序列中额外引入256个（16×16）元素。对于我们的Transformer解码器模型而言，过长的序列是不利的——因为计算成本会随输入序列长度的增加呈二次方增长。因此，如图1a所示，我们转而将图像的潜在表示（经图像嵌入层投影得到）与速度模型的潜在表示一起，添加到实际空间中与其对应的位置上。

四、数值实验

4.1 数据集

  我们采用公开可获取的OpenFWI数据集（Deng等人，2022），该数据集包含一系列多样的速度模型及其对应的炮集记录（shot gathers）。根据研究需求，我们仅提取其中的速度模型——这些模型被划分为8个类别，分别代表不同的地质特征与复杂程度，具体包括：FlatVelA（FVA，平层速度模型A）、FlatVelB（FVB，平层速度模型B）、CurveVelA（CVA，弯曲层速度模型A）、CurveVelB（CVB，弯曲层速度模型B）、FlatFaultA（FFA，平层含断层速度模型A）、FlatFaultB（FFB，平层含断层速度模型B）、CurveFaultA（CFA，弯曲层含断层速度模型A）、CurveFaultB（CFB，弯曲层含断层速度模型B）（见图3.1）。所有速度模型均从原始的70×70尺寸裁剪为64×64，每个像素的速度值范围为1500至4500米/秒（m/s）。此外，我们新增了两个速度模型类别：CurveVelS（CVS，平滑弯曲层速度模型）与FlatVelHR（FVHR，高分辨率平层速度模型）。其中，CVS通过对CurveVel（弯曲层速度模型）类别进行平滑处理得到，代表平滑型速度模型（例如可作为反演任务的初始模型）；FVHR则通过从FlatVel（平层速度模型）类别中随机选取4个模型，经组合与压缩后形成单个样本，代表高分辨率速度模型（例如可用于储层建模任务）。最终，我们共获得38.1万个速度样本，将其分别划分为32.4万个训练样本与5.7万个验证样本。同时，我们将所有速度值归一化到[-1, 1]区间内。

  在图像方面，我们采用典型的叠后图像设定，选择峰值频率为20赫兹（Hz）、采样率为16毫秒（ms）的雷克子波（Ricker wavelet）作为震源子波。随后，通过公式11与公式12，可生成每个速度模型对应的合成地震图像。图像示例如图5第一列所示。

图 3. OpenFWI 数据集中各类型速度模型示例。CurveVelS 与 FlatVelHR 为自定义类别：前者通过对 CurveVelA 样本做平滑得到，后者则从 FlatVelA 与 FlatVelB 样本中随机拼接而成。

4.2 向量量化变分自编码器（VQ-VAE）训练

  VQ-VAE的架构如第3节所述，由4个卷积层（2个编码器层、2个解码器层）组成，中间设有一个嵌入层。本研究共构建两个VQ-VAE模型：一个用于对速度模型进行编码（记为VQ-VAEₐ），另一个用于对叠后图像进行编码（记为VQ-VAEᵢ）。对于速度模型编码的VQ-VAEₐ，其嵌入维度设为128，输入与输出卷积层的通道数为128，中间卷积层的通道数为68，模型总参数量为440,641。而对于叠后图像编码的VQ-VAEᵢ，其嵌入维度同样设为128，所有卷积模块的通道数均为32，模型总参数量为79,841。

  本研究选择L2损失作为重建损失，且参数β设为0.1。在其他训练参数设置方面：

• VQ-VAEₐ采用的批大小（batch size）为1024，优化器选用Adam，学习率设为4e-3；
• VQ-VAEᵢ采用的批大小为128，优化器同样选用Adam，学习率设为1e-3。

  训练过程中引入早停（early stopping）机制，其容忍 epoch 数（patience）设为10——当验证损失在10个epoch内未出现下降时，早停机制将触发并终止训练。

  在上述设置下，使用单块NVIDIA RTX 8000 GPU进行训练：

• VQ-VAEₐ共训练142个epoch，耗时8.3小时；
• VQ-VAEᵢ共训练82个epoch，耗时1小时。

  图4展示了VQ-VAEₐ对两类速度模型样本的重建结果。该模型能够以极小误差重建速度模型（见图4最右侧列），仅在速度突变的界面处可观察到少量损失。VQ-VAEᵢ的图像重建结果可通过图5的两个样本观察到：具体而言，该模型在还原极端构造时表现稍显不足，但整体图像的重建效果良好。

图 4. 向量量化变分自编码器（VQ-VAE）训练完成后，对速度模型进行重建的两个示例。左侧对应 FlatVelA 类别 (a)，右侧对应 CurveFaultA 类别 (b)。

图 5. VQ-VAE 训练完成后，对叠后地震图像进行重建的两个示例：(a) 对应 FlatVelA 类别，(b) 对应 CurveFaultA 类别。

4.3 GPT训练

  描述Transformer架构配置的标准符号包括：$L$（解码器层数）、$H$（隐藏/潜在维度大小）和$A$（注意力头数量）。在本研究的测试中，GPT模型配置为$L=8$、$H=128$、$A=4$，总可训练参数为1,680,768个。模型采用Transformer解码器的层前归一化（pre-layer-normalization）变体（Xiong等人，2020），其中层归一化模块置于自注意力模块和前馈模块之前，如图1b所示。

  网络训练采用交叉熵损失，计算输出与输入速度模型离散向量的概率分布之间的差异（公式10），这本质上是一种自监督训练。我们采用了教师强制（teacher-forcing）训练策略（Toomarian和Barhen，1992），将标签作为网络输入，这是现代自回归模型训练中的常用方法。网络参数优化使用Adam优化器，学习率为4e-3，批大小为1,024。结合早停机制，训练共进行102个epoch，在单块NVIDIA RTX 8000 GPU上耗时约10.5小时。

4.4 推论

  VelocityGPT的采样机制通过从条件分布$q(e_t|e_{<t},c)$中迭代抽取标记集来实现自回归序列生成，其中$c$代表上下文（即条件和历史标记）。在每个时间步$t$，模型生成对数几率$y_t$，并将其转换为概率分布$q(e_t)=\text{softmax}(y_t)$。随后通过多项式采样从$q(e_t)$中抽取标记$e_t$，并将该标记追加到序列中，用于约束后续预测。这一过程不断迭代，直至生成长度为$T$的完整序列，最终得到分布$q(e)$的样本。

  本节将测试VelocityGPT生成新速度模型样本的能力。我们首先进行无条件采样测试，让网络自由生成速度模型；随后逐步添加不同条件（类别、井信息和图像）以约束生成过程，观察网络对这些条件的处理效果。

  我们首先测试训练后网络的无条件采样能力：仅向网络输入浅地表信息作为起点，让网络基于该输入生成不同的实现结果。如图6所示，网络能够生成多样化的速度模型样本，且这些样本与输入的浅层部分保持一致。

图 6. 以验证集中随机抽取模型的浅部（右列所示）为输入，在各类别条件下生成的样本。

  接下来测试类别条件采样：固定输入的浅层速度信息，改变输入序列前附加的类别标签。在所有条件采样示例中，我们以CurveFaultB类别的一个样本（见图7最右侧列）作为参考，称之为"代表性样本"。图7a展示了类别条件采样的结果，可见网络生成的速度模型既符合浅层结构，也与指定类别一致。需注意的是，生成的速度模型深部不受约束（图7b）。

图 7. (a) 以第一列为输入，按速度类别进行条件采样；各类别标注于对应子图顶部。 (b) 对应 (a) 中自左至右各列在 X=32 位置的一维速度剖面；蓝色线为生成速度，橙色线为井速度。

  随后，我们在输入序列前附加代表性样本的井剖面（取自X=32位置），以测试网络的井条件采样能力。从图8a和8b可以观察到，网络现在生成的深部速度模型与井速度信息一致，尤其在井位处表现明显。这表明网络已学会适应给定的先验信息（浅层部分、类别和井数据），而模型的其他结构在不同样本中仍保持多样性。

图 8. (a) 在第一列为输入的条件下，同时以类别和井资料为约束进行采样。 (b) 对应 (a) 中自左至右各列在 X=32 位置的一维速度剖面；图例含义与图 7 相同。

  最后，我们添加代表性样本的图像以向速度模型注入构造信息。图9c（最右侧列）显示了作为网络输入的图像样本。如图9a-9c所示，生成的样本不仅符合前文提到的所有条件（浅层部分、类别和井数据），其构造也受到有效约束（图9c）。对于某些类别（尤其是含平层的类别），由于类别信息与施加的构造信息存在冲突，预测结果可能出现失真，导致网络无法生成地质上合理的速度模型（详见第4节）。然而，当浅层、类别、井数据和图像信息完全一致时，网络生成的速度模型及其对应图像与真实值非常接近（对比图9a第二行第四列与第九列的速度模型，以及图9c第二行第三列与第五列的对应图像）。这表明该模型在生成速度模型样本时能够处理多条件约束。

图 9 (a) 以第一列为输入，在“类别 + 井资料 + 图像（c，最右列）”三重条件下进行采样。 (b) 对应 (a) 中自左至右各列在 X = 32 处的一维速度剖面；图例同图 7。 (c) 由生成速度计算得到的图像，对应 (a) 中自左至右各列；图例同图 7。

  利用VelocityGPT的自回归特性，可调整采样方式。在推理阶段的序列生成过程中，可任意更改前置的类别标签，迫使网络生成符合对应类别约束的下一个离散向量。因此，通过将序列位置与其在速度模型中的原始位置关联，我们能够对生成样本进行空间控制并融合不同类别，包括在特定位置（如储层区域）注入高分辨率类别，如图10所示。

图 10. 以第一列所示类别为空间条件生成的样本。

4.5 基于SEAM Arid模型的分布外案例研究

  在前一节中，我们阐述了如何利用训练好的VelocityGPT生成速度模型。然而，在实际场景中，地震速度模型的维度通常远大于我们训练时所用的模型。因此，下文将介绍一种方法，无需额外训练网络即可将预训练的VelocityGPT应用于更接近实际尺寸的模型，并将展示该网络的条件采样功能在生成过程中仍能有效发挥作用。

  为了验证这些特性，我们以SEAM Arid速度模型（Oristaglio，2015）作为参考，同时将其作为分布外测试案例。SEAM Arid模型是一个复杂的地球浅地表模型，模拟了沙漠环境。如图12a所示，该模型尺寸为832×300，远大于用于训练VelocityGPT的OpenFWI数据集（64×64）。我们再次考虑这样一种场景：模型的浅层部分已知（例如，通过近地表反演方法或先验知识获得），目标是基于近地表信息生成模型的深部部分。因此，为使预训练网络适应这些设置，我们采取以下步骤（如图11所示）：

1. 将速度值裁剪至预训练数据集的取值范围，并将速度归一化到[-1, 1]区间，以确保编码-解码过程的准确性；
2. 对SEAM模型进行带指定步长的重叠分块，同时保留从上到下的预测方案。因此，区块按从左上到右下的顺序排列。在我们针对SEAM模型的实验中，部分区块将完全填充先验速度信息（浅层），部分区块将部分填充（边界处），而大多数区块则无任何信息；
3. 使用预训练的VQ-VAE对所有区块进行编码，得到每个区块的离散潜在编码；
4. 从左上角开始，使用训练好的VelocityGPT模型补全在SEAM模型中遇到的第一个不完整区块；
5. 利用前一个区块中的预测速度（以其离散潜在形式）填补下一个区块（右侧）和下方区块重叠部分的缺失信息（见图11）；
6. 重复步骤4-5，直至所有区块都被处理完毕；
7. 对每个区块的离散潜在编码独立解码，将其转换回速度模型区块；
8. 重新排列区块以构建完整的速度模型，其中重叠区域进行相应的平均处理。

图 11. VelocityGPT 用于更大速度模型的推断流程示意图。按照从左上方到右下方的顺序，相邻 patch 的离散潜码根据步长及其在真实空间中的位置被复制到下一个 patch。当所有 patch 的离散潜码完成后，再将其变换回速度域。本图仅为演示，实际使用的步长比图中更小。

  我们选择64×64的区块大小（与OpenFWI速度模型尺寸一致）和16×16的步长。此外，我们仅将SEAM Arid模型的顶部64个样本作为先验，其余部分设为空（图12b）。最终，共获得833个SEAM Arid速度模型的区块。

  第一个测试仅使用类别作为每个区块的条件。图12c展示了所有区块均以FlatVelA类为条件生成的速度模型样本。可以观察到，层位连续性良好，且速度随深度总体呈增加趋势。不足之处是，由于部分特征不在模型学习的分布范围内，生成的速度模型分辨率有所降低，且出现了一些伪影。

  在第二个测试中，我们同时使用类别和井信息作为条件。假设在速度模型的X=416位置有一口井。为整合井速度信息，我们在井位处（X=416）将模型分割，翻转左侧部分，独立预测左右两部分，最后将它们拼接起来。之所以执行这些步骤，是因为我们方法的预测机制遵循特定方向；因此，为确保井信息得到充分利用和适当传播，我们从井位开始预测，然后向两个方向延伸。图12d展示了第二个测试的结果示例。可以看出，生成的速度模型在层理特征上更接近真实模型（图12a）。

  最后，在第三个测试中，我们使用与第二个测试相同的类别和井条件。但此次，我们对位于中间的区块使用不同类别，以展示我们仍能对生成模型中的地质特征进行局部控制。图12e展示了在这种设置下生成的速度模型示例，其中中间层以FlatFaultB类为条件，井位于X=0处。我们再次观察到，生成模型的层理特征与真实模型（图12a）吻合良好，且中间层（Z=100-200）还具有额外的断层特征。

图 12. (a) SEAM Arid 模型。(b) 用作先验的浅层区域。 (c) 以 FlatVelA 类别为条件生成的样本。(d) 以 FlatVelHR 类别并在 X=416 处给定一口井为条件生成的样本。(e) 以 FlatVelHR 类别（顶部与底部）、FlatFaultB 类别（中部）并在 X=0 处给定一口井为条件生成的样本。

五、讨论

  本文的主要目的是验证采用自上而下的GPT实现方式构建速度模型的可行性——该方式可将约束性更强的浅层部分作为先验来预测深部。作为概念验证性论文，我们使用了小型速度模型和简单示例。向量量化变分自编码器（VQ-VAE）作为一种压缩（潜在表示）工具，是将此方法扩展到更实用的速度模型乃至三维模型的关键。不过，我们现将该方法面临的一些挑战分享如下。从图6至图10的结果中，我们观察到了一些不足。一个显著缺陷是速度模型的生成效果不够理想。由于GPT仅处理速度模型的分布存储，任何与分辨率相关的问题都可归因于VQ-VAE的编码/解码过程。VQ-VAE通过混合损失函数（公式1）进行优化，该函数会对重建误差和潜在空间大小进行惩罚。因此，两者之间存在权衡关系，这取决于权重值（即β），而β的选择并非易事：当β值较高时（即对承诺损失赋予更大权重），重建质量可能下降；反之，当β值较低时（即对重建损失赋予更大权重），潜在维度可能任意增长。这些现象可能导致量化网络中著名的“码本坍塌”问题，即嵌入层未得到最优利用（Takida等人，2022）。在我们的应用中，β值是通过网格搜索法在VQ-VAE原始研究（Van Den Oord等人，2017）推荐的范围内（0.1–2.0）确定的，最终发现β=0.1已足够。已有研究提出了多种优化策略来解决重建-正则化平衡及β值手动调参的问题（如Mentzer等人，2023；Takida等人，2022），这些策略可用于未来版本的VelocityGPT。

  此外，当前VQ-VAE模型的架构较为简单，仅由堆叠的卷积层组成，这也是导致速度模型重建不准确的原因之一。图13a至13e展示了同一速度模型样本在不同深度范围的重建结果。在所有情况下，VQ-VAE模型的整体重建误差均可忽略不计，但在界面处除外。模型右侧背斜峰的重建效果存在明显差异——该特征越靠近边缘，重建效果越差。这种边缘效应会导致底部特征的离散表示不一致（图13f），进而可能在采样阶段造成边界附近的误差（图14）。不过，得益于VelocityGPT的自回归机制（可捕捉所有先前上下文），边缘效应并未降低生成样本的整体质量和连续性。这一问题可通过采用重建能力更强的先进量化网络（如Yu等人，2021）来缓解。与自然图像的应用场景不同，在地震领域，边缘（即层界面）的保留至关重要，因其能提供地质体深度的准确信息。探索可保留速度模型和叠后图像界面细节的技术（如Lin等人，2023提出的方法），也是改进VelocityGPT的一个潜在方向。

图 13. 使用不同深度范围的 VQ-VAE 重建与编码测试。(a–e) 深度范围分别为 Z = 8、Z = 16、Z = 32、Z = 48、Z = 64；(f) 对应 (a–e) 从左至右的离散潜码可视化。

图 14. 无条件采样测试：采用不同的起始深度。自上而下各行对应起始深度分别为 Z = 8、Z = 16、Z = 32 与 Z = 48。

  尽管VQ-VAE存在上述局限性，但其显著优势之一是对噪声的鲁棒性。我们进行了相关实验：向输入图像添加不同强度（0.25σ、0.5σ和σ）的高斯白噪声（图15）和高斯色噪声（图16）。图15b至15e及图16b至16e展示了这些图像对应的生成样本。可以观察到，即使输入图像的噪声强度高达0.5σ（部分结构因噪声过大而难以分辨，见图15a和16a），生成模型的结构仍具有合理的可信度，且与清晰图像的结构一致。需要说明的是，在训练图像VQ-VAE（VQ-VAEᵢ）时，我们既未引入噪声，也未对图像进行数据增强。这种鲁棒性源于量化过程：编码器将连续输入映射到有限的离散潜在编码集，从而有效滤除高频噪声。通过在压缩的潜在空间中运算，模型对不符合所学分布的变化敏感性降低。此外，解码器的重建能力（经训练可基于这些离散潜在编码生成输出）确保了最终输出中噪声引起的伪影被最小化。这一特性使VQ-VAE特别适用于处理含噪声的地震数据——这是实际陆地和海洋数据采集场景中常见的挑战。

图 15. 使用叠加白高斯噪声图像的采样测试。(a) 从左至右依次为：干净图像、以及叠加 0.25σ、0.5σ 和 σ 白高斯噪声的图像；(b)–(e) 分别对应使用 (a) 中各图像生成的样本。

图 16. 使用叠加色高斯噪声图像的采样测试。(a) 从左至右依次为：干净图像、以及叠加 0.25σ、0.5σ 和 σ 色高斯噪声的图像；(b)–(e) 分别对应使用 (a) 中各图像生成的样本。

  另一个问题是如何控制每个先验/条件的贡献度，这将为速度模型的生成控制提供更大灵活性。例如，当VelocityGPT被用作迭代反演过程的先验时，此功能十分有用：在迭代反演方案（如全波形反演FWI）中，早期迭代通常形成速度模型的整体结构，此时主要以类别作为条件较为合适；而在后期更新中，需恢复更精细的结构和层理细节，此时井条件可能发挥更大作用。我们的方法在多模态条件生成能力方面面临的一个挑战是：当采样过程中施加的条件存在冲突时，模型难以生成可解释的速度模型。例如，如图9所示，当地质类别条件指定为平层（如FlatVelA或FlatVelB），但施加的地震图像却显示具有倾斜和断层特征的复杂结构时，或者当井速度剖面与整体类别或结构信息相矛盾时，模型难以平衡这些不相容的约束。然而，图17中使用单一条件（类别、井或图像）的测试表明，生成的速度与各自条件吻合良好，体现了模型遵循单一约束的能力。这主要是由于训练数据中缺乏此类组合样本，限制了模型对这些场景的泛化能力。尽管模型会尝试协调这些条件，但往往生成的样本偏离地质合理性或所施加的先验。解决这一局限性需要要么在训练数据集中纳入冲突条件以增强模型的泛化能力，要么在推理过程中引入一种机制，基于预定义的优先级对某些条件进行加权。探索此类策略有助于提高未来版本模型的稳健性。

图 17. 单独条件采样测试 。(a) 无条件采样； (b) 类别条件采样； (c) 井条件采样（井位于 X=32）； (d) 图像条件采样。 所有条件均以右列所示的示例为基准。

六、结论

  我们开发了一个名为VelocityGPT的框架，该框架能够以自回归方式从浅地表到深部生成速度模型，并以浅层部分作为引导。通过使用VQ-VAE，VelocityGPT得以在速度模型的降维表示上进行运算。我们对VelocityGPT进行了训练，使其在生成过程中能够接受包括地质类别、井信息和可用图像在内的条件约束。因此，本文展示了如何在训练过程中整合速度模型的不同先验信息——包括地质特征（作为类别）、测井速度和地质构造信息（作为叠后图像），从而在推理阶段实现条件采样。我们验证了该网络生成符合注入先验信息的速度模型的能力。此外，我们还提出了一种通过基于区块的生成方式将VelocityGPT扩展到更大规模速度模型的机制，这使我们能够处理超出训练数据尺寸的分布外样本。

七、数据可用性声明

  本研究使用的数据集来源于Deng等人（2022）的研究，该数据集可公开获取。用于复现本研究实验的数据及代码，可参见Harsuko等人（2024）的相关资源。