0 引言

地震波传播的正演数值模拟是地震勘探开发中重要研究的课题。波动方程法能完整表征波在传播过程中振幅、相位、频率的变化，同时保持了地震波的运动学与动力学特征，被广泛应用于地震正演模拟。有限差分法由于适用于大规模并行计算和易编程等优势，是求解波动方程的常用方法。然而用离散的差分算子代替连续的微分算子，产生的误差会存在于整个数值模拟记录中，导致数值频散，影响正演模拟精度和逆时偏移成像效果。因此，频散压制是提高数值模拟精度的重要手段。目前数值频散压制方法从原理上大体分为以下三种。

(1) 从数值频散产生的机理出发，使用较长的差分算子或精细的差分网格进行频散压制。Alford等^[1]和Dablain^[2]对声波二阶差分数值频散进行分析，指出网格尺寸和地震波传播方向是影响频散的两个因素；董良国等^[3-4]指出高阶有限差分法可以较好解决数值频散问题；裴正林等^[5]利用一阶速度—应力弹性波方程在计算域内采用高阶有限差分法高精度模拟了弹性波场；Boore^[6]的研究表明，当震源信号波长与网格间距的比值超过10时，频散显著减弱；孙林洁等^[7]根据频散关系改变网格剖分方式，提高了有限差分法模拟的精度。此类方法会使计算复杂度和内存需求显著增加。

(2) 通量校正传输(Flux-corrected Transport，FCT)法。Fei等^[8]首次在弹性波的数值模拟中运用FCT技术，发现该方法对数值频散有较好压制效果；Book等^[9]在求解声波方程时应用FCT法有效压制了数值频散；中国学者也对FCT法消除数值频散现象方面进行了研究、分析^[10-13]。FCT技术压制数值频散的基本思想是先认为全部的极值点是数值频散形成的，然后通过对全部的网格节点都实行漫射校正，从而实现数值频散的压制。此类方法虽然减轻了计算量，但会降低有效信号的分辨率，造成振幅损失。

(3) 优化波动方程或差分系数。Etgen^[14]使用高斯—牛顿迭代方法确定有限差分系数；张志禹等^[15]在常规差分方程中增加了频散校正项，有效衰减高波数成分、抑制频散；Liu^[16]应用最小二乘法获取给定波数范围内的优化有限差分系数。还有许多优化差分系数求取方法^[17-21]，在压制数值频散方面取得了较好的效果。此类方法是多参数优化问题，需要凭借丰富的经验选取所优化的参数值，用较长计算时间确定系数且不能保证所得参数为全局最优解。

上述频散压制方法在一定精度要求下，往往导致计算复杂度大幅度增加，耗时过长，因此探究新的数值频散压制方法具有重要意义。随着硬件能力的提高和大数据技术的成熟，深度学习神经网络广泛应用到各个领域。深度学习可以从大数据中自动学习样本特征，从而表现出显著优势，在地球物理领域也得到了广泛的应用^[22-28]。Kaur等^[29]首次提出利用cycleGAN网络去除有限差分波外推中的频散伪影。该方法通过深度学习创建一个传播后置频散校正滤波器，能够以自适应的方式消除频散伪影，且与传统方法相比效果明显较好，但GAN网络复杂度较高，训练难度大，同时仅考虑到了空间域特征提取，忽略了波场数据的稀疏特性，波场的主要特征提取不充分。

本文结合波场数据的特点，参考卷积神经网络在地震资料处理方面的成功应用，提出了一种基于卷积神经网络压制频散的网络模型。该模型主要特征如下：

(1) 对不同速度模型进行正演模拟，分别使用低阶和高阶有限差分对波动方程求解，得到大量含频散输入数据集和干净波场数据集作为网络的训练输入和标签；

(2) 构建联合学习网络模型，通过空间—波数域联合损失误差约束进行多维度的特征提取和参数学习，增加描述波场特征的准确性，提高网络模型的波场特征学习能力；

(3) 考虑到波场数据的稀疏特征，在网络模型中引入稀疏正则化约束，将频散压制问题转化为L₁范数约束的最优化问题求解，降低网络模型复杂度，提高网络模型的泛化能力；

(4) 利用卷积神经网络压制数值频散，避免传统方法的计算困难和复杂度，自适应提取波场数据浅层的像素级特征和深层的语义级特征，与FCT方法相比，无需人为设定漫射系数等相关参数，技术人员无需具有丰富的地震波动力学理论基础，可直接应用预训练好的模型，实现“端到端”的处理模式。

本文首先通过均匀速度模型验证网络内部各结构的合理性；其次，通过Marmousi模型测试该算法在复杂模型中的频散压制能力；最后，通过与迁移学习结合的方式验证该算法对于新模型正演数据的频散压制泛化性能。结果表明所提算法能够保护有效信号，降低计算成本，有效地压制有限差分正演模拟过程中产生的数值频散。

1 数值频散分析 1.1 数值频散

在有限差分法求解波动方程的过程中，用离散的差分算子代替连续的微分算子，使波传播的相速度减小，不同波数分量的波场以不同传播速度进行传播，波前形状发生变化并逐渐扩散导致波形失真，造成频散假象，影响数值模拟的精度。

以各向同性介质中的二维波动方程

$\frac{1}{{{v^2}}}\frac{{{\partial ^2}U}}{{\partial {t^2}}} = \frac{{{\partial ^2}U}}{{\partial {x^2}}} + \frac{{{\partial ^2}U}}{{\partial {z^2}}}$

(1)

为例分析数值频散产生的机制。式中：U为随空间位置和时间变化的波场；v为速度；t为时间；x和z分别为空间上的水平、垂直方向坐标。

式(1)的平面波解可表示为

$p(x, z, t) = \exp [{\rm{i}}(\omega t - kx\cos \varphi - kz\sin \varphi )]$

(2)

式中：k为波数；φ为入射角；ω=kv₀为角频率，其中v₀为数值相速度。通过推导可得v₀与v的关系为

$\begin{array}{l} \frac{{v_0^2}}{{{v^2}}} = - \frac{1}{{{k^2}{h^2}}} \times \\ \;\;\;\;\;\;\;\;\;\;\sum\limits_{m = - M}^M {{c_m}} [\cos (kmh\cos \varphi ) + \cos (kmh\sin \varphi )] \end{array}$

(3)

式中：M为有限差分阶次的一半；c_m为有限差分系数；h=Δx=Δz，为网格间距。

由式(3)可以看出：①频散误差可基于波传播的相速度和群速度加以描述，波数越大，相速度滞后越多；②差分精度越高，v₀²/v²比值越接近1，波场外推精度越高。

以均匀介质中的单炮波场快照为例，监测区域横、纵方向均设置50个网格，地下介质速度为2000m/s，在区域中心选取主频为40Hz的Ricker子波作震源，网格尺寸为8m×8m，不同精度差分算子模拟的1000ms波场快照如图 1所示。由于差分阶数较低，数值频散造成的干扰波逐渐向有效波场内部扩散，使波场快照的分辨率和信噪比降低(图 1a)；随着差分算子精度的提高，数值频散逐渐减弱(图 1b、图 1c)。地震有效波形所占检测区域的比重较少，因此在空间域内波场数据具有一定的稀疏性。

图 2为图 1在x=200m处抽取的单道波形数据。低阶差分算子模拟的波形曲线中出现严重的振铃效应，干扰了有效波形信息；随着差分算子精度的提高，所得波形趋于光滑，数值频散造成的伪振动被压制，更接近于输入的Ricker子波波形。可见通过提高差分算子精度的方式可有效压制波场模拟中的数值频散。

1.2 波场数据的波数域特征

傅里叶变换是地震数据分析的主要数学工具。对有限差分法波场数值模拟结果，利用傅里叶变换将波场快照从空间域变换至波数域，降低了波场信号的相关性，能量集中在少数波数成分上，实现对信号的稀疏表示，使得在波数域中更易找出数值频散的模式。二维傅里叶变换的形式为

$\widetilde U\left( {{k_x}, {k_z}} \right) = \iint{{{U(x, z)}\exp \left[ { - {\rm{i}}2\pi \left( {{k_x}x + {k_z}z} \right)} \right]{\rm{d}}x{\rm{d}}z}}$

(4)

式中：k_x、k_z分别为x和z方向的波数；$ \widetilde U\left( {{k_x}, {k_z}} \right)$为变换后的波数域波场。

图 3为不同阶次差分算子模拟的波数域波场。由于波数较高的地震波相速度较慢，滞后于群速度形成“波拖尾”现象，因此当不同差分精度算子所得波场快照变换至波数域时，含频散波场数据傅里叶变换域内的系数向高波数分量扩散，而干净波场数据波数域内信息分布则更为集中。

2 数值频散压制网络模型 2.1 网络模型

本文应用残差学习的思想设计基于监督学习的深度卷积网络。训练数据以输入—标签对的形式输入网络

$\left\{ {\left( {{{\mathit{\boldsymbol{x}}}_1}, {{\mathit{\boldsymbol{y}}}_1}} \right), \left( {{{\mathit{\boldsymbol{x}}}_2}, {{\mathit{\boldsymbol{y}}}_2}} \right), \cdots , \left( {{{\mathit{\boldsymbol{x}}}_{{N_{\rm{d}}}}}, {{\mathit{\boldsymbol{y}}}_{{N_{\rm{d}}}}}} \right)} \right\}$

(5)

式中：N_d为波场数据总量；x为含频散波场数据；y为对应的干净波场数据，即标签数据。通过残差学习策略提取数值频散产生的伪波纹理及分散的特征。神经网络的输出为

$\widehat {\mathit{\boldsymbol{y}}} = {\mathit{\boldsymbol{x}}} - {\mathop{\rm Res}\nolimits} ({\mathit{\boldsymbol{x}}};{\mathit{\boldsymbol{\theta }}})$

(6s)

式中：$ \widehat {\mathit{\boldsymbol{y}}}$表示网络压制频散后波场数据；Res(·)代表残差学习；θ表示网络参数。通过训练使网络预测结果$ \widehat {\mathit{\boldsymbol{y}}}$与样本标签y近似。

网络模型结构如图 4所示。本文网络结构由多层的卷积单元(Convolutional Layer，Conv)、批归一化(Batch Normalization，BN)单元以及激活函数单元组成。选择线性整流(Rectifier Linear Unit，ReLU)单元作为激活函数进行非线性映射，可表示为

$f(l) = \max (0, l)$

(7)

式中l表示激活函数的输入。ReLU函数具有单侧饱和特性，在负半区的导数为零，减少激活的矩阵中非0元素个数，从而在网络内部引入稀疏性，避免无用信息的干扰。

第1层为输入层，由Conv +ReLU单元组成，其中Conv用于提取含频散数据的特征，ReLU用于执行非线性映射，去除含频散数据中的冗余特征；第2~第19层为网络的隐藏层，由交替的Conv +BN+ReLU单元组成，用于完成数据处理后特征提取的非线性映射；前20层卷积核大小设为3×3，步长为1，每一层卷积后得到64张特征图；最后一层为Conv单元，将上一层完成的非线性映射输出作为输入进行特征提取，在前向传递过程中每进行一次卷积后均进行补零操作，保证训练过程中每个特征图与输入数据尺寸一致。

图 5为网络压制频散的过程，在空间域训练含频散波场数据后，对输出的波场数据与标签分别进行傅里叶变换，得到波场数据的波数域特征；对空间域和波数域上的波场信息分别进行误差计算后得到联合误差函数，通过使联合误差函数最小化进行网络参数更新。

2.2 稀疏正则化联合损失函数的建立 2.2.1 联合损失函数

针对频散特性，使用联合损失函数从空间—波数联合域上对预测波场与目标波场的差异进行多维度描述，联合损失函数的形式为

$L = \frac{\mu }{{{N_{\rm{d}}}}}\sum\limits_{n = 1}^{{N_{\rm{d}}}} {\left\| {{{\mathit{\boldsymbol{y}}}_n} - {{\widehat {\mathit{\boldsymbol{y}}}}_n}} \right\|_P^P} + \frac{\upsilon}{{{N_{\rm{d}}}}}\sum\limits_{n = 1}^{{N_{\rm{d}}}} {\left\| {{\mathit{\boldsymbol{y}}}_n^{\rm{f}} - \widehat {\mathit{\boldsymbol{y}}}_n^{\rm{f}}} \right\|_P^P} $

(8)

式中：$ {\mathit{\boldsymbol{y}}}_n^{\rm{f}}$为波数域样本标签；$ \widehat {\mathit{\boldsymbol{y}}}_n^{\rm{f}}$为波数域输出；μ、υ为空间、波数域的平衡因子，通过调整这两个参数可以控制不同域之间误差的权重；‖·‖_P^P表示范数约束，若P=1，使用L₁损失函数；若P=2，则用L₂损失函数。

在波场数据频散压制的任务中，数值频散形成的伪波纹理并不是均匀分布在整个波场快照中，而是集中在有效波形内部附近，且强度较小。神经网络回归问题中常用L₂损失函数，但通常情况下适用于对高斯噪声进行建模，对强度较大的噪声抑制效果较好，不适用于对波场数据强度较小的频散纹理进行压制；L₁损失函数对估计值$ \widehat {\mathit{\boldsymbol{y}}}_n$与真实值y_n之差取绝对值，与L₂损失函数相比更具鲁棒性。本文选取L₁范数作为损失函数度量空间域和波数域的误差。目标函数的梯度为

$\frac{{\partial L}}{{\partial {{\widehat {\mathit{\boldsymbol{y}}}}_n}}} + \frac{{\partial L}}{{\partial \widehat {\mathit{\boldsymbol{y}}}_n^{\rm{f}}}} = \mu {\mathop{\rm sign}\nolimits} \left( {{{\mathit{\boldsymbol{y}}}_n} - {{\widehat {\mathit{\boldsymbol{y}}}}_n}} \right) + \upsilon{\mathop{\rm sign}\nolimits} \left( {{\mathit{\boldsymbol{y}}}_n^{\rm{f}} - \widehat {\mathit{\boldsymbol{y}}}_n^{\rm{f}}} \right)$

(9)

式中sign(·)为符号函数。由于L₁范数在0处不可导，设置sign(0)=0，代替不连续区域的梯度值。

2.2.2 稀疏约束正则项

本文联合损失函数仅可保证模型取得较小的训练误差，网络的复杂程度和泛化能力需要引入稀疏约束项。稀疏约束项可以理解为对网络复杂度的约束，以减少网络中的参数冗余，提高模型的泛化性能。由于波场数据在空间域和波数域均具有稀疏性，样本中有效信号的比重较小，引入稀疏正则约束后则可以较好收缩模型、减小解空间，提高网络预测的准确性和鲁棒性。因此本文通过在联合损失函数中加入模型的稀疏约束，达到利用先验信息压制波场频散的目的。稀疏约束项为

$R({\mathit{\boldsymbol{\theta }}}) = ||{\mathit{\boldsymbol{\theta }}}||{_P}$

(10)

式中‖θ‖_P表示网络权值的L_P正则项。

L₀正则项从直观上是对网络非零参数的个数进行约束，可实现网络参数对模型的有效表达得到最稀疏的解，但L₀范数非凸且不可微，很难对其进行优化求解。L₁正则化是L₀正则化的最优凸近似，L₁范数作为惩罚项在训练过程中对权重的更新进行稀疏约束，使模型的大量参数近似为0，实现网络模型参数的稀疏化和特征选择优化。对于卷积核的权重θ′={θ′ ₁, θ′ ₂, …，θ′ _J}，其中J为卷积核中参数的数量，稀疏约束表达式为

$R\left( {{\mathit{\boldsymbol{\theta }}^\prime }} \right) = \sum\limits_{r = 1}^J {\left| {\mathit{\boldsymbol{\theta }}_r^\prime } \right|} $

(11)

则添加正则约束的最终损失函数为

$\begin{array}{l} L = \frac{\mu }{{{N_{\rm{d}}}}}\sum\limits_{n = 1}^{{N_{\rm{d}}}} {{{\left\| {{\mathit{\boldsymbol{y}}_n} - {{\mathit{\boldsymbol{\widehat y}}}_n}} \right\|}_1}} + \\ \;\;\;\;\;\;\frac{\upsilon }{{{N_{\rm{d}}}}}\sum\limits_{n = 1}^{{N_{\rm{d}}}} {{{\left\| {\mathit{\boldsymbol{y}}_n^{\rm{f}} - \mathit{\boldsymbol{\widehat y}}_n^{\rm{f}}} \right\|}_1}} + \lambda R\left( {{\mathit{\boldsymbol{\theta }}^\prime }} \right) \end{array}$

(12)

式中λ为稀疏约束项系数。当λ→0时，表示不对网络权值进行约束，只通过样本数据学习权值；当λ→∞时，表示只根据惩罚项确定网络的权值。

为了平衡空间域与波数域损失的权重，设置μ∶v=1∶5，λ=5×10^-5，采用自适应矩估计(Adaptive Moment Estimation，Adam)算法对损失函数进行最优化学习。

3 模型实验及结果

分别选取均匀速度模型以及Marmousi速度模型对本文算法进行训练和测试。对于每一个速度模型，设置不同的网格尺寸、主频、采样间隔等正演参数，选取Ricker子波作为震源，分别采用低阶差分和12阶差分进行波场模拟，所得到的波场数据作为含频散的训练数据集和标签数据集。

频散压制效果可采用信噪比(SNR)衡量，其数学定义为

${\rm{SNR}} = 10\lg \frac{{\sum\limits_{i = 1}^{{N_x}} {\sum\limits_{j = 1}^{{N_y}} {{{\left( {{y_{i, j}}} \right)}^2}} } }}{{\sum\limits_{i = 1}^{{N_x}} {\sum\limits_{j = 1}^{{N_y}} {{{\left( {{y_{i, j}} - {{\hat y}_{i, j}}} \right)}^2}} } }}$

(13)

式中N_x、N_y为波场切片的x和z方向上的样点数。

实验选用GPU为NVIDA GTX-2080的平台，操作系统为64位Ubuntu 18.04，软件环境为Python 3.8，联合深度学习模型使用Pytorch1.6框架搭建，CUDA版本为11.0。

3.1 均匀速度模型实验

为了验证模型的泛化性能，通过不同正演模拟参数数值模拟构建数据集A：监测区域横向和纵向分别设置为200个网格，采用完全匹配层(Perfectly matched Layer，PML)吸收边界条件在波场外围布置100个网格的吸收衰减层避免边界反射，时间步长为0.25ms。训练数据集A所用正演参数如表 1所示，对每种正演参数数值模拟获得含频散数据和样本标签数据，对每一个样本、标签进行傅里叶变换，得到波数域的样本、标签。将数据集A分别按照8∶1∶1的比例划分为训练集、验证集和测试集，对网络模型各部分的有效性进行验证。

3.1.1 学习率更新对网络的影响

学习率是深度学习中重要的超参数之一，由于L₁损失在网络更新时的梯度始终相同，在接近最优值处可能维持较大的梯度而错过最优值，因此本文在训练迭代过程中自动调节学习率：当验证集的损失函数在3次迭代后不再下降时，将学习率调整为当前学习率的0.5倍。在训练过程中初始学习率设为0.003，批大小为8，迭代次数为400，学习率的更新与SNR的变化关系如图 6所示，分别在迭代64、92、124、156、200、254、285、316、348、380次时更新学习率，使目标函数重新收敛到全局最小值，得到最优解。

3.1.2 稀疏约束正则项效果分析

将不含稀疏约束项的网络模型记为算法Ⅰ，本文算法和算法Ⅰ对验证集的SNR变化曲线如图 7所示，可以看出两个网络模型在验证集性能上存在较大差别，原因在于算法Ⅰ在训练过程中网络过分拟合训练数据的特征，导致在未知样本集上表现不佳；本文算法在验证集上的SNR接近于训练集，提升了网络的泛化能力。

在神经网络提取特征的学习过程中，浅层卷积核提取纹理等细节特征，较深层次卷积核提取波场数据的轮廓、形状等抽象语义特征。为了进一步说明稀疏约束项的效果，对比两种算法隐藏层中较深层次卷积核权重分布如图 8所示。由图 8a可以看出，结合稀疏约束项的卷积核所学习的参数特征分布基本一致，表明本文算法与算法Ⅰ相比，在特征学习时可较为稳定地提取波场数据的深层语义信息；算法Ⅰ中卷积核学习的权重较为分散(图 8b)，由于对卷积核施加L₁正则约束，本文算法卷积核权重分布服从拉普拉斯分布，更多较小权重被惩罚近似为0，实现特征优化选择。综上分析可知加入稀疏约束的网络模型具有一定的合理性。

3.1.3 空间—波数域联合误差效果分析

将只含空间域误差与稀疏约束项的模型记为算法Ⅱ，只含波数域误差与稀疏约束项的模型记为算法Ⅲ。分别使用算法Ⅱ、算法Ⅲ和本文算法对波场数据进行频散压制，均进行400次迭代，不同算法网络性能随迭代次数的变化如图 9所示。从SNR随迭代次数变化曲线可知，算法Ⅱ由于未结合波数域误差，在网络学习过程中未能充分提取特征，结果地SNR低于本文算法。算法Ⅱ、算法Ⅲ和本文算法进行400次迭代后SNR的值分别为22.87dB、6.72dB与34.82dB，本文算法的SNR高于使用单一空间域或波数域误差作为损失函数的算法(图 9)。从算法Ⅲ中第400次迭代中任取一个样本展示其压制频散效果，由于仅含波数域计算目标函数，能量没有被恢复，所提取特征有限导致无法根据现有特征重构出完整波形，频散纹理仍存在于有效波形附近且波场快照分辨率较低(图 10)。

3.1.4 残差学习效果分析

本文所提出模型在训练过程中对频散形成的伪波纹理及伪振动所造成的能量损失进行学习，之后由输入数据与残差数据相减得到压制频散后的干净波场。

将本文方法与不含残差学习策略的模型(记为算法Ⅳ)进行对比。由二者SNR随迭代次数变化的曲线(图 11a)可见，本文算法比算法Ⅳ更容易优化，且网络收敛更快(图 11b)。从总体压制频散效果来看，引入残差学习策略的模型优于直接映射训练的网络模型。

3.1.5 不同损失函数效果分析

损失函数的选择对网络性能有很大影响。将使用L₂范数为损失函数的模型记为算法Ⅴ，将该算法与本文算法进行对比测试。训练过程中其余参数保持不变，两种算法均迭代400次。由图 12可以看出，算法Ⅴ训练时的SNR曲线较为平滑，收敛时的SNR远低于本文模型，这是由于算法Ⅴ使用L₂损失度量预测值与真实值差异时将误差平方化，在进行频散压制的任务中会赋予伪波纹理较大的权重，以牺牲有效信号的精度为代价，朝着减少伪波纹理的方向更新网络，降低了模型的整体性能。

3.1.6 不同网络结构模型测试

将上述训练好的网络模型保存，使用与训练数据不同震源位置的波场切片对上述网络模型进行测试。图 13为不同算法测试后的波场快照SNR统计，可见本文所提出的算法在测试数据上具有更高的信噪比，使用结合稀疏约束的联合损失函数有效降低了传统损失函数对波场重构带来的负面影响，在频散压制效果上有明显优势。

Kaur等^[29]使用生成对抗网络(GANS)中的cycleGAN网络压制有限差分法中的数值频散。将含频散的波场数据记为源域，无频散伪影的波场数据记为目标域。该网络的目的是将数据由源域的特征转换为目标域特征，采用了循环一致性损失等多种损失函数来增强由源域至目标域的映射约束。为了比较两种网络对数值频散的压制效果，选用与本文均匀速度模型相同的数据集进行对比实验。cycleGAN网络参数为：将cycleGAN网络生成器部分调整为9个残差块，网络初始学习率设为0.0002，每迭代50次调整学习率，批大小数量设为8。两个网络模型均迭代400次。为了进行两个模型泛化性能的对比，将训练好的网络模型保存，使用与训练数据震源位置不同的波场切片对两种网络模型进行测试，结果如图 14所示。

由两种网络模型测试结果(图 14)可以看出，频散现象均被压制，但cycleGAN网络压制结果的波形边缘不清晰，影响了波场快照质量。由于cycleGAN网络在压制过程中仅仅关注于波场数据的空间域特征，无法获取更多维度特征去加强对目标域的映射，转化后的波场快照质量与样本标签存在一定差距；从运行时间上看，由于GAN网络复杂度较高，因此训练难度较大，迭代400次所需时间为137h42min，而本文算法训练耗时为6h3min，本文网络模型的运行时间远小于cycleGAN。

3.1.7 模型的适用性分析

为了验证本文模型的适用性能，选用不同正演参数的波场数据进行测试(表 2)。SNR统计结果表明本文网络对未经训练过的数据频散压制效果明显。测试数据集d和测试数据集e的处理结果如图 15所示，可见，数值频散形成的伪波纹理被充分压制，得到较为清晰的波形边缘，说明本文网络具有一定的泛化能力和适用性，不仅仅在训练数据上有着很好的拟合效果，当模型不变，正演数值模拟参数在一定范围内发生变化时，无需重新训练即可获取高质量的波场数据。

3.2 复杂Marmousi模型

选取Marmousi速度模型(图 16)验证本文方法在复杂波场情况下的频散压制能力。模型网格数为298×300，PML吸收边界厚度设为100个网格点。网格尺寸为8m×8m，时间步长为0.25ms。由于Marmousi模型中包含丰富的反射界面与速度信息，为使训练数据包含波场信息更丰富，分别在(800m，800m)和(1200m，1200m)处设置主频为40Hz的Ricker子波作为震源，两处震源的波场快照按1∶1的比例制作成训练数据集B，并按照8∶1∶1的比例划分为训练集、验证集与测试集。

使用本文算法对训练数据集B进行压制频散处理，400次迭代后一个样本的压制频散结果如图 17所示，可以看出：训练输入样本出现很多同相轴，频散形成的伪波与反射波混叠，波形边缘模糊且分辨率很低，数值频散现象严重；经过频散压制后波场快照的分辨率有所提升，去除了大量的伪波纹理的同时保留了有效信号，具有较高的保真性。

抽取图 17在x=720m处波形曲线(图 18a)，并计算波数谱(图 18b)。由图 18a可以看出，频散波波形经过网络频散压制后，输出的波形与样本标签接近一致，有效去除了频散干扰。由图 18b可以看出，波场中低波数部分的能量对应着真正的波场值，数值频散形成的能量集中在高波数区域，输入数据频散现象较为严重，在没有振动的高波数成分上也产生振动，形成了虚假绕射波；频散压制后数据的高波数分量上的幅值低于输入数据，意味着频散能量被有效压制。

将训练好的网络模型保存，保持其他参数不变，使用与训练数据不同时刻的波场切片对模型进行测试，为对比测试前后的SNR，同样利用正演记录获得测试数据的样本标签，频散压制前、后的空间、波数域如图 19所示。

由图 19a可以看出，频散现象较为严重，信噪比较低；经过本文模型频散压制后(图 19b)，波场切片变得清晰，对于差分精度低造成的能量损失也得到了较精确的补偿，消除了频散带来的干扰，在空间域内频散所形成的伪波被去除、在波数域内高波数分量的系数被压制，能量更为集中。

为进一步对比本文所提方法的压制频散效果，将网络模型与本文相同、使用MSE为损失函数且不含波数域误差和稀疏约束项的模型记为算法Ⅵ。图 20a为算法Ⅵ的压制频散结果，图 20b为本文算法的输出结果，图 20c和20d分别为两种算法输出结果与标签数据的残差。由于算法Ⅵ仅提取了波场的空间域特征进行学习，且未引入稀疏约束项，因此频散伪波纹理学习不够充分，少量频散伪波仍存在于压制后的数据中，波形整体边缘模糊，频散压制效果较差；从残差剖面上看，本文压制频散结果与标签数据的残差更小，对频散造成的波形和能量损失学习得较为充分。抽取两种方法频散压制结果在z =600m处进行单道分析(图 21)。由图 21a可以看出，测试样本的振动曲线波动剧烈，算法Ⅵ压制结果仍存在锯齿状的波动，本文算法频散压制后可得到完整的Ricker子波波形。由图 21b可以看出，本文算法频散压制后，与样本标签相比振幅误差约为0.04，远小于算法Ⅵ的预测误差，说明本文算法对振幅补偿效果较好。由图 21c可以看出，本文算法压制后的波数谱曲线高波数能量低于算法Ⅵ压制结果，且充分保留了低波数分量能量，说明本文方法有效压制了数值频散且没有破坏有效信号，验证了本文算法处理复杂波场数据的有效性。

3.3 断陷模型

为了测试本文方法对不同地下介质速度模型压制频散的能力，设计一个断陷速度模型(图 22)。该速度模型横向网格数为400，纵向网格数为300，网格尺寸为10m×10m。在(2000m，1500m)处激发40Hz的Ricker子波作为震源，采用时间2阶、空间2阶和12阶的有限差分数值算法完成该断层速度模型地震波场正演数值模拟，得到新的数据集C。

结合迁移学习的思想应用该断陷模型数据集对网络进行再训练，具体实现方式如下：首先，利用多种地下介质速度模型所得正演波场数据作为源域进行训练，得到预训练网络模型；选取数据集C中的20%作为目标域的训练集，将预训练模型隐层参数赋值到新数据集模型相应层中进行权重初始化，训练迭代300次；使用目标域的80%剩余数据进行测试，验证经过迁移学习训练的网络模型效果。

图 23为测试数据中一个样本及经网络频散压制后的结果，测试样本的SNR为0.9357dB，经频散压制后SNR达到26.5088dB，可见经频散压制后，有效波形突显出来，可清晰观察到反射波波形。

将在其他正演数据集上预训练过的网络模型作为初始网络，从而避免了模型随机初始化参数时的逐步调优过程，使用少量的新的模型数据对网络进行训练，就可得到高质量波场数据，可见通过迁移学习的方式节约了训练成本、提升了频散压制的效率。

4 结论

本文提出了一种基于联合学习与稀疏约束的波场数值频散压制方法，该方法通过对波场数据的空间域和波数域特征分析定义了联合损失函数，并在网络中增加稀疏约束，达到频散压制目的。通过均匀模型、复杂的Marmousi模型以及断陷模型验证了本文所提方法的有效性、适用性与鲁棒性。获得的结论如下：

(1) 与传统数值频散压制方法相比，可以较低的计算代价完成频散压制，避免了FCT方法所带来的振幅损失或计算量过大的问题；

(2) 联合空间域和波数域的信息定义损失函数，能够更充分抑制频散形成的高波数分量，使频散特征提取更充分、对频散压制更彻底；

(3) 选取合理的稀疏正则化约束参数对卷积核权重进行约束，使网络频散压制效果提升且对于未经训练的未知样本频散压制效果较好，提升网络模型的泛化性能；

(4) 由于边缘保持能力及鲁棒性，L₁范数损失函数更适用于频散压制，结果更清晰，边缘更明显；

(5) 本文算法与迁移学习结合，对新模型的正演波场进行频散压制可取得较好的实验结果，在节约训练成本的同时使频散压制效率有所提升。

三个不同速度模型的实验结果表明，本文方法在有效压制频散的同时对有效波形损伤较小，有效保护了波形的完整性和边缘的光滑性。