三维弹性波方程的修正时空优化保辛数值求解方法

王健， 贺茜君， 董兴朋， 杨顶辉*， 李静爽， 黄雪源， 周艳杰

1 清华大学数学科学系， 北京 100084 2 北京工商大学数学与统计学院应用统计系， 北京 100048 3 中国矿业大学(北京)理学院， 北京 100083

0 引言

地震波场正演模拟是指利用数值方法求解地震波在已知地下介质中传播的技术，它是基于波动方程反问题研究的基础.此外，正演模拟对于解释地震信息，评估观测系统等也具有十分重要的意义.地震波正演模拟有很多实现方法，包括有限差分方法(finitie difference, FD, Dablain, 1986)、有限元法(finite element, FE, Lysmer and Drake, 1972; Smith, 1975; He et al., 2020)、伪谱法(pseudo-spectral, PS, Kosloff and Baysal, 1982)、谱元法(spectral-element, SEM, Priolo and Seriani, 1991)等.这些方法各有优缺点及适用范围，其中，有限差分方法是目前应用最为广泛的正演模拟算法.

有限差分方法将计算区域划分为有限个规则的或不规则的网格点，利用泰勒展开的方式将波动方程中的微分算子转化为网格点上函数值的差分形式，通过求解网格点上的函数值近似得到原波动方程的解.不同的差分形式对应着不同的有限差分方法，其具有易于实现、存储量小、运算速度快、易于并行等优点.Alterman 和Karal(1969)最早将该方法引入到地震波场的正演模拟中，并在后续得到了广泛应用(如, Moczo et al., 2002).然而，有限差分方法面临着可能产生严重数值频散的问题(Fei and Larner, 1995).这是由于利用差分形式近似微分算子时，不同频率波的相速度产生了差别，特别是当网格步长很大或波场梯度很大时，这种现象尤为明显.因此当波传播时间较长时，不同频率的波由于相速度不同，会发生分离.这显然与实际物理定律相违背，故而产生了虚假信息，影响了正演模拟的精度.

为解决这一问题，Yang 等(2003)将近似解析离散化(nearly analytic discrete, NAD)算子引入波动方程的空间偏导数离散中，同时利用网格点上的位移及梯度的差分来近似微分算子.由于利用了网格点上的梯度值，因此与相同数值精度的传统有限差分方法相比，NAD算子长度更小、算子更为紧致.同时，数值实验表明，在相同网格步长的情况下，NAD算子的数值频散明显低于传统的有限差分方法，能够在粗网格情形下有效压制数值频散(Yang et al., 2003).近年来，基于NAD算子的思想，一系列改进算法被提出，包括优化近似解析离散化方法(optimal nearly analytic discrete, ONAD, Yang et al., 2006)、带权重的近似解析离散化方法(weighted nearly analytic discrete, WNAD, 王磊等, 2009)、近似解析保辛分部龙格-库塔方法(nearly analytic symplectic partitioned Runge-Kutta, NSPRK, 马啸等, 2010; Ma et al., 2011)、改善近似解析离散方法(refined nearly analytic discrete, RNAD, Yang et al., 2010)等，并已被应用于全波形成像(Wang et al., 2019).

压制数值频散的另一种思路是提高离散格式相速度的准确性.为达此目的，可以构造优化有限差分格式，将数值波数和真实波数的误差作为优化函数，从而得到最优的差分离散系数(Haras and Ta'asan, 1994; Kim and Lee, 1996; Lee and Seo, 2002).其中，Zhang和Yao(2013)将数值波数和真实波数的无穷范数误差作为目标函数，通过模拟退火法(Kirkpatrick et al., 1983; Sen and Stoffa, 2013)优化目标函数，最终得到最优的差分离散系数.数值结果表明，使用优化的高阶有限差分方法可以大大节省内存需求和计算代价.

波动方程可以被视为一个无限维线性可分的哈密尔顿系统(Hamilton system)(Schmid, 1987)，其相空间存在着辛几何结构，可以表示物体的运动规律.在哈密尔顿体系下构造的辛算法(symplectic method)可以保持哈密尔顿系统的辛几何结构，具有重要意义.Feng和 Qin(1987)根据生成函数理论，给出了任意阶精度的隐式保辛数值格式的构造方法，并据此提出了多种不同类型的辛格式，包括保辛龙格-库塔(Runge-Kutta, RK)格式、蛙跳格式等(冯康等, 2003).Okunbor和Skeel(1992)给出了对于可分哈密尔顿系统的分部龙格-库塔(patitioned Runge-Kutta, PRK)格式的构造方法.Qin和Zhang(1990)构造了一至四阶精度、一至四级的PRK格式.Ma 等(2011)将保辛格式与NAD算子相结合，给出了NSPRK方法，其数值频散误差更小.Liu 等(2017)提出了一种修正的保辛分部龙格-库塔(symplectic partitioned Runge-Kutta, SPRK)格式，使一个二级格式达到了三阶时间精度.Ma和Yang(2017)提出了优化SPRK格式，具有保相位的优点.Ma等(2018)提出了时空优化保辛(time-space optimized symplectic method, TSOS)方法，将优化SPRK格式与优化有限差分格式结合，同时具有保相位和低数值频散的优点，更易于处理非均匀介质.

本文将修正的SPRK方法和优化有限差分方法思想结合，发展了用于求解三维非均匀介质中弹性波方程的修正时空优化保辛方法(modified time-space optimized symplectic method, MTSOS)，并分析了该方法的稳定性及数值频散特性，数值实验表明该方法能有效压制数值频散，同时可以准确、高效地模拟弹性波在地下介质中的传播.

1 修正时空优化保辛(MTSOS)方法

1.1 时间格式

一般的2n维哈密尔顿系统可表示为：

(1)

(2)

且f和g分别是关于u和v的线性函数，则称系统(1)为线性可分的哈密尔顿系统.Qin和Zhang(1990)指出，波动方程是一个线性可分的哈密尔顿系统.对于弹性波方程：

(3)

可将其表示为线性可分的哈密尔顿系统的形式：

(4)

其中

(5)

(6)

(7)

空间算子L的其他项类似可得.

由Sanz-Serna(1988),s阶显式SPRK格式可表示为：

(8)

其中，un和vn表示该系统在第n个时间步的数值解，c=(c1,c2,…,cs)和d=(d1,d2,…,ds)是该格式的系数向量，Δt表示该格式的时间步长.

为保证该格式具有s阶精度，需通过泰勒展开的方式确定系数向量c和d.但要想使格式达到s阶精度，至少需要s级格式才能满足要求.Liu 等(2017)提出了一种修正SPRK格式，具体形式为：

(9)

因为格式(9)中引入了修正项Lv1，可以通过泰勒展开的方式确定系数向量c和d，使一个二级格式具有三阶时间精度.

经计算，修正SPRK格式的系数向量为：

(10)

1.2 空间格式

为了利用修正SPRK格式得到波动方程的数值解，还需要离散式(4)中的空间算子L.有限差分法方法易于编程及并行，是一种常见的数值离散方法.一般的有限差分格式对函数的空间高阶偏导数利用该点及周围点的函数值进行离散，即

(11)

(12)

(13)

(14)

(15)

(16)

1.3 等效介质参数与边界条件

(17)

(18)

可以看出，在均匀介质中，密度和弹性参数的等效介质参数值即为其真实值.在非均匀介质中，只需在网格附近利用算术或调和平均即可得到等效介质参数值.Ma 等(2018)通过数值算例证明了等效介质参数处理可以更好地模拟间断面处地震波的传播.

自由地表边界条件是指地表处的应力为0(兰海强等, 2012).为了模拟自由地表边界条件，Nilsson等(2007)提出了边界修正法，该方法是一种显式方法，具有二阶空间精度，并具有很好的数值稳定性.本文以三维弹性波方程为例说明该方法.

假设z=0处为自由地表，满足自由地表边界条件.为下文表述方便，记(u1,u2,u3)T=(u,v,w)T.在三维弹性波方程中，自由地表边界条件可表示为：

(19)

设z=0对应k=1，边界修正法首先引入一个虚拟层k=0，通过数值离散自由地表边界条件(19)来得到虚拟层的波场.具体来说，用如下格式来离散自由地表边界条件(Nilsson et al., 2007)：

(20)

(21)

(22)

对于其他的k=1处关于z方向的混合偏导数，可以类似地得到其数值离散格式.Nilsson 等(2007)指出，上述的边界修正法是一种二阶方法，且当格式满足CFL条件(Courant-Friedrichs-Lewy condition)时是稳定的.

需要注意的是，由于边界修正法是一种二阶方法，所以当计算区域内部的空间精度超过二阶时，会造成边界与计算区域内精度不匹配的现象，影响计算区域内的空间精度.为解决这一问题，可以在计算区域内从边界层开始逐层提高空间精度，直至达到计算区域内的空间精度，以保证数值格式的稳定性.

由于计算时间和规模的限制，一般只能模拟地震波在给定区域内的传播情况.在这种情况下，需要在计算区域的人工边界添加吸收边界条件(absorbing boundary conditions)，来消除人工边界产生的反射波对计算区域的影响.近年来很多吸收边界条件方法被提出，其中包括旁轴近似方法(Clayton and Engquist, 1977)，完全匹配层吸收边界条件(Berenger, 1994)等.本文采用Martin等(2010)提出的辅助微分方程-完全匹配层(auxiliary differential equation-perfectly matched layer, ADE-PML)吸收边界条件来处理人工边界问题.

将1.1节提出的修正SPRK时间格式与1.2节提出的优化有限差分格式相结合，将空间格式拓展到三维，同时利用等效介质参数，边界修正法和完全匹配层吸收边界条件，即获得了三维非均匀介质弹性波方程正演模拟的修正时空优化保辛(MTSOS)方法.该方法是一种保辛方法，且由于使用了优化有限差分格式和等效介质参数，更适用于非均匀介质中弹性波方程的求解.在第2节，将对三维MTSOS方法的稳定性和数值频散进行分析.

2 三维MTSOS方法的稳定性和数值频散分析

2.1 稳定性条件

本节首先利用傅里叶分析的方法，求出时间离散算子的谱半径，进而得到MTSOS方法的稳定性条件.对于三维均匀介质中的声波方程，将谐波解

(23)

代入声波方程，其中ωnum表示数值角频率，k=(kx,ky,kz)表示波矢量，Δt表示时间步长，h=Δx=Δy=Δz表示空间步长.可以得到：

(24)

其中增长矩阵G为：

(25)

I表示单位算子，L表示空间离散算子.设矩阵G的特征值为ψ，算子L的特征值为.由于波动方程是双曲型偏微分方程，故其空间离散算子<0.由式(25)，ψ的特征多项式f(ψ)和满足关系：

f(ψ)=ψ2-tr(G)ψ+det(G)=ψ2-(2+Δt2

(26)

为保证格式稳定，需要满足|ψ|≤1，即：

(27)

解上述不等式，可以得到：

-12≤Δt2≤0.

(28)

(29)

式(29)对算子L的任意特征值均成立，因此要保证格式稳定，最大库朗数αmax需满足：

(30)

式(30)即为三维MTSOS方法的稳定性条件，可以根据不同精度的空间算子L的最小特征值得到与之对应的最大库朗数.在表1中给出了三维情形下不同精度的MTSOS及TSOS方法的最大库朗数.可以看出，MTSOS方法的稳定性相比TSOS方法有了明显的提升.

表1 三维情形下不同空间精度的MTSOS方法与TSOS方法的最大库朗数Table 1 Maximum Courant numbers of MTSOS method and TSOS method with different spatial accuracies in 3D case

2.2 数值频散分析

为了对格式进行数值频散分析，首先定义数值频散比率为如下形式(Yang et al., 2006)：

(31)

将谐波解式(23)代入式(24)，可以得到：

(32)

由于方程必有非零解，因此有：

det(exp(iωnumΔt)I-G)=0,

(33)

即：

ωnumΔt=arccos(Re(ψ)).

(34)

又由式(26)可知：

(35)

因此，可得：

(36)

式(36)被称为数值频散关系式，它表示了数值频散比率R与库朗数、采样率及空间算子的特征值之间的关系.R越接近于1，说明数值频散越小.接下来，本小节将比较三维情形下MTSOS方法与SPRK方法的数值频散比率，以说明MTSOS方法的优势.

在图1和图2中分别给出了在α=0.1,Sp=0.4时不同空间精度的MTSOS方法和SPRK方法在不同方向上的频散误差|1-R|.从图中可以看出，MTSOS方法与SPRK方法的数值频散误差都表现出了各向异性的特点：数值频散误差的最大值都出现在沿坐标轴的方向，而在与坐标轴呈45°角的方向，数值频散误差最小.同时，随着空间精度的提升，数值频散误差也呈现出下降的趋势.

图1 不同空间精度的MTSOS方法在α=0.1,Sp=0.4时不同方向的数值频散误差(a) 六阶； (b) 八阶； (c) 十阶； (d) 十二阶.Fig.1 The numerical dispersion error of the MTSOS method with different spatial accuracies in different directions when α=0.1,Sp=0.4(a) Sixth order; (b) Eighth order; (c) Tenth order; (d) Twelfth order.

图2 不同空间精度的SPRK方法在α=0.1,Sp=0.4时不同方向的数值频散误差(a) 六阶； (b) 八阶； (c) 十阶； (d) 十二阶.Fig.2 The numerical dispersion error of the SPRK method with different spatial accuracies in different directions when α=0.1,Sp=0.4(a) Sixth order; (b) Eighth order; (c) Tenth order; (d) Twelfth order.

从图中还可以看出，从总体上看，MTSOS方法的数值频散误差要小于SPRK方法.为了量化比较两种方法的数值频散误差，在表2中分别给出了不同空间精度的MTSOS方法与SPRK方法沿各个方向的最大数值频散误差.可以看出，不同空间精度的MTSOS方法的最大频散误差均小于相同空间精度的SPRK方法的最大频散误差，而且随着空间精度的提高，MTSOS方法对于最大频散误差的改进更加明显.

表2 不同空间精度的MTSOS方法及SPRK方法在α=0.1,Sp=0.4时的最大数值频散误差

3 数值算例

本节将MTSOS方法与边界修正法、吸收边界条件等结合，来计算三个三维弹性波模型，以验证MTSOS方法的有效性，以及边界修正法和吸收边界条件对计算区域边界的处理效果.

3.1 均匀介质模型

首先，选取三维均匀介质模型，利用MTSOS方法来进行波场模拟.计算区域为40 km×40 km×40 km，介质密度为2660 kg·m-3，P波速度为5.8 km·s-1，S波速度为3.199 km·s-1.震源位于(20 km,20 km,18 km)处.震源取为主频2 Hz的点震源，震源函数为雷克子波

×exp[-8(0.6f0t-1)2],

(37)

震源位于区域中心.时间步长为2.5 ms，空间步长为0.2 km.计算区域的上方为自由地表边界条件，其他五个方向采用吸收边界条件，吸收层取为15层.采用六阶空间精度的MTSOS方法进行波场模拟.

在图3中给出了在t=3.0 s时刻位移场的不同分量分别在震源所在的XY、XZ和YZ平面的波场快照.从图中可以清晰的发现直达P波及直达S波，且没有可见的数值频散.

图3 MTSOS方法计算得到的t=3.0 s时刻位移场的不同分量分别在震源所在的XY、XZ和YZ平面的波场快照Fig.3 Snapshots of different components of the displacement field at the time t=3.0 s calculated by the MTSOS method in the XY, XZ and YZ planes where the seismic source is located

为了验证MTSOS方法的准确性，在(23 km,24 km,18 km)处设置一个接收器，并将接收器处的波形记录与解析解进行对比.图4是在0～5 s内位移的三分量与解析解的比较图，其中实线表示解析解，虚线表示数值解.从图中可以看出，接收器处的波形记录与解析解能够很好的吻合，说明数值模拟结果是可信的.

图4 六阶MTSOS方法计算得到的位于(23 km,24 km,18 km)处的接收器的波形图(a) x分量； (b) y分量； (c) z分量.Fig.4 Waveforms of the receiver located at (23 km,24 km,18 km) calculated by the sixth-order MTSOS method(a) x component; (b) y component; (c) z component.

3.2 双层介质模型

选取一个三维双层介质模型，来检验MTSOS方法在利用等效介质参数法处理模型间断以及吸收边界条件的有效性.计算区域为20 km×20 km×20 km，上层介质密度为2600 kg·m-3，P波速度为5.8 km·s-1，S波速度为3.2 km·s-1；下层介质密度为3723 kg·m-3，P波速度为9.134 km·s-1，S波速度为4.932 km·s-1，间断面位于z=20 km处.震源取为主频8 Hz的点震源，震源函数为雷克子波.震源位于(10 km,10 km,8 km)处.时间步长为1 ms，空间步长为0.08 km.计算区域的上方为自由地表边界条件，其他五个方向采用吸收边界条件，吸收层取为15层.采用六阶空间精度的MTSOS方法进行波场模拟.

在图5中给出了在t=1.75 s时刻位移场的不同分量分别在震源所在的XY、XZ和YZ平面的波场快照.从图中可以清晰的发现直达波、反射波及透射波的各个震相以及产生的转换波.在XY平面，由于不存在间断面，波场呈同心圆结构，不同的圆表示了以不同速度传播的波场.其中，最外层为直达P波(P)，向内依次为反射P波(PP)、S波反射后转换成的P波(SP)、直达S波(S)、P波反射后转换成的S波(PS)、反射S波(SS)等.在XZ及YZ平面，可以看到经间断面得到的反射、透射和转换波，而且没有可见的数值频散.图6中选取了t=1.75 s时位移场的x分量在XZ平面的波场快照，并标注出所有产生的震相.

图5 六阶MTSOS方法计算得到的t=1.75 s时刻位移场的不同分量分别在震源所在的XY、XZ和YZ平面的波场快照Fig.5 The wave field snapshots of different components of the displacement field at the time t=1.75 s calculated by the sixth-order MTSOS method in the XY, XZ, and YZ planes where the seismic source is located

图6 六阶MTSOS方法计算得到的t=1.75 s时刻位移场的x分量在震源所在的XZ平面的波场快照Fig.6 The wave field snapshots of the x component of the displacement field at the time t=1.75 s calculated by the sixth-order MTSOS method in the XZ plane where the seismic source is located

在图7中以位移场的x分量在XZ平面为例给出了从t=0.5 s至t=2.5 s时刻的波场快照.从图中可以看出波场随时间稳定的传播：最初时刻仅有直达P波及直达S波.当波传播到间断面，产生了各种反射波、透射波及转换波.各种震相均稳定且无可见的数值频散.约t=1.5 s时，透射P波(Pp)传播到边界；约t=2.0 s时，其他波的震相也陆续传播到边界.其中，自由地表处的直达P波经过自由地表的反射形成反射P波和转换S波，而在其他的边界处，波场均被吸收边界条件很好地吸收，未产生可见的反射波.该数值实验很好地说明等效介质参数可以准确地处理模型的间断，从而更适用于求解非均匀介质弹性波方程，而边界修正法与吸收边界条件可以很好地与数值算法结合，给出准确的边界条件.

图7 六阶MTSOS方法计算得到的位移场的x分量在XZ平面上从t=0.5 s至t=2.5 s时刻的波场快照Fig.7 The wave field snapshots of the x component of the displacement field on the XZ plane from t=0.5 s to t=2.5 s calculated by the sixth-order MTSOS method

3.3 含水裂隙介质模型

选取一个三维含水裂隙介质模型来考察MTSOS方法对于细微的非均匀结构的分辨能力.含水裂隙介质是在石油勘探领域经常会接触到的地质结构.这是由于石油公司会向井中注水，用以驱油并减少开采时间.因此，模拟弹性波在含水裂隙介质模型中的传播对于石油勘探具有重要的意义.

在这个例子中，将计算区域取为2 km×2 km×2 km，背景介质密度为2170 kg·m-3，P波速度为2.618 km·s-1，S波速度为1.421 km·s-1.震源取为主频40 Hz的点震源，震源函数为雷克子波.震源位于(1 km,1 km,1 km)处.时间步长为0.2 ms，空间步长为5 m.计算区域的上方为自由地表边界条件，其他五个方向采用吸收边界条件，吸收层取为15层.在(0.9 km,0.9 km,1 km)至(1 km,1 km,1.1 km)有一个宽度为5 m的含水裂隙，裂隙的两端恰好分别位于震源所在的XY和YZ、XZ平面内.在含水裂隙内部，密度为1000 kg·m-3，P波速度为1.5 km·s-1，S波速度为0 km·s-1.采用六阶空间精度的MTSOS方法进行波场模拟.在t=0.25 s时刻位移场不同分量在震源所在的XY、XZ和YZ平面的波场快照如图8所示.在图8中，可以清晰地发现由于含水裂隙产生的散射波.为了更准确地显示散射波的形成过程，分别在(0.925 km,0.92 km,1 km)和(0.885 km,0.88 km,1 km)各放置一个接收器.注意到第一个接收器和震源位于含水裂隙的同侧，而第二个接收器和震源分别位于含水裂隙的两侧.图9和图10分别给出了波形图，其中实线为不含含水裂隙时的波形图，虚线为含含水裂隙时的波形图.右侧为波形图的局部放大图.可以看出，对于位于震源同侧的接收器，由于含水裂隙的存在，在直达S波后可清晰地发现散射波所形成的一系列尾波.而对于位于震源异侧的接收器，散射P波与直达S波几乎同时到达，使得直达S波的振幅与没有含水裂隙时的情形有了明显的区别.该数值实验说明MTSOS方法可以准确模拟勘探尺度下模型内的细微结构变化对于波形的影响.

图8 六阶MTSOS方法计算得到的t=0.25 s时刻位移场的不同分量分别在震源所在的XY、XZ及YZ平面的波场快照Fig.8 The wave field snapshots of the different components of the displacement field at the time t=0.25 s calculated by the sixth-order MTSOS method in the XY, XZ, and YZ planes where the seismic source is located

图9 六阶MTSOS方法计算得到的位于(0.925 km,0.92 km,1 km)接收器的波形图(a) x分量； (b) x分量局部放大图； (c) y分量； (d) y分量局部放大图； (e) z分量； (f) z分量局部放大图.Fig.9 Waveforms at the receiver located at (0.925 km,0.92 km,1 km) calculated by the sixth-order MTSOS method(a) x component; (b) Partial enlarged view of x component; (c) y component; (d) Partial enlarged view of y component; (e) z component; (f) Partial enlarged view of z component.

图10 六阶MTSOS方法计算得到的位于(0.885 km,0.88 km,1 km)接收器的波形图(a) x分量； (b) x分量局部放大图； (c) y分量； (d) y分量局部放大图； (e) z分量； (f) z分量局部放大图.Fig.10 Waveforms at the receiver located at (0.885 km,0.88 km,1 km) calculated by the sixth-order MTSOS method(a) x component; (b)Partial enlarged view of x component; (c) y component; (d) Partial enlarged view of y component; (e) z component; (f) Partial enlarged view of z component.

4 结论

本文发展了适用于三维非均匀介质弹性波方程的正演模拟新方法.该方法在时间上源自Liu等(2017)提出的修正SPRK格式，利用二级格式达到了三阶时间精度；在空间上采用Ma 等(2018)提出的优化有限差分格式，并将其推广到三维情况，同时与等效介质参数、吸收边界条件等处理技术相结合，获得了用于求解三维非均匀介质弹性波方程的MTSOS新方法.该方法尤其适用于求解非均匀介质弹性波方程.我们进一步给出了MTSOS方法的稳定性条件，并进行了数值频散分析.理论分析和数值结果均表明，新方法的数值频散误差低于相同空间精度的SPRK方法，而且随着空间精度阶数的提高，这种优势也越来越明显；同时新方法具有较高的数值稳定性，能有效提高波场模拟的计算效率.波场模拟结果进一步验证了MTSOS方法能够有效压制数值频散，清晰地给出弹性波传播过程中产生的各种波震相，可以模拟非均匀结构对于波场的影响，证实了MTSOS新方法的有效性和准确性.对于含水裂隙模型，MTSOS方法能够准确地反映弹性波经过含水裂隙所形成的散射波，接收器处接收到的波形能充分显示含水裂隙的存在对波形的影响，表明MTSOS方法能有效模拟非均匀结构对波场的影响.进一步地，在未来的研究中，可将本研究提出的MTSOS新方法与三维全波形反演方法相结合，得到基于三维弹性波方程的新的全波形反演算法，并应用于实际地震数据的地震成像研究中，以获得地球内部高分辨率的成像结果.

附录 优化有限差分格式二阶空间偏导数离散格式

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)