0 引言

随着计算机硬件水平和地震观测技术的不断发展，越来越多的研究者开始关注基于波动方程的高分辨率地震成像方法(Tape et al., 2009; Lei et al., 2020; Liu et al., 2018a)，例如双程波逆时偏移(刘洪等，2009；Yan et al., 2015)和波动方程伴随层析成像(Tarantola, 1984; Pratt and Worthington, 1990; Tromp et al., 2004).相比于传统射线走时地震成像技术，一方面，基于波动方程的成像方法可同时利用地震观测数据中的走时信息和波形信息，这为地下结构反演提供更多的有力约束(Tarantola, 1986; Tape et al., 2007)；另一方面，波动方程成像方法以声波或者弹性波方程作为数学模型，能够更加贴切地模拟实际地震波在地下传播的动力学特征(Liu et al., 2014b).

基于波动方程的地震成像既可在“时间-空间”域进行，也可在“频率-空间”域开展(Pratt et al., 1998; Chen and Cao, 2018).相比于时间域情形，频率域计算主要具有以下四方面优势：(1)由于频率域波动方程只涉及空间偏导数项并且关于各个频率点是相互解耦的，所以频率域计算具有天然的并行特性；(2)可以避免数值离散误差的累积(Pratt, 1999)；(3)通过引入复速度或者复频率(Song and Williamson, 1995)，能够更方便地刻画衰减效应(Wang and Rao, 2009)；(4)频率域反演成像只需选择少量频率点计算(Sirgue and Pratt, 2004)，消耗的计算代价较小.因此，研究频率域全波形成像方法具有重要的科学意义与实际应用价值(Chen and Cao, 2016).

类似于时间域情形，频率域成像过程的主要计算量也集中于对地震波传播规律的正演模拟.因此，正演算法的优劣直接决定着全波形成像结果的精度和整体计算效率.目前用于地震波场模拟的主流正演算法包括有限差分方法(Alford et al., 1974; Virieux, 1984; Moczo et al., 2007)、有限元方法(Marfurt, 1984; Liu et al., 2014a；Liu et al., 2017c)、伪谱法(Fornberg, 1975; Kosloff and Baysal, 1982; 黄继伟和刘洪，2020)、谱元法(Komatitsch et al., 2005; Liu et al., 2017a)等.有限差分方法因其原理直观、格式构造简单、实现容易、便于并行等优点，已成为地震波模拟和成像的有力工具之一(Yang et al., 2006; 刘少林等，2013).然而传统有限差分方法容易引起数值频散效应，当遇到复杂地质结构(例如强间断面)时，往往需要加密离散网格以提升数值精度，这将导致计算效率降低(Moczo et al., 2000; Liu et al., 2017c).近似解析离散化(Nearly Analytic Discrete，NAD)方法作为一种新型差分方法，由杨顶辉等(Yang et al., 2003)首先引入到时间域地震波场数值计算中.该方法的基本构造思想是同时采用波场值以及波场梯度值来近似高阶偏导数项(Yang et al., 2004)，这不仅继承了普通有限差分方法的优点，同时也能在粗网格下有效地压制数值频散，进而得到准确的数值结果，因此该方法可以提高地震波数值模拟的计算效率(Tong et al., 2013; Liu et al., 2015).

求解离散频率域波动方程形成的大型稀疏线性代数方程组是频率域全波形成像的核心问题之一.该方程组的系数矩阵(也叫做“阻抗矩阵”(Pratt, 1999))通常情况下具有不定性，并且在某些常用地震波频率取值下病态程度严重，这为有效求解相应线性方程组带来了挑战.目前主流的求解算法有两种：直接法和迭代法(Ernst and Gander, 2012).直接法是基于对阻抗矩阵进行带有排序(George, 1977)的三角分解，它的优点是对于多震源情形，分解过程只需进行一次，可以减少计算代价.然而直接法通常只适用于二维小尺度区域地震波模拟与波形成像计算.这是由于该方法不能充分利用阻抗矩阵的稀疏分块结构，导致分解矩阵产生大量非零元素，引起内存和计算量的急剧增加，计算效率偏低.对于大尺度区域模型频率域地震波场模拟，迭代法是较为合适的选择.这类方法主要包括基于代数预处理子的Krylov子空间方法(Saad, 2003)、多重网格(Multigrid，MG)方法(Plessix, 2007)以及区域分解方法(Domain Decomposition Methods，DDM)(Gander et al., 2007)等.迭代法在计算过程中只涉及若干次稀疏矩阵与向量乘积，采用稀疏处理技术可以避免内存和计算操作数的大量消耗(Bai, 2015)，从而更加适用于大规模计算.

为了降低阻抗矩阵的病态程度从而达到快速稳定求解线性方程组的目的，Lang和Ren(2015)提出采用不精确旋转分块三角预处理子(Inexact Rotated Block Triangular Preconditioners, IRBTP)结合广义极小残量(Generalized Minimum Residues, GMRES)迭代方法(Saad and Schultz, 1986).其中“不精确”是指采用不精确预处理的思想，在保证准确性的前提下可充分利用相关矩阵的稀疏特性，节约计算成本.近期的研究成果(Lang et al., 2020)进一步完善了IRBTP预处理矩阵的特征值性质，可以从理论上保证相应预处理迭代方法的收敛特性.通过采用数值试验与其他经典迭代方法进行比较，IRBTP预处理迭代方法在加速频率域声波波场模拟方面的优势得到了充分体现(Lang and Yang, 2017).然而IRBTP预处理方法对于频率域弹性波场模拟的数值效率尚不明确，因此有必要进一步研究.

由于弹性波方程比声波方程更能准确地刻画地震波在地下介质中的传播规律，为了贴近实际情形，通常采用弹性波动方程作为模拟地震波传播的数学模型.然而弹性波方程的结构较为复杂，相应求解过程比声波方程消耗的计算代价更多，因此构造高效的频率域弹性波动方程数值模拟方法至关重要.本文主要介绍频率域弹性波方程四阶NAD方法的详细离散过程，同时也分别采用四阶普通有限差分(Ordinary Finite Difference，OFD)方法(Charl-Hyun et al., 1996)和四阶交错网格(Staggered Grid，SG)(Moczo et al., 2000)方法进行数值离散，在检验IRBTP预处理迭代方法的数值效率之后将其运用于相应的弹性波场计算中.分别通过对各种介质模型进行波场模拟、数值频散分析以及与解析解的波形对比，验证NAD方法较其他两种典型数值方法在提高波场求解准确度和计算效率方面的优势.

1 频率域弹性波方程的NAD离散过程

考虑各向同性二维频率域弹性波动方程：

(1)

其中u, v分别表示弹性波场的水平和垂直分量，ρ表示介质密度，λ和μ是拉梅系数，ω=2πf表示角频率(f是频率)，s₁和s₂分别是震源的水平和垂直分量.若引入矩阵和向量记号

(2)

则(1)式可以表示为矩阵分块形式

(3)

当使用NAD方法进行数值离散时，需要对(3)式分别沿x和z方向求偏导数(Liu et al., 2017b, 2018b)，则可得到微分方程组

(4)

与此同时，为了消除人工边界处产生的反射波，需要采用吸收边界条件.本文选取完美匹配层(Perfect Matched Layer, PML)边界条件(Komatitsch and Tromp, 2003).具体地，首先引入复坐标(以x方向为例)

(5)

其中表示虚单位，表示衰减函数，δ表示所考虑方向吸收层的厚度，w是计算区域和PML区域交界面处的波速，r表示PML区域内的网格点到达计算区域的投影距离(有正负差异)，R是离散后的理论反射系数，通常可取作10^-3.为了表述方便，引入记号，则实坐标x与复坐标有如下偏导数关系：

(6)

z方向也完全类似.将(4)式中的实坐标换成复坐标，再代入(5)和(6)式，则形成带有PML吸收边界条件的微分方程组

值得注意的是，(7)式可看作是(4)式的推广：当网格节点位于PML区域之外时，d_x=d_z=0，则(7)可退化为(4)式.

若令Δx和Δz分别为x和z方向上的空间离散步长，采用四阶NAD方法所对应的网格差分模板(具体推导过程可见附录A)进行离散，则(7a)式可以离散为

(8)

为了简化表达式，再引入记号，则(7b)式离散为

(9)

(7c) 式离散为

(10)

其中I₂表示二阶单位矩阵.

若将二维计算区域中的网格节点按行排列，在每一行中按照从左向右的次序并且同一节点处的波场及其梯度值按照如下规则放置：

(11)

根据(8)-(10)式，可形成线性方程组

(12)

阻抗矩阵C具有如下稀疏分块结构：

(13)

其中B_i^(j)(i=1, …, 9, j=1, …, N)表示6×6非零子块矩阵，未列出部分均为零元素，N=n_x×n_z表示所有的离散网格节点个数，n_x, n_z分别表示x和z方向上的节点数，‘ ’表示两个非零子块之间有(n_x-3)个6×6零子块矩阵.根据(13)式可以看到矩阵C每个分块行的非零子块个数不超过9个，进一步地，C的所有非零元素个数大致为2(67(n_x-1)×(n_z-1)+19(n_x+n_z-4)+7)=Ο(134N).子块矩阵B_i^(j)的具体元素表达式见附录B.

值得注意的是，C是复值矩阵并且元素取值依赖于网格规模N，空间离散步长Δx, Δz，以及频率f等因素.若设置参数N=41×41，Δx=Δz=0.02 km，f=15 Hz，则矩阵C的特征值分布如图 1所示(图中虚线分别表示实轴和虚轴所在位置)，可看到C的特征值位于复平面内虚轴的两侧(实部有正有负)并且某些特征值非常接近于0.从数值代数的观点来看，C是不定矩阵并且接近于奇异矩阵，病态程度较高.对于地震波常用频率取值范围内的其他参数情形，阻抗矩阵C的这些特点依然存在.这就为快速而稳定地求解线性方程组(12)带来了困难，进而影响频率域反演的整体计算效率.因此，构造高效的求解算法对于频率域弹性波场模拟和全波形反演至关重要.

2 不精确旋转分块三角预处理迭代方法 2.1 不精确旋转分块三角预处理子的构造过程

为了快速有效地求解复线性方程组(12)，本文采用一种基于代数预处理子的预处理Krylov子空间迭代方法.首先分别提取方程(12)中矩阵和向量的实部与虚部，则有

(14)

其中C_r, x_r, b_r分别表示C, x, b的实部，C_i, x_i, b_i表示虚部.通过基本代数运算，将等式两端的实部与虚部分别对应，则方程组(14)可以等价地转化为实值线性方程组

(15)

针对(15)式中实矩阵A所具有的特殊分块‘2×2’结构，可以构造预处理子M，使得预处理矩阵AM^-1接近于单位矩阵并且条件数远远小于阻抗矩阵C，从而关于AM^-1的线性方程组求解过程稳定并且需要较小的计算代价，达到加速求解的目的.此外，M^-1在迭代计算过程中不必显式求出，而只需求解以M为系数矩阵的广义残量方程组，这也称作“预处理过程”(Saad, 2003).根据以上分析，预处理子M应当是A的一个好的近似，并且具备某种特殊结构，由此可利用A的子块矩阵C_r和C_i分别构造“不精确旋转分块下三角(Inexact Rotated Block Lower Triangular, IRBLT)”预处理子

(16)

“不精确旋转分块上三角(Inexact Rotated Block Upper Triangular, IRBUT)”预处理子

(17)

以及“不精确旋转分块三角因子(Inexact Rotated Block Triangular Factor, IRBTF)”预处理子

(18)

其中

(19)

表示Givens正交旋转矩阵，I_N是N阶单位矩阵，α>0为常数(通常可以取α=1).(αC_r+C_i)_approx表示子块矩阵αC_r+C_i的某个可逆稀疏近似，可以通过Chebyshev半加速迭代公式(Bai et al., 2013)或者不完全LU分解(ILU)(Saad, 2003)来构造，预处理子命名中“不精确”一词也是由此产生.

对于系数矩阵AM^-1非对称的一般型线性方程组，能够有效求解的Krylov子空间迭代方法主要包括GMRES方法和BiCGSTAB方法(Sleijpen and Fokkema, 1993).相比于后者，GMRES方法只针对(15)式中系数矩阵A本身而不涉及法方程组，计算过程中无需A^T的信息，从而具有更广泛的适用性.此外，在这类三角分块预处理子中，IRBTF预处理子所对应的预处理迭代方法具有最好的收敛特性，获得给定精度解所需的迭代步数最少(Lang et al., 2020).然而IRBTF预处理子结构较为复杂，每一步预处理过程消耗过多的计算时间，从而导致求解所需的总时间高于IRBLT与IRBUT预处理子，实际计算效率偏低；通过数值算例测试，对于频率域波动方程求解问题而言，IRBLT的计算效率略优于IRBUT预处理子.因此，在下文的数值试验中统一采用IRBLT预处理子结合GMRES方法(简记为IRBLT-GMRES)求解所对应的线性方程组.关于这类预处理迭代方法的具体实现步骤以及“不精确”预处理过程的处理方法可参见文献(Lang and Ren, 2015).

2.2 不精确旋转分块三角预处理子的数值评估

我们采用IRBLT-GMRES方法求解线性方程组(15)，再与基于复转换(complex-shifted)预处理子M₀(Laird, 2001)以及不完全LU分解(Gander and Nataf, 2005)的预处理GMRES方法(分别简记为M₀-GMRES和ILU-GMRES)直接求解方程组(14)进行对比，考查各种方法在频率域弹性波场模拟方面的数值效率.

考虑计算区域0≤x, z≤5 km，网格规模n_x×n_z=201×201, PML层数为20层.在此区域内P波速度v_P=5 km·s^-1, S波速度v_S=3.0 km·s^-1.实验过程中，沿水平和垂直方向设置相等的空间离散步长，记作h=Δx=Δz=0.025 km; 同时统一采用标准Ricker子波在垂直方向上施加震源作用力(刘少林等，2014)，它的频率域表达式为(Lang and Yang, 2017)

(20)

其中Amp表示振幅，f₀表示震源主频，可取作20 Hz. 根据以上参数设置，可对阻抗矩阵以及震源右端项取值进行计算.在线性方程组的迭代求解过程中，初始迭代解均设置为0.停止准则是一旦当前迭代步残量的欧氏范数小于或者等于初始残量范数的10^-6，即达到迭代终止条件.

表 1列出了三种迭代方法在不同频率取值下计算弹性波场所需要的时间，其中后两列括号中列出的加速比定义为M₀-GMRES或ILU-GMRES方法所需计算时间与IRBLT-GMRES计算时间的比值.可以看出所有的加速比均大于1，IRBLT-GMRES方法消耗的计算时间最少.相对于M₀-GMRES和ILU-GMRES方法，IRBLT-GMRES方法最多可分别加速3.98倍和15.79倍.此外，通过按列观察M₀-GMRES方法和ILU-GMRES方法的计算时间，可以看到随着频率取值的增加，这两种方法所需计算时间急剧增长，而关于ILU-GMRES方法的计算时间只随着频率增长而缓慢增加，计算过程最为稳定.

3 波场模拟

我们运用四阶NAD方法、四阶OFD方法以及四阶SG方法分别在均匀介质模型、双层介质模型和Marmousi模型中进行波场模拟来对比各种方法的数值效率.由于从频率域波场结果中较难判断所得数值解的精确程度，需要对单频波波场作离散Fourier逆变换(IDFT)合成时间域波场，通过观察数值频散现象来衡量这三种数值方法的准确性.根据对Ricker子波震源的频谱分析(Liu et al., 2018b)，可采用0~2.5f₀这一频段的波场数据进行时间域波场计算.此外，本文还通过与解析解的波形对比和数值频散分析进一步考查NAD方法在压制数值频散方面的优势.

为了定量刻画离散网格的疏密程度，这里定义网格频率(grid frequency, G_f)为每个波长WL中的最小网格点数，可表示为(Liu, 1997)

(21)

在实际地震波的传播过程中，最小波速v_min通常取作S波速度v_S.

3.1 均匀模型

首先进行均匀介质中的波场模拟，基本计算参数设置见表 2，P波速度v_P=5.6 km·s^-1, S波速度v_S=3.0 km·s^-1. 在此参数设置下，可计算出网格频率G_f=3.

图 2显示由四阶NAD方法计算得到的单频波波场快照，其中第一行子图表示水平分量，第二行是垂直分量，每列子图分别对应于频率f=10, 20, 30 Hz. 由此可以合成时间域波场，相应的计算过程也对四阶OFD方法和四阶SG方法同样进行.在t=0.25 s时刻由三种方法计算得到的时间域波场结果如图 3所示，可以看到NAD方法计算的波场快照没有明显数值频散，而对于另外两种方法，则有比较清楚的频散现象，其中OFD方法对应的频散情况最为严重.

为了进一步验证所讨论数值方法的准确性，我们放置接收器于(1.0 km, 2.0 km)处，合成对应于三种方法的正则化解与解析解(王美霞等，2012)的波形对比如图 4所示.可以看出由NAD方法计算的数值解与解析解的波形匹配程度最高，而对于另外两种方法，则存在比较明显的偏差.

若使OFD和SG方法能够压制数值频散，需要加密相应的离散网格.当OFD方法的网格参数设置为h=0.015 km, n_x=n_z=267，SG方法的参数设置为h=0.017 km, n_x=n_z=233，相应的波场快照如图 5所示，其中位于左列的子图对应于OFD方法，右列对应于SG方法，第一行为水平分量，第二行为垂直分量，通过观察未发现明显的数值频散.比较三种方法在成功压制数值频散时所需的计算时间(见表 3)，可以看到NAD方法用时最少，相比于次之的OFD方法可以节省8.9%，比SG方法减少20.4%.

3.2 数值频散分析

我们分析三种方法的数值频散特征，其中频率域四阶NAD方法的频散分析过程可见附录C，而对于另外两种数值方法，可参见文献(Charl-Hyun et al., 1996; Chen and Cao, 2016).由于S波的频散情况更具有代表性(刘少林等，2015)，本文列出了对于不同波传播角度(θ)和纵横波速比(r)取值下的S波数值速度与实际波速的比值v_S^num/v_S随网格频率的倒数1/G_f变化的S波频散曲线，如图 6所示.通过对比可发现在各种参数情形下，NAD方法的数值速度最接近于实际波速，频散误差最小，其次是SG方法.

3.3 双层介质模型

接下来考虑双层介质模型，弹性波速度在该介质计算区域的上下两半部分取值不同.具体地，P波速度，S波速度 ，其他计算参数的设置情况见表 2.

图 7显示了双层介质中由四阶NAD方法计算得到的单频波波场快照，其中第一行子图为水平分量，第二行为垂直分量，每列子图分别对应频率取值为f=15, 25, 35 Hz.由此合成对应于三种方法在t=0.25 s时刻的时间域波场快照如图 8所示.类似于均匀介质情形，四阶NAD方法的波场快照结果最为准确，而另外两种方法则存在明显的数值频散，SG方法的频散情况略好于OFD方法.

同样地，为了使四阶OFD方法和SG方法能够成功地压制数值频散，分别加密离散网格至参数设置为h=0.0147 km, n_x=n_z=272以及h=0.0169 km, n_x=n_z=237, 相应的波场快照见图 9，其中左右两列子图分别对应于OFD和SG方法，第一行子图为水平分量，第二行为垂直分量，通过观察未看到明显的数值频散.比较三种方法在得到准确波场结果的前提下所花费的最少计算时间(如表 3所列)，可以看出NAD方法用时最少，相比于OFD方法和SG方法，可以分别加速9.1%和19.6%.

3.4 Marmousi模型

最后考察四阶NAD方法在复杂介质模型模拟中的数值特性，我们选取Marmousi模型(见图 10)(Versteeg, 1994), 其S波速度变化范围是1.5~5.5 km·s^-1, P波速度设置为·v_S, 震源主频下调至f₀=10 Hz，其他参数设置见表 2.

根据四阶NAD方法计算得到的单频波波场快照如图 11所示，其中左列子图为水平分量，右列为垂直分量，每一行分别对应频率取值为6 Hz、12 Hz、18 Hz. 若在地表层每个网格节点处放置接收器，所得理论地震图如图 12所示.此外，在t=0.4 s、0.8 s、1.2 s时刻的时间域波场快照见图 13.通过观察图 12和13，可以看到四阶NAD方法成功压制数值频散，波场结果符合实际情形.

4 结论

在NAD方法成功运用于频率域声波方程数值模拟以及全波形反演之后，本文进一步研究了这种方法在频率域弹性波场模拟方面的数值效率.对于各向同性介质中的频率域弹性波动方程，我们详细推导了四阶NAD方法的离散过程并得到大型线性代数方程组.针对阻抗矩阵稀疏分块结构与数学性质的深入分析结果，揭示了对该线性方程组进行快速有效求解的本质困难.为此，本文引入IRBLT-GMRES预处理迭代求解算法，并将其运用于弹性波场计算.相比于经典的M₀-GMRES和ILU-GMRES预处理迭代方法，IRBLT-GMRES方法最高可分别加速频率域弹性波场模拟大致3.98倍和15.79倍.接着我们采用四阶NAD方法分别在均匀介质和双层介质中进行波场模拟、数值频散分析以及与解析解的波形对比，并与四阶OFD和四阶SG方法进行比较.各种数值结果均显示了NAD方法在压制数值频散和提高计算效率方面的优势，相比于OFD方法和SG方法，可分别加速约9%和20%.此外，四阶NAD方法还用于Marmousi模型的弹性波场模拟并得出无明显数值频散的地表观测记录与波场快照，这进一步说明了NAD方法在复杂介质模型数值模拟中依然可靠.由此可见，IRBLT-GMRES迭代方法配合NAD方法可作为频率域弹性波数值模拟与全波形反演的有力工具.

附录A 四阶NAD方法网格差分模板的具体构造过程

以x方向为例，若对波场u在(i, j)点处的二阶和三阶偏导数项进行离散，NAD格式所对应的网格差分模板可表示为

(A1)

为了计算(A1)式中的模板系数，需要将u_{i-2, j}, u_{i-1, j}, u_{i+1, j}, 以及u_{i+2, j}在(i, j)处沿x方向作五阶Taylor公式展开，并且将和作四阶展开，则可形成线性方程组

(A2)

和

(A3)

求解(A2)和(A3)式，可得网格差分模板系数

(A4)

将其代入(A1)式，即可得到四阶NAD方法的离散格式，z方向也可类似推导.

附录B 阻抗矩阵的非零子块元素表达式

根据(8)-(10)式，公式(13)中非零子块矩阵B_i^(j)的具体元素构成可以表示为

(B1)

(B2)

(B3)

(B4)

(B5)

(B6)

(B7)

(B8)

(B9)

其中j为分块行号，对应于计算区域内网格节点的不同位置，上述表达式中d_x和d_z的取值也随之变化.

附录C 关于频率域弹性波方程四阶NAD方法的频散分析

类似于声波情形(Lang and Yang, 2017)，考虑齐次频率域弹性波动方程及其偏导形式

(C1)

采用四阶NAD方法所对应的网格差分模板对(C1)式进行数值离散，再代入平面谐波解U=U₀e^{-i(k_x·x+k_z·z)}，其中U₀=[u₀, v₀]^T表示常向量，k_x, k_z分别是波数k在x和z方向上的分量，则有

(C2)

矩阵

(C3)

的元素可以表示为

(C4)

(C5)

(C6)

(C7)

(C8)

(C9)

(C10)

(C11)

(C12)

其中E_x⁺=e^ik_xh, E_x^-=e^-ik_xh, E_z⁺=e^ik_zh, E_z^-=e^-ik_zh.若考虑S波频散情形，需要在(C2)式的两边同时除以ρv_S²，则有

(C13)

其中矩阵，再利用波速与拉梅系数的关系式和，以及

(C14)

其中r=v_P/v_S表示纵横波速比.将(C4)-(C12)式中的子块矩阵P₁, P₂, P₃分别替换为Q₁, Q₂, Q₃，即可得到矩阵A_h的表达式.根据(C13)式知，为A_h的特征值，记作λ(A_h)，由此便可得到S波数值速度v_S^num和理论波速vS的比值表达式

(C15)

事实上，对于P波的频散分析完全类似，所得结果只相差一个因子r².