引言

自古以来，地震都是一种难以避免的自然现象（龙锋等，2006）。中国作为全球地震频发的地区之一，迫切需要先进的地下探测技术来刻画岩石圈的精细结构，从而更好地了解地震的孕育与发展规律，尽可能减少其对人类社会的危害（王椿镛等，2016）。自20世纪80年代以来，全波形反演逐渐兴起（Tarantola，1984；Jeong et al，2012），进一步丰富了常规速度建模方法的理论体系（Virieux，Operto，2009）。这种反演方法凭借高精度、多参数、多维度的建模优势，引起了许多研究者的关注（张文生等，2015；Virieux et al，2017）。

作为一种高分辨率地震反演方法，全波形反演通过波场模拟与最优化迭代来极小化衡量合成波形与地震观测数据偏差的目标函数，获得最优地下介质参数并进行构造成像（Bornstein et al，2013；Yang et al，2020）。在研究的早期阶段，全波形反演主要是在时间域上进行（Tarantola，1984；Tessmer et al，1992；Hestholm et al，1999）。随后，Pratt （1990）使用频率域波动方程进行逐个频率点的迭代计算，将反演理论推广到频率域。相比之下，频率域反演方法易于并行可减轻计算负担且容易刻画衰减效应（Pratt et al，1998；Shin et al，2001）。Sirgue和Pratt （2004）提出只需要选取少量频率在保证波数充分覆盖的前提下即可刻画反演模型的细节结构，从而在不损伤反演精度的条件下进一步降低频率域全波形反演的计算代价。殷文等（2006）研究得出若频率域各频率点之间相互独立计算，则各个频率点误差不会产生累积，可极大提高模拟精度。此外，还有研究发现，频率域波场模拟对于多震源和增加衰减因子更加有效（Ben-Hadj-Ali et al，2011）。由此可见，在频率域上进行全波形反演具有重要意义（Geng et al，2018）。

随着地震学研究的不断深入与计算机硬件水平的快速提升，传统声波方程较难满足地下精细探测的需求（张霖斌，姚振兴，2000；Jeong et al，2012）。全波形反演逐步在向基于各向同性、各向异性、黏弹性介质波动方程等的多弹性参数反演方向发展（宋海斌，张关泉，1998；Operto et al，2007）。从物理的角度来讲，在短时间震源激发下，离震源作用地点较远处的地下岩层可近似地看成弹性体，此时地震波在地球介质中的传播可以看成弹性波（孙成禹，李振春，2011）。相较于声波，弹性波包括了P波与S波的传播规律，包含的波场信息更加丰富。另外，在地表区域较短的传播时间导致体波和面波不容易分离，此时过滤或减弱表面波并不容易，两种类型的波都需要参与成像。因此，弹性波方程能够更加准确地模拟地震波在地下介质中的传播情况（Brossier et al，2009；Zhao et al，2020）。相较于声波情形，弹性波全波形反演能同时获取P波、S波速度等参数，但由于不同波形之间会相互干扰，因此会产生严重的非线性问题（Virieux，Operto，2009）。到目前为止，对于弹性波反演出现的串扰问题已有很多学者进行了研究，所使用较多的方法有辐射模式法（Prieux et al，2013）、Hessian矩阵法（Wang et al，2016）、纵横波分离法（黄少华等，2019）等，其中纵横波分离法早期被运用到弹性波逆时偏移来压制串扰噪声，近几年很多学者将其运用到弹性波全波形反演中。针对多参数反演的梯度有很强的串扰问题，Ren和Liu （2016）首先提出了在波数域进行P波与S波分离的方法以降低梯度串扰从而提高反演精度，而Qu等（2018）应用矢量纵横波分离方法直接对纵横波速度进行扰动来求取P波梯度和S波梯度也取得不错的效果。

由于全波形反演的主要计算量集中在正演过程，反演的计算效率和精度在很大程度上取决于正演模拟方法的优劣。同时考虑到弹性波方程本身结构较为复杂，数值求解需要更多的计算代价，因此采用高效的正演方法至关重要（刘璐等，2013；Li et al，2018）。目前频率域计算普遍采用的正演方法包括有限差分方法（Moczo et al，2007）、有限元方法（Liu et al，2017）、谱元法（Komatitsch et al，2005）等。其中有限差分法凭借原理直观、离散算子结构简单及易于实现并行等优点而被广泛应用。相比而言，传统差分格式精度较低且容易引起数值频散假象，为克服这一缺点通常需要加密网格，这势必导致计算效率的明显下降。

近似解析离散化（nearly analytic discrete，缩写为NAD）算子是近些年来兴起的一类重要差分算子（Yang et al，2003），构造思想是利用节点位移及其梯度值的交叉组合来逼近高阶偏导数（Yang et al，2009）。这同时具备较高精度，可显著提升计算效率。Yang等（2010）、Tong等（2011）及Li等（2013）将NAD算子用于时间域声波和弹性波方程的数值离散过程并进一步优化了算子结构。Lang和Yang （2017）采用NAD算子进行了频率域声波方程正演模拟，并通过与其它经典数值方法（例如九点有限差分法（OFD）和交错网格（SG）方法）的对比，显示了该算子在压制数值频散与提高计算效率方面的优势。随后，Lang等（2020）与郎超等（2022）将NAD算子用于频率域声波全波形反演和逆时偏移计算中，均得到了可靠的成像结果。然而该算子在频率域弹性波方程反演中的讨论还处于空白阶段，相应算法的适用性尚不明确。

为此，本文试图搭建基于NAD算子的频率域弹性波全波形反演框架；利用阻抗矩阵的稀疏分块结构来推导反演目标函数的梯度计算公式，并且将分离纵横波的思想贯彻整个求解过程，充分缩减梯度计算的规模并压制纵横波之间的串扰从而进一步提高反演计算效率；最后利用三种经典介质模型进行数值实验，结合不同频率的反演误差曲线趋势图验证所提出弹性波反演方法的有效性和正确性.

1.   基于NAD算子的频率域弹性波正演模拟

在各向同性介质中的二维频率域弹性波方程可表示为

式中，$ x $表示地面横向坐标，$ {\textit{z}} $表示地下纵向坐标，$ u $，$ v $分别定义为弹性波场的水平和垂直分量，$ \mu $和$ \lambda $为拉梅系数，$ {s}_{1} $和$ {s}_{2} $为震源分量，$ \rho $和$ \omega $分别为介质密度和角频率。由于在影响地震成像相关参数的分析中，速度参数的作用更为关键（张文生，张丽娜，2020），所以在下文的正演模拟与反演迭代计算中，本文假定反演的密度值为某一常数且不再变动，只反演P波与S波速度。

为后续论述清晰明了，可将式（1）归结为矩阵与向量形式，即

其中

式中${\boldsymbol{w}} $和${\boldsymbol{s}} $分别为波场向量和震源向量。

在使用NAD算子对式（2）进行离散时，需要对其两端分别沿着$ x $和$ {\textit{z}} $轴求偏导，可得到如下偏微分方程组：

1.1   吸收边界条件

当使用计算机进行地震波正演模拟时，由于计算空间的局限性，研究者会引入人工截断边界，这导致波传播至人工边界时会引发非物理的反射而产生虚假反射波，因此必须对边界进行处理。本文采用完全匹配层（perfectly matched layer，缩写为PML）吸收边界条件，它的构造思想是在边界处加一个吸收层，数学上可采用衰减函数来表征这一衰减效应（Komatitsch，Tromp，2003）。

具体地，以$ x $方向为例，首先在PML区域内引入复坐标$ \tilde x $替换式（4）中的实坐标$ x $，即

关于$ x $求偏导，可得

式中：$ \mathrm{i}＝\sqrt{-1} $；衰减函数的具体表达式为${d_x} （ x ） ＝ - \dfrac{{3{v_x}}}{{2\delta }}\ln \left(R\right){\left( {\dfrac{x}{\delta }} \right)^2}$，其中$ \delta $表示$ x $方向PML层厚度，$ {v}_{x} $为PML区域和计算区域交界处的波速，$ x $为PML区域内的网格点到计算区域的投影距离，$ R $为反射系数（一般取10⁻³）。

将式（4）中的实坐标替换为复坐标，再将式（5）、（6）代入后（$ {\textit{z}} $方向同理可得），可得带有PML边界条件的方程组

1.2   四阶NAD算子

若式（7）中波场${\boldsymbol{w}} $在点（i，j）处的二阶、三阶偏导数离散，则所采用的$ x $方向上的单方向四阶NAD网格差分模板为（Yang et al，2003，2004）

$ {\textit{z}} $方向可相应推导.

通过单方向上的NAD模板再借助旋转坐标轴的思想（童平，2012）可推导出双方向四阶NAD差分格式为

式中$ h $表示空间离散步长（$ x $和$ {\textit{z}} $方向上的相等）。将差分模板式（8）—（11）代入式（7）进行离散，可形成大型稀疏线性方程组（郎超等，2021）

式中C为阻抗矩阵。

1.3   阻抗矩阵的稀疏分块结构

线性方程组式（12）中的矩阵$ \boldsymbol{C} $称为阻抗矩阵（Pratt，1999），$ {{\boldsymbol{w}}} $表示弹性波场的水平与垂直分量交叉排列所组成的波场向量，$ \boldsymbol{s} $为右端震源向量，$ \boldsymbol{C} $，$ \boldsymbol{w} $，$ \boldsymbol{s} $随反演模型向量$ \boldsymbol{P} $的更新而改变，同时也依赖于角频率$ \omega $。为了选择合适的算法以快速而稳定地求解式（12）并在反演过程中推导出目标函数梯度计算的有效方法，阻抗矩阵$ \boldsymbol{C} $的稀疏分块结构必不可少，其示意图如下所列（郎超等，2021）

式中，非零子块矩阵$ {\boldsymbol{C}}_{i}^{ （ j ） } （ i＝1 \text{，} \cdots \text{，} 9 \text{，} j＝1 \text{，} \cdots \text{，} N ） $的规模为$ 6{\text{×}} 6 $，其余空白位置未列出的均为零元素，$ N＝{n}_{x}{\text{×}}{n}_{{\textit{z}}} $表示正演离散网格点总数，$ {\boldsymbol{C}}_{i}^{ （ j ） } $的具体元素表达式参考郎超等（2021）。利用式（13）的稀疏结构特点，本文中相应的反演计算依然选择以往研究工作中（Lang et al，2020）所采用的三角分块预处理GMRES迭代方法对线性方程组（12）进行求解，该方法配合频率域NAD离散已展现其在提高计算效率方面的优势（郎超等，2021）。

2.   基于NAD算子的频率域弹性波反演原理

全波形反演的基本思路是通过合成波场与观测波场数据残差的欧式范数建立反演目标函数，并利用最优化方法（例如：梯度类方法）更新模型参数向量，实现观测波场与合成波场的最佳匹配，从而达到重构速度场的效果（祝贺君等，2023）。对于极小化反演目标函数的计算过程，本文采用梯度类迭代方法来实施（张文生，张丽娜，2020）。

常用梯度类算法主要包括最速下降法（Jeong et al，2012）、非线性共轭梯度法（Tong et al，2014）、高斯牛顿法（Pratt，Worthington，1988）、拟牛顿法（例如：BFGS）（Tromp et al，2005）等。在上述所有算法中，最速下降法的优点是迭代流程直观，计算操作简单，缺点是收敛速度较慢，计算效率偏低；而牛顿类方法作为二阶方法，所得反演结果的精度较高，但计算量较大，不适用于大规模计算（任志明，2016）。为了平衡精度和计算量，本文采用非线性共轭梯度方法（Saad，2003）。

根据非线性共轭梯度方法的原理，求解目标函数的梯度值是整个迭代算法不可缺少的环节（卞爱飞等，2010）。因此，本节着重推导基于NAD算子的目标函数梯度计算公式，由此实现相应的频率域全波形反演流程。相应的具体算法流程如图1所示，虚线轮廓斜体文字部分是本研究的重要环节。

图  1  频率域NAD全波形反演流程图

Figure  1.  Flow chart of NAD full waveform inversion in frequency-domain

下载: 全尺寸图片 幻灯片

在本文的讨论中，反演目标函数的建立采用不同频率与震源设置下所对应的观测波场与合成波场差值的欧氏范数来表示，记作

式中，$ { （ \cdot ） }^{\mathrm{T}}{\mathrm{和} （ \cdot ） }^{\mathrm{*}} $分别表示向量的转置与复共轭，$ \boldsymbol{p}\in {R}^{2m {\text{×}} 1} $为反演模型的参数向量，$ {{{\boldsymbol{\delta}} }\mathit{d}}^{ （ {n}_{{\mathrm{s}}} ） } （ \boldsymbol{p} \text{，}\omega ） $表示在角频率$ \omega $取值下第$ {n}_{{\mathrm{s}}} $个震源的合成波场与观测波场的残差向量。假设在地表有$ {N}_{r} $个接收器，则$ {{{\boldsymbol{\delta}} }\boldsymbol{d}}^{ （ {n}_{{\mathrm{s}}} ） } $的分量表达形式为

式中，${\boldsymbol{w}}_{i}^{ （ {n}_{{\mathrm{s}}} ） } （ \boldsymbol{p} \text{，} \omega ） \text{，} {\boldsymbol{d}}_{i}^{ （ {n}_{{\mathrm{s}}} ） } （ \boldsymbol{p} \text{，} \omega ） $表示由第$ {n}_{{\mathrm{s}}} $个震源激发并且在第$ i $个接收器上的合成波场和观测波场数据，$ {\boldsymbol{\delta }\boldsymbol{d}}_{i}^{ （ {n}_{{\mathrm{s}}} ） } （ \boldsymbol{p}\text{，}\omega ） $为相应的残差数据分量。

梯度类优化迭代格式的一般形式为

式中，$ k $表示迭代步指标，$ {\mathrm{\alpha }}_{k} $为搜索步长，${d_k} ＝ - {\nabla _{ {\boldsymbol{p}}}}E （ {{\boldsymbol{p}}_k} ） ＋ {\beta _{k － 1}}{d_{k － 1}}$表示第$ k $步的共轭梯度搜索方向，参数${\beta _k} ＝ {{{\nabla _{ \boldsymbol{p}}}E （ {{\boldsymbol{p}}_{k ＋1}} ） \cdot {\nabla _{ \boldsymbol{p}}}E （ {{\boldsymbol{p}}_{k ＋ 1}} ） }}{ {/}} ［ {{{\nabla _{ \boldsymbol{p}}}E （ {{\boldsymbol{p}}_k} ） \cdot {\nabla _{ \boldsymbol{p}}}E （ {{\boldsymbol{p}}_k} ） }} ］ $。

梯度迭代算法的关键环节在于计算反演目标函数的梯度，接下来详细列出这一过程，将式（14）两端分别关于模型参数$ \boldsymbol{p} $进行求导可得（Pratt et al，1998）

其中${\mathrm{Re}} （ \cdot ） $表示取实部，$ \boldsymbol{J} $ 表示雅克比矩阵，即

在实际计算中，显式地求解雅克比矩阵需要大量的正演过程，所消耗计算代价过大，应当尽可能避免。具体的化归思路可叙述为：对线性方程组式（12）两端同时关于模型参数$ \boldsymbol{p} $进行求导，则有

将式（19）代入式（18）得

记$\boldsymbol{y}＝{\boldsymbol{C}}^{-\mathrm{T}}{\boldsymbol{\delta }\boldsymbol{d}}^{*}$，同时将式（20）代入梯度表达式式（17）得

观察式（21）可知，梯度值的计算可避免直接求解雅克比矩阵，而只需求解一次关于${\boldsymbol{C}}^{\mathrm{T}}$的线性方程组获取虚拟震源激发的反传波场（Lang et al，2020）。$ {v}_{{\mathrm{P}}} $和$ {v}_{{\mathrm{S}}} $与拉梅系数$ \lambda $和$ \mu $有着如下一一对应关系

接下来需求解阻抗矩阵的转置对模型参数的偏导数，首先将式（22）中P，S波速度和拉梅系数的关系代入式（3），可得

然后分别对$ {\boldsymbol{P}}_{1} $，$ {\boldsymbol{P}}_{2} $，$ {\boldsymbol{P}}_{3} $关于$ {v}_{{\mathrm{P}}}，{v}_{{\mathrm{S}}} $求偏导，结果如下

式中$ {\boldsymbol{Q}}_{1} $，$ {\boldsymbol{Q}}_{2} $，$ {\boldsymbol{Q}}_{3} $，$ {\boldsymbol{Z}}_{1} $，$ {\boldsymbol{Z}}_{2} $，$ {\boldsymbol{Z}}_{3} $均为对称矩阵，则阻抗矩阵式（13）对$ {v}_{{\mathrm{P}}} $求导得如下稀疏矩阵

同理，阻抗矩阵$\boldsymbol{C}$对$ {v}_{{\mathrm{S}}} $求导得到与式（25）稀疏结构完全相同，非零子块不同的稀疏矩阵，其子块表示为${\boldsymbol{F}}_{i}^{ （ j ） } （ i＝1 \text{，}\cdots 9\text{，}j＝1\text{，}\cdots \text{，}N ） \mathrm{。}{\boldsymbol{F}}_{i}^{ （ j ） }$具体表达式见附件。

为叙述方便，本文将阻抗矩阵式（13）改写为如下形式：

根据模型参数向量排列规则，若将${\boldsymbol{p}} $写作分量形式则为${\boldsymbol{p}} ＝ （ {{p_1}}\text{，} {{p_2}}\text{，} \cdots \text{，} {{p_{2m}}} ） _{2m {\text{×}}1}^{{\mathrm{T}}}$。为缓解P波与S波之间的串扰，缩减梯度计算规模，本文将P波与S波速度进行分离，令${\boldsymbol{p}} ＝ （ {{q_1}}\text{，} {{q_2}} \text{，} \cdots \text{，} {{q_{2i}}} ） ^{{\mathrm{T}}}$，P波速度出现在模型参数分量$2i－1$位置处的取值，而S波速度为模型参数分量$ 2i $处的取值，其中$i＝\mathrm{1\text{，} 2} \text{，} \cdots \text{，} m$。

将阻抗矩阵对模型参数分量的偏导数分解为两部分讨论，一部分为阻抗矩阵对$ {v}_{\mathrm{P}} $分量的偏导数，另一部分为对$ {v}_{\mathrm{S}} $分量的偏导数也可相应推导，则有

记$\boldsymbol{w} ＝ \left( {\begin{array}{*{20}{c}} {{\boldsymbol{w}_1}} \\ {{\boldsymbol{w}_2}} \\ \vdots \\ {{\boldsymbol{w}_N}} \end{array}} \right) \text{，} {\text{ }}{\boldsymbol{y}} ＝\left( {\begin{array}{*{20}{c}} {{{\boldsymbol{y}} _1}} \\ {{{\boldsymbol{y}} _2}} \\ \vdots \\ {{{\boldsymbol{y}} _N}} \end{array}} \right)$，$ {\boldsymbol{w}}_{i} $，$ {{\boldsymbol{y}}}_{i}\in {\boldsymbol{C}}^{6 {\text{×}} 1} $，$i＝\mathrm{1 \text{，} 2} \text{，} \cdots N$，则

为了便于计算，根据式（27）—（30）可将阻抗矩阵对$ {v}_{\mathrm{P}} $的偏导数式（25）改写为如下形式

同理，阻抗矩阵对$ {v}_{\mathrm{S}} $的偏导数改写为

按照式（13）$ \boldsymbol{C} $的稀疏结构，子块矩阵${\boldsymbol{C}}_{ij}＝ {\text{∂}} {\boldsymbol{C}}_{ij}/{ {\text{∂}} {q}_{2i-1}}$和${\boldsymbol{C}}_{ij}＝{ {\text{∂}} {\boldsymbol{C}}_{ij}}/{ {\text{∂}} {q}_{2i}} （ i＝\mathrm{1\text{，} 2} \text{，} \cdots \text{，} N ） $大部分为零矩阵，非零取值为式（13）中的子块${\boldsymbol{C}}_{i}^{ （ j ） }$分别对$ {v}_{\mathrm{P}} $和$ {v}_{\mathrm{S}} $的偏导数.事实上，这两个矩阵是通过阻抗矩阵$ \boldsymbol{C} $分别对$ {v}_{\mathrm{P}} $与$ {v}_{\mathrm{S}} $求偏导数获得，由于矩阵大量元素值为零，实际计算中所需的计算代价均非常小。通过观察式（21）的表达式，进一步计算过程为

所以

代入到${\nabla _{ {\boldsymbol{p}}}}{\boldsymbol{E}} ＝\left(\cdots \text{，} {{\boldsymbol{w}^{\text{T}}}\dfrac{{ {\text{∂}} {\boldsymbol{C}^{\text{T}}}}}{{ {\text{∂}} {p_{2i － 1}}}}{{\boldsymbol{y}}} } \text{，} {{\boldsymbol{w}^{\text{T}}}\dfrac{{ {\text{∂}} {\boldsymbol{C}^{\text{T}}}}}{{ {\text{∂}} {p_{2i}}}}{{\boldsymbol{y}}} }\text{，} \cdots\right)^{{\mathrm{T}}}$中，把式（35）中的两式交叉组合即可求得梯度向量。

经过上述过程求出反演目标函数梯度值后，根据共轭梯度法原理，搜索方向也可随之确定。在反演迭代过程中，另一个需要确定的关键参数是先搜索步长$ {\alpha }_{k} $。为保证良好的迭代效果且不至于消耗过多的正演计算次数，可采用二次抛物线拟合策略（Tape et al，2007；Lang et al，2020）。

3.   数值实验

为验证所研究频率域弹性波全波形反演方法的可行性，本节分别选取异常块模型、双层介质模型和Marmousi模型进行数值实验。反演迭代过程采取逐频反演策略（Bunks et al，1995）：上一个频率的迭代结果作为下一频率的初始模型，即先利用低频数据勾勒出模型的大致结构，再利用高频反演细化内部结构。震源项统一选取雷克子波施加在垂直方向，其频率域表达式为

式中，$ A $为振幅，$ {f}_{0} $为主频，$ f $为频率。在本文的所有数值实验中，振幅和主频统一设置为1.0×10⁵和10 Hz，PML层数均设置为10层。三次实验的部分参数设置如表1所示。

表  1  三种介质模型的反演参数设置

Table  1.  Inversion parameter settings of three media models

	正演网格点数（nx×nz）	空间步长h/km	计算区域/km2	介质密度ρ/（kg·m−3）
异常块模型	101×81	0.025	1.75×2.25	2.63
双层介质模型	101×81	0.025	1.75×2.25	2.63
Marmousi模型	233×78	0.015	5.775×1.9	2.63

下载: 导出CSV 

| 显示表格

3.1   异常块模型

反演的真实模型如图2所示：图2a为P波速度模型，背景速度为3 km/s，中间是边长为0.4 km的矩形异常块，其$ {v}_{\mathrm{P}} $取值为3.5 km/s；图2b为S波速度模型，背景速度为2 km/s，异常块区域的速度为2.5 km/s。初始模型设置为真实模型的背景速度，即$ {v}_{\mathrm{P}}＝ $3 km/s，$ {v}_{\mathrm{S}}＝ $2 km/s.

图  2  P波速度v_P （a）和S波速度v_S （b）异常块反演真实模型

Figure  2.  Real model of anomaly block inversion for P-wave velocity v_P （a） and S-wave velocity v_S （b）

下载: 全尺寸图片 幻灯片

接收器和震源的放置方式如图3所示。根据上文所介绍的频率域全波形反演原理可知，增加接收器数目不会引起计算代价的增加，而震源的数目过多会增大计算量（尽管这一因素可通过震源间并行计算得到缓解）。因此可选择适当的震源数量来保证反演成像的准确性和清晰度。在整个反演阶段，接收器（图中蓝点）在初始模型PML边界逐个排列，上下分别安置98个，共计196个；震源（红色稀疏点）在接收器靠内一层排列，共计2×27＝54个。

图  3  震源和接收器位置图

Figure  3.  Location maps of seismic source and receiver

下载: 全尺寸图片 幻灯片

异常块反演共采用20个频率，频率取值区间为1—30 Hz，间隔为1.5 Hz。第1—10个频率最大迭代步数为120步，剩余频率的最大迭代步数为100步。低频时$ f $＝10.5 Hz的反演成像如图4a所示，可以看出反演模型的大致轮廓开始显现并接近真实模型，但异常块和背景区域的界面还较模糊，此时反演误差较大。为进一步提高成像效果，迭代过程继续向高频阶段进行，当迭代至$ f $＝30 Hz时，相应的反演成像如图4b所示，相比于图4a，异常块和背景的界限逐步趋于清晰，对应的反演结果基本接近真实模型，成像效果较好。

图  4  P波速度v_P （左）和S波速度v_S （右）异常块模型第一阶段（a）及第二阶段（b）反演结果图

Figure  4.  The first-stage （a） and the second-stage （b） inversion results of the anomaly block model for P-wave velocity v_P （left） and S-wave velocity v_S （right）

下载: 全尺寸图片 幻灯片

当挑选某一频率考察反演误差下降情况时，频率取值$ f $＝9 Hz所对应的反演误差曲线如图5所示：相对误差随着迭代步数的增加总体呈下降趋势，并逐步接近于零；在迭代次数小于20次时，相对误差下降较快，然后逐步平稳并最终下降至0.1 （一个量级）以内，这与大多数成功反演情形所表现出的趋势类似（张文生等，2015）。

图  5  反演误差曲线

Figure  5.  Inversion error curve

下载: 全尺寸图片 幻灯片

众所周知，当异常块与背景速度差异较大时，反演的难度也随之上升。为了进一步测试NAD方法弹性波场模拟的有效性，本文将P波速度调整至4 km/s，S波速度调整至3 km/s，背景速度和震源接收器的放置位置均不变。此时的成像结果如图6所示，可以观察到异常块清晰可见，整体反演效果明显。

图  6  P波速度v_P （a）和S波速度v_S （b）异常块模型反演结果图

Figure  6.  The inversion results of the anomaly block model for P-wave velocity v_P （a） and S-wave velocity v_S （b）

下载: 全尺寸图片 幻灯片

3.2   双层介质模型

双层介质反演的真实模型如图7所示：图7a为P波速度的双层介质模型，其中上层$ {v}_{\mathrm{P}} $为3 km/s，下层$ {v}_{\mathrm{P}} $为3.4 km/s；图7b为S波速度模型，其中上层$ {v}_{\mathrm{S}} $为1.732 km/s，下层$ {v}_{\mathrm{S}} $为1.963 km/s，反演初始模型均选取真实模型的上层速度值。

图  7  P波速度v_P （a）和S波速度v_S （b）双层介质真实模型

Figure  7.  Real model of double-layer medium for P-wave velocity v_P （a） and S-wave velocity v_S （b）

下载: 全尺寸图片 幻灯片

接收器和震源采用井间数据（cross-hole）方式放置（Wang，Rao，2006），如图8所示，震源（红色离散点）在反演区域的左右两侧间隔放置，共计40个，接收器（蓝色点）左右放置，不间隔，共计240个。

图  8  震源和接收器位置图

Figure  8.  Location maps of seismic source and receiver

下载: 全尺寸图片 幻灯片

双层介质反演共采用30个频率，频率取值区间为1—30 Hz，间隔为1 Hz。每个频率最大迭代步数为130步。第一阶段迭代至$ f $＝10 Hz的反演模型如图9a所示，可以看出整体的双层结构较为准确，但上下层的分界面比较模糊，分辨率较低。当迭代至30 Hz时所对应的反演结果如图9b所示，相较于图9a，上下层的分界面趋于清晰，图9b对模型的细节刻画更加准确，反演结果与真实模型基本相似，精度较高。

图  9  P波速度v_P （左）和S波速度v_S （右）双层介质模型第一阶段（a）和第二阶段（b）反演结果图

Figure  9.  The first-stage （a） and the second stage （b） inversion results of the two-layer medium model for P-wave velocity v_P （left） and S-wave velocity v_S （right）

下载: 全尺寸图片 幻灯片

同样地，提取频率取值为2 Hz时双层介质模型反演的误差曲线如图10所示（实际上，其它频率迭代情况也基本类似），当迭代步数在40以内时，相对误差下降速度较快，随后相对误差下降速度变慢，逐步趋于平缓，并且接近于0。

图  10  反演误差曲线

Figure  10.  Inversion error curve

下载: 全尺寸图片 幻灯片

接下来探索只浅部接收的方式生成合成数据时的反演情况，震源接收器的放置如图11所示，接收器放置于右侧，共计59个，震源放置于左侧，共计30个，其它参数设置均不变。最终成像效果如图12所示，可以看到S波上下层界面依旧清晰，P波界面较为模糊。因为本文反演实验主要用到直达波震相，所以其整体反演效果略逊色于图9b深浅两侧接收数据时的成像。

图  11  震源和接收器位置图

Figure  11.  Location maps of seismic source and receiver

下载: 全尺寸图片 幻灯片

图  12  P波速度v_P （a）和S波速度v_S （b）双层介质模型反演结果图

Figure  12.  The inversion results of the two-layer medium model for P-wave velocity v_P （a） and S-wave velocity v_S （b）

下载: 全尺寸图片 幻灯片

3.3   Marmousi模型

Marmousi真实模型如图13a所示，P波和S波满足一定的比例关系，即$ {v}_{\mathrm{P}}＝\sqrt{3}{v}_{\mathrm{S}} $ （齐诚等，2006），S波速度变化范围为1.5—5.5 km/s （Youn，Zhou，2001）。图13b显示的初始模型形态是真实模型经过高斯平滑化得到的，可看到只具有大致的层状结构。

图  13  P波速度v_P （左）和S波速度v_S （右） Marmousi真实模型（a）和初始模型（b）

Figure  13.  Marmousi real model （a） and initial model （b） for P-wave velocity v_P （left） and S-wave velocity v_S （right）

下载: 全尺寸图片 幻灯片

接收器和震源的放置方式如图14所示，为获得高分辨率成像结果，震源（红色离散点）分上下两行且只间隔一个节点放置，各安置71个，共计142个，接收器（蓝色点）逐个放置，共计2×230＝460个。

图  14  震源和接收器位置图

Figure  14.  Location maps of seismic source and receiver

下载: 全尺寸图片 幻灯片

迭代过程选取1—30 Hz范围内的频率值，均分为20个频率，间隔为1.5 Hz。在频率取值1—20 Hz内的最大迭代步数为150步，20—30 Hz内的最大迭代步数为120步。第一阶段反演至$ f $＝15 Hz的反演成像如图15a所示，可以看出整体的结构轮廓与真实模型大体一致，但是中间层的细节结构清晰度较差，不易区分，深层的高速体描绘刻画不够理想。迭代至第二阶段$ f $＝30 Hz的反演成像如图15b所示：与图13b相比，细节刻画更加准确，滑脱背斜结构显现更加清晰，深层高速体呈现效果较好，反演结果与真实模型基本一致。整体迭代过程再一次诠释了这一规律：低频时勾勒出了大体轮廓，高频反演则反映更多细节结构。

图  15  P波速度v_P （左）和S波速度v_S （右） Marmousi模型第一阶段（a）和第二阶段（b）反演结果图

Figure  15.  The first-stage （a） and the second stage （b） inversion results of Marmousi model for P-wave velocity v_P （left） and S-wave velocity v_S （right）

下载: 全尺寸图片 幻灯片

Marmousi模型反演实验中，在频率取值为4.5 Hz时的误差曲线如图16所示，整体的趋势与前两个模型很类似：先快速下降，后平缓接近于零。

图  16  反演误差曲线

Figure  16.  Inversion error curve

下载: 全尺寸图片 幻灯片

为了进一步验证近似解析离散化方法的有效性，本文对Marmousi模型的观测数据加入40%的高斯白噪声（Pyun et al，2009）后进行反演。其它实验参数设置均不变，最终成像结果如图17所示，可以观察到滑脱背斜结构显现清晰，深层高速体呈现效果较好，反演结果比较接近真实模型，分辨率较高。

图  17  P波速度v_P （a）和S波速度v_S （b） Marmousi模型反演结果图

Figure  17.  The inversion results of Marmousi model for P-wave velocity v_P （a） and S-wave velocity v_S （b）

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文在频率域弹性波方程NAD正演模拟的基础上，进一步构造相应的全波形反演方法。首先通过NAD算子离散频率域弹性波方程形成大型稀疏线性方程组；然后充分利用相应阻抗矩阵的稀疏分块结构，研究如何花费较小计算代价得到阻抗矩阵对模型参数向量的导数信息与反演目标函数梯度的计算方法，进而建立基于NAD算子的频率域弹性波反演框架；最后，通过数值实验分别对异常块模型、双层介质模型和Marmousi模型进行反演计算。在这三种实验的绝大多数频率迭代过程中，通过反演误差趋势图可发现相对误差均下降一个量级（10%以内），最终也得到了准确可靠的反演结果，从而验证了所提出频率域弹性波反演方法的有效性。

数值频散是由波动方程数值离散引起的非物理现象。减少数值频散的一种有效方法是加密网格，但此方法会带来计算成本的显著增加。由于NAD方法相较于传统的有限差分方法和交错网格方法可以更加有效地压制数值频散并且具有较高的稳定性，所以NAD的总CPU时间较短，在大规模地震图计算中效率更高，此方法在弹性波反演中具有一定的优势。全波形反演作为一种高精度成像方法，目前还有两个难点需要进一步解决：第一，正演过程中离散弹性波方程会得到大型稀疏线性方程组，如何高效求解该方程组是全波形反演面临的难题；第二，全波形反演流程主要用到了直达波震相，对于反射波的利用相对较少。对于以上两个问题，作者今后将构造更加高效的求解线性方程组的方法，同时将进一步研究利用反射波震相的全波形反演算法。