地球物理学报  2021, Vol. 64 Issue (11): 4150-4165   PDF    
三维弹性波方程的修正时空优化保辛数值求解方法
王健1, 贺茜君2, 董兴朋1, 杨顶辉1, 李静爽3, 黄雪源2, 周艳杰2     
1. 清华大学数学科学系, 北京 100084;
2. 北京工商大学数学与统计学院应用统计系, 北京 100048;
3. 中国矿业大学(北京)理学院, 北京 100083
摘要:正演计算是反演研究的基础,为了实现基于三维弹性波方程的全波形反演成像,发展准确、高效、低数值频散的三维正演模拟方法至关重要.为此,本文将修正保辛分部龙格-库塔格式与优化有限差分算子结合,发展了用于数值求解三维弹性波方程的修正时空优化保辛方法(MTSOS).新方法使用二级龙格-库塔格式达到了三阶时间精度,且更适用于求解非均匀介质情况下的弹性波方程,数值频散误差小于同精度保辛分部龙格-库塔(SPRK)方法的误差,提高了计算精度.波场模拟结果表明,三维MTSOS方法可以精确给出数值模拟结果,能够清晰模拟地震波传播过程中产生的各种震相、有效压制数值频散.
关键词: 数值模拟      三维弹性波方程      非均匀介质     
A modified time-space optimized symplectic method for solving 3D elastic wave equations
WANG Jian1, HE XiJun2, DONG XingPeng1, YANG DingHui1, LI JingShuang3, HUANG XueYuan2, ZHOU YanJie2     
1. Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China;
2. School of Mathematics and Statistics, Beijing Technology and Business University(BTBU), Beijing 100048, China;
3. School of Science, China University of Mining and Technology(Beijing), Beijing 100083, China
Abstract: Forward modeling is the basis of inversion research. In order to achieve full waveform inversion imaging based on 3D elastic wave equation, it is necessary to develop accurate and efficient 3D forward modeling methods with low numerical dispersion. For that reason, this paper combines the optimized finite difference operator with the modified symplectic partitioned Runge-Kutta (SPRK) scheme and extends it to 3D to develop a modified time-space optimized symplectic method (MTSOS) for solving the elastic wave equation in 3D inhomogeneous media. The method uses a second-order scheme to achieve third-order time accuracy, and is more suitable for solving the inhomogeneous medium model. The numerical dispersion error is smaller than the SPRK method, which improves the calculation accuracy. The numerical results show that the MTSOS method can give numerical simulation results accurately. The various phases in the inhomogeneous medium can be clearly seen, and there is no visible numerical dispersion, which shows the effectiveness of the new method.
Keywords: Numerical simulation    Three-dimensional elastic wave equation    Inhomogeneous medium    
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)

其中uRunvRvn分别表示广义动量和广义坐标,H为系统的哈密尔顿函数.如果有

(2)

fg分别是关于uv的线性函数,则称系统(1) 为线性可分的哈密尔顿系统.Qin和Zhang(1990)指出,波动方程是一个线性可分的哈密尔顿系统.对于弹性波方程:

(3)

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

(4)

其中

(5)

(6)

(7)

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

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

(8)

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

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

(9)

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

经计算,修正SPRK格式的系数向量为:

(10)

1.2 空间格式

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

(11)

其中Δx表示离散网格的空间步长.可利用不同的格式系数bj(m)发展不同的有限差分离散方法.常见的离散方法利用u(x0+jΔx)在x=x0处的Taylor级数展开,保证格式具有事先给定的精度,通过求解线性方程组确定格式系数bj(m).但这类方法在粗网格下可能会产生强烈的数值频散(Fei and Larner, 1995).如果利用减小空间步长的方法来压制数值频散,又会增加计算时间.为解决这一问题,一系列的优化有限差分方法被提出(Yang et al., 2003, 2006, 2010; Ma et al., 2011, 2014, 2018; Zhang and Yao, 2013).其中,Zhang和Yao(2013)将数值波数与真实波数的无穷范数误差作为目标函数.对于函数的一阶空间偏导数,定义误差函数为:

(12)

其中,kc是给定误差限后的最大真实波数.通过模拟退火法(simulated annealing algorithm, Kirkpatrick et al., 1983; Sen and Stoffa, 2013),得到最优系数bn(1).对于二阶空间偏导数,类似地定义目标函数为:

(13)

通过模拟退火法,得到最优系数bn(2).这样,就得到了式(11)所需的格式系数.Zhang和Yao(2013)证明了,使用如上定义的优化的格式系数所形成的有限差分格式,可以减少经典高阶显式有限差分格式的数值频散.

但在非均匀介质中离散空间算子(5)时,需要考虑介质的速度或弹性参数,对于弹性波方程,还要考虑混合偏导数.因此,不能直接利用式(11)对空间算子L进行离散.为此,Ma等(2018)提出了一种用于二维非均匀介质中波动方程空间算子离散的优化有限差分格式.对于声波方程,以为例,其离散格式定义为(Ma et al., 2018):

(14)

其中,bi0(2)取值参见Zhang和Yao(2013)中的Table 2.

对于弹性波方程,需要处理单向二阶空间偏导数和混合二阶空间偏导数.对于单向二阶空间偏导数,以为例,定义为(Ma et al., 2018):

(15)

其中,b(2)i0的取值参见Zhang和Yao(2013)中的Table 2.对于混合二阶空间偏导数,以为例,定义为(Ma et al., 2018):

(16)

其中,b(1)i0b(1)j0取值参见Zhang和Yao(2013)中的Table 1.对于其他的二阶空间偏导数,可以类似得到它们的离散格式,进而得到空间算子L的离散格式,即离散非均匀介质波动方程中高阶空间偏导数的优化有限差分格式,具体可参见附录.

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

在利用有限差分方法离散空间算子时,如果研究区域内存在速度间断,特别是当不连续面与离散网格不匹配时,会影响波动方程数值模拟的准确性(Moczo et al., 2002).为解决这一问题,Moczo等(2002)提出了等效介质参数法,对于密度ρ采用算术平均,对于弹性参数λμ和(λ+2μ)采用调和平均处理,以提高数值模拟的准确性.具体来说,对于点(xi, yj, zk)处的密度值ρi, j, k,其等效介质参数值ρi, j, kA为以其为中心的立方体内密度值的算术平均,即(Moczo et al., 2002):

(17)

而对于点(xi, yj, zk)处的弹性参数λi, j, k, μi, j, k和(λ+2μ)i, j, k,它们的等效介质参数值λi, j, kH, μi, j, kH和(λ+2μ)i, j, kH分别为以它们为中心的立方体内弹性参数值的调和平均,即(Moczo et al., 2002):

(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)

用该格式得到虚拟层的波场后,即可利用传统的方法离散空间算子中的各项高阶空间偏导数.需要注意对于k=1处关于z方向的混合偏导数,可利用偏心格式进行离散.具体来说,以为例,其数值离散格式分别为(Nilsson et al., 2007):

(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表示时间步长,hxyz表示空间步长.可以得到:

(24)

其中增长矩阵G为:

(25)

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

(26)

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

(27)

解上述不等式,可以得到:

(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

对于非均匀情形,可采用冻结系数法(陆金甫和关治,2004)的思路,将库朗数定义为,其中cmax表示研究区域中的最大波速.而对于弹性波情形,可以将库朗数定义为,其中vp表示P波的最大波速.

2.2 数值频散分析

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

(31)

其中,cnum为数值波速,λgrid为数值波长,ωnum为数值角频率,为空间采样率,且由Nyquist-Shannon采样定理可知有,Sp≤0.5时格式才有意义(Dablain, 1986).

将谐波解式(23)代入式(24),可以得到:

(32)

由于方程必有非零解,因此有:

(33)

即:

(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时的最大数值频散误差 Table 2 The maximum numerical dispersion error of MTSOS method and SPRK method with different spatial accuracies when α=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的点震源,震源函数为雷克子波分量分别在震源所在的XYXZYZ平面的波场快照.从图中可以清晰的发现直达P波及直达S波,且没有可见的数值频散.

(37)

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

图 3中给出了在t=3.0 s时刻位移场的不同

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

图 5 六阶MTSOS方法计算得到的t=1.75 s时刻位移场的不同分量分别在震源所在的XYXZYZ平面的波场快照 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的含水裂隙,裂隙的两端恰好分别位于震源所在的XYYZXZ平面内.在含水裂隙内部,密度为1000 kg·m-3,P波速度为1.5 km·s-1,S波速度为0 km·s-1.采用六阶空间精度的MTSOS方法进行波场模拟.在t=0.25 s时刻位移场不同分量在震源所在的XYXZYZ平面的波场快照如图 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时刻位移场的不同分量分别在震源所在的XYXZYZ平面的波场快照 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新方法与三维全波形反演方法相结合,得到基于三维弹性波方程的新的全波形反演算法,并应用于实际地震数据的地震成像研究中,以获得地球内部高分辨率的成像结果.

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

以下给出式(5)中各项的2N阶优化有限差分格式离散形式.其中bi(1)bi(2)(i=1, …, N)的取值分别参见Zhang和Yao(2013)中的Table 1Table 2.

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

References
Alterman Z, Karal F C Jr. 1969. Propagation of elastic waves in layered media by finite difference methods. Bulletin of the Seismological Society of America, 59(1): 471. DOI:10.1785/BSSA0590010471
Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-200. DOI:10.1006/jcph.1994.1159
Clayton R, Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations. Bulletin of the Seismological Society of America, 67(6): 1529-1540. DOI:10.1785/BSSA0670061529
Dablain M A. 1986. The application of high-order differencing to the scalar wave equation. Geophysics, 51(1): 54-66. DOI:10.1190/1.1442040
Fei T, Larner K. 1995. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport. Geophysics, 60(6): 1830-1842. DOI:10.1190/1.1443915
Feng K, Qin M Z. 1987. The symplectic methods for the computation of Hamiltonian equations. //Zhu Y I, Guo B eds. Numerical Methods for Partial Differential Equations. Berlin: Springer, 1-37.
Feng K, Qin M Z. 2003. Symplectic Geometric Algorithms for Hamiltonian Systems (in Chinese). Hangzhou: Zhejiang Science and Technology Press.
Haras Z, Ta'asan S. 1994. Finite difference schemes for long-time integration. Journal of Computational Physics, 114(2): 265-279. DOI:10.1006/jcph.1994.1165
He X J, Yang D H, Ma X, et al. 2020. A modified numerical-flux-based discontinuous Galerkin method for 2D wave propagations in isotropic and anisotropic media. Geophysics, 85(5): T257-T273. DOI:10.1190/geo2019-0485.1
Kim J W, Lee D J. 1996. Optimized compact finite difference schemes with maximum resolution. AIAA Journal, 34(5): 887-893. DOI:10.2514/3.13164
Kirkpatrick S, Gelatt C D, Vecchi M P. 1983. Optimization by simulated annealing. Science, 220(4598): 671-680. DOI:10.1126/science.220.4598.671
Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method. Geophysics, 47(10): 1402-1412. DOI:10.1190/1.1441288
Lan H Q, Zhang Z, Xu T, et al. 2012. Effects due to the anisotropic stretching of the surface-fitting grid on the traveltime computation for irregular surface by the coordinate transforming method. Chinese Journal of Geophysics (in Chinese), 55(10): 3355-3369. DOI:10.6038/j.issn.0001-5733.2012.10.018
Lee C, Seo Y. 2002. A new compact spectral scheme for turbulence simulations. Journal of Computational Physics, 183(2): 438-469. DOI:10.1006/jcph.2002.7201
Liu S L, Yang D H, Ma J. 2017. A modified symplectic PRK scheme for seismic wave modeling. Computers & Geosciences, 99: 28-36.
Lu J F, Guan Z. 2004. Numerical Solution of Partial Differential Equations (in Chinese). 2nd ed. Beijing: Tsinghua University Press.
Lysmer J, Drake L A. 1972. A finite element method for seismology. Methods in Computational Physics: Advances in Research and Applications, 11: 181-216.
Ma J, Yang D H, Tong P, et al. 2018. TSOS and TSOS-FK hybrid methods for modelling the propagation of seismic waves. Geophysical Journal International, 214(3): 1665-1682. DOI:10.1093/gji/ggy215
Ma X, Yang D H, Zhang J H. 2010. Symplectic partitioned Runge-Kutta method for solving the acoustic wave equation. Chinese Journal of Geophysics (in Chinese), 53(8): 1993-2003. DOI:10.3969/j.issn.0001-5733.2010.08.026
Ma X, Yang D H, Liu F Q. 2011. A nearly analytic symplectically partitioned Runge-Kutta method for 2-D seismic wave equations. Geophysical Journal International, 187(1): 480-496. DOI:10.1111/j.1365-246X.2011.05158.x
Ma X, Yang D H, Song G J, et al. 2014. A low-dispersive symplectic partitioned Runge-Kutta method for solving seismic-wave equations: Ⅰ. scheme and theoretical analysis. Bulletin of the Seismological Society of America, 104(5): 2206-2225. DOI:10.1785/0120120210
Ma X, Yang D H. 2017. A phase-preserving and low-dispersive symplectic partitioned Runge-Kutta method for solving seismic wave equations. Geophysical Journal International, 209(3): 1534-1557. DOI:10.1093/gji/ggx097
Martin R, Komatitsch D, Gedney S D, et al. 2010. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML). Computer Modeling in Engineering and Sciences, 56(1): 17.
Moczo P, Kristek J, Vavryčuk V, et al. 2002. 3D heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities. Bulletin of the Seismological Society of America, 92(8): 3042-3066. DOI:10.1785/0120010167
Nilsson S, Petersson N A, Sjögreen B, et al. 2007. Stable difference approximations for the elastic wave equation in second order formulation. SIAM Journal on Numerical Analysis, 45(5): 1902-1936. DOI:10.1137/060663520
Okunbor D, Skeel R D. 1992. Explicit canonical methods for Hamiltonian systems. Mathematics of Computation, 59(200): 439-455. DOI:10.1090/S0025-5718-1992-1136225-3
Priolo E, Seriani G. 1991. A numerical investigation of Chebyshev spectral element method for acoustic wave propagation. //13th World Congress Computation and Applied Mathematics, 2: 551-556.
Qin M Z, Zhang M Q. 1990. Multi-stage symplectic schemes of two kinds of Hamiltonian systems for wave equations. Computers & Mathematics with Applications, 19(10): 51-62.
Sanz-Serna J M. 1988. Runge-Kutta schemes for Hamiltonian systems. BIT Numerical Mathematics, 28(4): 877-883. DOI:10.1007/BF01954907
Schmid R. 1987. Infinite dimensional Hamiltonian systems. //Monographs and Textbooks in Physical Sciences. Biblionopolis.
Sen M K, Stoffa P L. 2013. Global Optimization Methods in Geophysical Inversion. Cambridge: Cambridge University Press.
Smith W D. 1975. The application of finite element analysis to body wave propagation problems. Geophysical Journal International, 42(2): 747-768.
Wang J, Yang D H, Jing H, et al. 2019. Full waveform inversion based on the ensemble Kalman filter method using uniform sampling without replacement. Science Bulletin, 64(5): 321-330. DOI:10.1016/j.scib.2019.01.021
Wang L, Yang D H, Deng X Y. 2009. A WNAD method for seismic stress-field modeling in heterogeneous media. Chinese Journal of Geophysics (in Chinese), 52(6): 1526-1535. DOI:10.3969/j.issn.0001-5733.2009.06.014
Yang D H, Teng J W, Zhang Z J, et al. 2003. A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media. Bulletin of the Seismological Society of America, 93(2): 882-890. DOI:10.1785/0120020125
Yang D H, Peng J M, Lu M, et al. 2006. Optimal nearly analytic discrete approximation to the scalar wave equation. Bulletin of the Seismological Society of America, 96(3): 1114-1130. DOI:10.1785/0120050080
Yang D H, Song G J, Hua B L, et al. 2010. Simulation of acoustic wavefields in heterogeneous media: A robust method for automatic suppression of numerical dispersion. Geophysics, 75(3): T99-T110. DOI:10.1190/1.3428483
Zhang J H, Yao Z X. 2013. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm. Journal of Computational Physics, 250: 511-526. DOI:10.1016/j.jcp.2013.04.029
冯康, 秦孟兆. 2003. 哈密尔顿系统的辛几何算法. 杭州: 浙江科技出版社.
兰海强, 张智, 徐涛, 等. 2012. 贴体网格各向异性对坐标变换法求解起伏地表下地震初至波走时的影响. 地球物理学报, 55(10): 3355-3369. DOI:10.6038/j.issn.0001-5733.2012.10.018
陆金甫, 关治. 2004. 偏微分方程数值解法. 2版. 北京: 清华大学出版社.
马啸, 杨顶辉, 张锦华. 2010. 求解声波方程的辛可分Runge-Kutta方法. 地球物理学报, 53(8): 1993-2003. DOI:10.3969/j.issn.0001-5733.2010.08.026
王磊, 杨顶辉, 邓小英. 2009. 非均匀介质中地震波应力场的WNAD方法及其数值模拟. 地球物理学报, 52(6): 1526-1535. DOI:10.3969/j.issn.0001-5733.2009.06.014