2018, Vol. 61
2. 长江大学地球物理与石油资源学院, 武汉 430100;
3. 胜利油田勘探开发研究院西部分院, 山东东营 257001
2. College of Geophysics and Petroleum Resources, Yangtze University, Wuhan 430100, China;
3. West Branch of Shengli Oilfield Exploration and Development Research Institute, Shandong Dongying 257001, China
地震波场数值模拟是应用数值方法模拟地震波在地下介质中的传播过程，计算在地面或地下各观测点地震记录的一种数值模拟方法，是地球物理勘探和地震学的重要基础.目前，常用的数值模拟方法主要有射线追踪法和波动方程法，其中波动方程法有伪谱法、有限元法、边界元法、谱元法和有限差分法(Alterman and Karal, 1968；Graves，1996；董良国等，2000a；Saenger et al., 2000, 2004).随着数值模拟技术的发展和生产实践的要求，围绕着提高有限差分计算效率(黄超和董良国，2009；唐佳等，2016)、模拟精度(Liu，2013；李振春等，2016)、算法稳定性(Virieux, 1986；董良国等，2000；杜启振等，2015)、压制数值频散、处理复杂介质(Fletcher et al., 2008；王璐琛等，2016；程玖兵等，2013)和吸收边界条件(Berenger，1994；Pan et al., 2012)等方面，研究人员已经发展了多种方法.特别是在如何压制数值频散和提高计算效率方面，杨顶辉等在这方面做了大量相关工作(Yang et al., 2009, 2012a, b，2014；Yang and Wang, 2010, Ma and Yang, 2017；Tong et al., 2013；Zhou et al., 2015；He et al., 2015；Liu et al., 2017；Jing et al., 2017)，取得了许多有意义研究成果.
提高差分格式数值计算精度最直接的方法就是在差分计算时增加网格节点数，但这也增加了计算量和所需的存储空间.紧致差分方法恰好能够较好地解决这个矛盾，同时紧致差分是一种隐式差分格式，具有较好的稳定性，这些优势也使得它成为目前研究较多的有限差分方法之一.1989年，Dennis和Hudson(1989)针对Navier-Stokes方程提出了空间四阶的紧致格式，能够使用粗网格获得较高的精度.1992年，Lele(1992)构造了求解一阶导数和二阶导数的紧致差分格式.Adams和Shariff(1996)提出了紧致ENO格式，用于求解激波湍流相互作用问题.Chu和Fan(1998)提出了三点组合型紧致差分格式，并将其用于求解对流扩散方程.随后人们针对不同的问题和方程，提出了许多种不同的紧致差分格式，广泛用于空气动力学、流体力学和电磁波方程等方面.
在地震波场数值模拟方面，王书强等(2002)和Zhou和Greenhalgh(2003)先后将Lele(1992)提出的紧致差分格式用于弹性波方程的数值模拟中.Yang等(2003)年给出各向异性情况下的紧致差分方法以及声波和弹性波数值模拟结果，并给出相应的吸收边界条件.Du等(2009a, b)利用紧致交错网格差分方法对横向各向同性介质地震波场进行了数值模拟，Liu和Sen(2009)研究了任意阶空间导数的紧致格式和一阶空间导数的交错紧致格式，分析了它们的精度和计算效率，并将其用于声波和弹性波数值模拟中.Kosloff等(2010)提出可以利用当前点两边各任意多个点计算导数值的一般紧致格式，并基于Remes方法求取差分系数.杨宽德等(2011)研究了Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟方法.此外，Chu和Stoffa(2010)、Liu(2014)还研究了频率域地震波数值模拟中的紧致差分方法.而组合型紧致差分格式在地震波场模拟中的应用，尚未见到文献报道.本文在声波和弹性波方程时间四阶离散格式的基础上，将组合型紧致格式应用到位移场空间偏导数的求取，实现了二维各向同性介质地震波场的数值模拟.
1 组合型紧致有限差分方法原理早期的紧致差分格式是基于Hermite多项式构造而来的，在Hermite多项式中，相邻三个节点的函数值、一阶和二阶导数值关系可以表示为
|
(1) |
其中fi表示节点函数值，f′i和f″i分别表示一阶和二阶导数值.
1992年，Lele(1992)对Hermite公式进行了扩展，构造了求解一阶导数和二阶导数的紧致差分格式，表示为
|
(2) |
其中a, b, c, α, β为差分系数，h表示网格间距.
对上述Lele差分格式进行泰勒公式展开和求解差分系数，可以得到一阶导数五点中心六阶精度差分格式：
|
(3) |
二阶导数五点中心六阶精度差分格式：
|
(4) |
由公式(3)和(4)可以看出，紧致差分(Compact Difference，以下简称CD)格式的特点就是用相邻节点的函数值求解若干节点上的导数值.而常规中心差分(Central Finite Difference，以下简称CFD)仅求解中心节点的导数值.另一方面，也可以认为常规中心差分是紧致差分的特例，即公式(2)中的α=β=0，例如一阶和二阶导数的六阶精度差分格式表示为
|
(5) |
其中b1=0.75, b2=-0.15, b3=0.0167, c0=-2.7222, c1=1.5, c2=-0.15, c3=0.0111.
常规中心差分如果要得到六阶精度的差分格式，需要用到七个节点，而紧致差分只需要五个点，这使得紧致差分格式比常规差分处理边界节点上更为方便.且一般情况下，在相同网格间距时，紧致差分格式比常规差分具有更高的精度，具有更小的数值频散(王书强等，2002).
Chu和Fan(1998)等构造了精度更高的三点六阶组合型紧致差分格式(Combined Compact Difference，以下简称CCD)：
|
(6) |
与公式(3)和(4)对比可以看出，上述组合型紧致差分格式只需要相邻的三个节点就可以同时求得一阶和二阶导数的六阶精度近似值，比常规紧致差分的节点数更少.公式(6)中的一阶导数和二阶导数是耦合的，既可以同时求出也有利于波形的保真性.
2 声波方程组合型紧致有限差分方法分析二维声波方程可以表示为
|
(7) |
式中u(x, z)为地震波位移，v(x, z)为地震波速度.利用截断的泰勒公式表示n+1和n-1时刻的位移场可以得到：
|
(8) |
|
(9) |
两式相加，略去高次项，并代入声波方程，得到位移场时间四阶精度的差分格式：
|
(10) |
公式(10)为位移场的三层显式差分格式，利用它就可以进行地震波场时间层的推进.在公式中含有位移对x和z的二阶和四阶导数，对这些导数将利用组合型紧致差分格式(6)进行求取.
假设模型纵向和横向节点数为m和n，h为空间网格大小.利用公式(6)求公式(10)中的空间偏导数∂2u/∂z2和∂2u/∂x2，矩阵表示为
|
(11) |
其中A和C为公式(6)左端的差分系数矩阵，大小分别为2m×2m和2n×2n；E和F为待求位移场空间一阶和二阶导数矩阵，大小分别为2m×n和m×2n；B和D为公式(6)右端的差分系数矩阵，大小分别为2m×m和n×2n；U为位移场矩阵，大小为m×n.这些矩阵分别表示为
|
(12) |
|
(13) |
|
(14) |
|
(15) |
|
(16) |
|
(17) |
|
(18) |
公式(11)可以用追赶法进行求解.由公式得出位移场ui, j可同时求得它的空间一阶和二阶导数∂ui, j/∂x、∂ui, j/∂z、∂2ui, j/∂x2和∂2ui, j/∂z2，表示为
|
(19) |
差分公式(10)中的空间四阶导数可以在求得二阶导数后，将二阶导数值当作U，再次使用公式(19)进行求取.
2.1 精度分析不论是利用CCD格式还是利用常规的七点六阶CFD或五点六阶CD格式，进行声波方程数值模拟时，它们在时间层推进方式上是相同的，即都是利用差分格式(10)，时间差分为四阶精度.CD、CFD和CCD不同之处在于它们分别利用公式(4)、(5)和(6)计算位移空间高阶导数.虽然在声波方程的差分格式中并没有空间一阶偏导数，但在下文中的第4节进行弹性波场数值模拟中，存在空间的混合偏导数，需要通过一阶导数的迭代求取，所以这里对这三种格式计算的一阶和二阶导数的截断误差进行对比，结果见表 1，α、β、a、b和c为公式(2)中的差分系数.
|
表 1 一阶和二阶导数不同差分方法的截断误差比较 Table 1 Truncation errors in various difference schemes for the first and second derivative calculations |
从表 1可以看出，虽然三种方法都能达到空间六阶精度，但截断误差有较大的差别.CFD和CD差分格式计算一阶偏导数的截断误差约是CCD的40倍和4.4倍，计算二阶偏导数的截断误差约是CCD的36倍和8.5倍，这说明CCD方法的具有更高的差分精度.
2.2 误差分析通过模拟二维平面谐波初值问题并计算模拟误差，定量分析CCD方法的数值模拟精度.二维平面波初值问题可以表示为
|
(20) |
其中v是平面波的波速，θ0是初始时刻波阵面法线方向(即传播方向)与x轴的夹角，f0是平面简谐波的频率.
上述初值问题的解析解为
|
(21) |
在二维波场数值模拟中，假设f0=20 Hz，θ0=π/4，均匀介质的波速为3600 m·s-1，模型长度和深度均为2000 m，纵横向网格长度相同，采样时间为1 s.在不同的空间网格长度和时间步长条件下，计算数值模拟的相对误差.相对误差定义为
|
(22) |
式中ui, jn表示数值解，u(tn, xi, zj)表示解析解.
研究表明(Yang et al., 2014)，优化的近似解析离散化方法(Optimal Nearly-Analytic Discrete，以下简称ONAD)也具有较高的数值模拟精度和计算效率，为了比较，在此计算CCD、CD、ONAD和CFD四种方法在不同模拟参数情况下的相对误差曲线见图 1.从图 1可以看出，随着空间网格长度和时间步长的增加，四种方法的相对误差均会逐渐增加，并且随模拟时间的增加而增加.
|
图 1 四种数值模拟算法的相对误差曲线 (a) CCD; (b) CD; (c) ONAD; (d) CFD. Fig. 1 Relative error-time curves of the numerical simulation using four different methods and models (a) CCD; (b) CD; (c) ONAD; (d) CFD. |
便于比较，将实验结果列于表 2.从表中可以看出，四种情况下，误差相对关系基本一致，例如当空间步长为20 m，时间步长1 ms时，CCD算法的相对误差仅为0.066%，CD次之，为0.089%，CFD相对误差为0.943%，这与精度分析的结论是一致的.需要说明的是，本次实验的ONAD算法相对误差比CFD方法大，这是由于ONAD只有空间四阶精度，而本文的CFD空间精度为六阶.当采用较细网格时，如空间步长10 m，时间步长0.5 ms时，四阶ONAD方法的误差要小于六阶CFD方法，这也说明了ONAD方法在提高模拟精度方面也有积极的意义.用较小的空间步长(10 m)和时间步长(0.5 ms)时，CCD算法的精度显著提高，相对误差仅有0.0008%，并且相对误差增长缓慢，这说明了采用CCD格式的声波模拟结果精度较高，能进行较长时间的地震波场数值模拟.
|
|
表 2 四种方法在不同情况下的最大相对误差(%)比较 Table 2 Relative errors in four different schemes and parameters |
在数值模拟时，如果空间网格长度使用过大，会造成较大的求解误差，并会产生所谓的数值频散现象，所以频散关系分析既是判断数值模拟方法优劣的重要方法，也是确定空间网格大小的重要依据.通过对声波方程的组合型紧致差分格式进行数值频散分析，进一步分析该方法的适用条件.
首先考虑x方向，令简谐波解ui, jn=Aexp×[I(ihkx+jhkz-kv′nΔt)]，且：
|
(23) |
如果数值模拟时不存在误差，则A, B, C应该满足以下关系：
|
(24) |
但在差分的实际计算中，所产生的数值波数(又称修正波数)与解析波数不同.令：
|
(25) |
将公式(23)和(25)代入差分格式(6)，可以得到方程组：
|
(26) |
其中：ϕx=hkx, ϕ′x=hk′x, ϕ″x=hk″x.解上述方程组可得：
|
(27) |
同理z方向可得：
|
(28) |
式中ϕx=ϕcosθ，ϕz=ϕsinθ，且ϕ=kh.θ是波的传播方向与x轴的夹角，k为波数，h为空间步长.
将公式(23)、(27)和(28)代入声波方程的差分格式(10)中，并令x=kv′Δt，α=vΔt/h.其中v′为数值波速，v是实际波速，Δt是时间步长，α称为Courant数，
|
(29) |
通过上述频散关系确定ϕ后可以解得对应的x.定义数值波速与真速度的比值为
|
(30) |
在压制数值频散方面，优化的近似解析离散化方法(ONAD)方法也具有较好的效果，所以选取CFD、CD和ONAD三种方法与之进行比较.在理想情况下，如果不存在数值频散则速度比γ恒等于1.γ偏离1越大，则说明该方法的数值频散越严重，反之则说明该方法能更好地压制数值频散.图 2显示了以上四种方法在不同的α和θ的时候的速度比曲线.
|
图 2 不同差分方法的速度比曲线 (a) θ=0, α=0.25; (b) θ=π/8, α=0.25; (c) θ=π/4, α=0.25; (d) θ=0, α=0.45; (e) θ=π/8, α=0.45; (f) θ=π/4, α=0.45. Fig. 2 Velocity ratio curves of the numerical simulation in different methods |
取ϕ∈[0, π]作为横坐标，它是波数与空间步长的乘积，单位波长内采样点数N=2π/ϕ，所以横坐标也可以看作N由∞逐渐减小至2.从图中的速度比曲线可以看出：①随着空间采样点数的减小，四种方法的频散现象逐步加剧.CCD、CD和ONAD方法的数值频散比CFD方法要小，其频散曲线更趋近于1，这也说明了紧致差分(CCD和CD)和近似解析离散类(ONAD)方法在压制数值频散方面的都具有较大的优势，都适合使用较大的空间步长，进而提高计算效率；②假设数值速度在理论速度的99.9%~100.1%以内表示不存在频散，则CCD、CD、ONAD和CFD对应的ϕ的极小值分别为1.12、0.98、1.06和0.91，所以它们在最小主波长内需要使用的网格点数分别为5.6、6.4、5.9和6.9个，即CCD和ONAD方法对于空间网格长度的要求最低，其次是CD方法，CFD方法最为严格；③CCD方法的空间网格长度过大时，速度比曲线上翘，数值速度大于真速度，这与另外三种方法相反.这一结论也能从图 3c中可以看出，当ϕ=π，即一个周期内取2个样点，蓝色的速度比曲线大部分要大于1.
|
图 3 不同差分方法的方位角-速度比曲线 (a) CFD; (b) CD; (c) CCD; (d) ONAD. Fig. 3 Azimuth-speed ratio curves of the numerical simulation in different methods |
为研究该差分方法的数值各向异性特征，即速度比与传播方向的关系，绘制极坐标速度比曲线，如图 3所示.图中的极角表示地震波传播方向与x轴的夹角，极径表示速度比.四种方法在空间采样(ϕ=π/3)满足频散要求时，都没有各向异性特征.在空间采样不足时(ϕ=7π/5，ϕ=9π/5)均会出现速度各向异性，0°(或90°)与45°方向上的速度差别最大.
利用二维模型来说明CCD方法在压制数值频散方面的优势.模型长度和深度均为4 km，纵波速度3600 m·s-1，震源子波为f(t)=sin(2πf0t)exp× (-π2f0t2/4)，主频20 Hz.根据频散分析结论，CCD方法在最小主波长内只需要使用约5.6个点，所以设置空间网格间距32 m，时间步长1 ms.图 4为四种方法数值模拟的波场快照，可以看出CCD和ONAD方法的结果清晰，没有频散现象，CD方法在0°和90°方向上存在微弱的频散现象，而CFD的频散非常严重，与前面理论分析结果一致.
|
图 4 空间网格32 m，时间步长1 ms，420 ms时刻波场快照 (a) CFD; (b) CD; (c) CCD; (d) ONAD. Fig. 4 Snapshots at 440ms: Δx=32 m, Δt=1 ms |
稳定性条件也是有限差分数值模拟中一个非常重要的问题，是影响差分方法计算效率的重要因素.采用Fourier方法(张文生，2006)对声波方程的组合型紧致差分格式进行稳定性分析.定义：
|
(31) |
其中
|
(32) |
令θ1=kxh，-π/2≤θ1≤π/2，并利用欧拉公式，将上式化简可得：
|
(33) |
解方程可得：
|
(34) |
其中
|
(35) |
同理，定义：
|
(36) |
且令θ2=kzh，-π/2≤θ2≤π/2，代入公式(6)可以得到：
|
(37) |
其中
|
(38) |
公式(10)是一个三层显式差分格式，为了分析其增长矩阵，将其改写为：
|
(39) |
式中gi, jn=∂neIh(kxi+kzj).
将公式(34、35、37和38)代入(∂2u/∂x2)i, jn=ζ1nexp[Ih(kxi+kzj)]和(∂2u/∂z2)i, jn=ζ2nexp[Ih(kxi+kzj)]，并代入声波差分格式(39)中，化简并写成矩阵格式为
|
(40) |
其中
|
(41) |
为增长矩阵，其中α=vmaxΔt/h为该差分格式的Courant数.
为使差分格式(40)满足稳定性条件，需要其增长矩阵的谱半径小于等于1.直接求解该增长矩阵的特征值较为复杂，用二分法求该问题的解(宋国杰，2008)，计算得到的近似稳定性条件为
|
(42) |
采用2.3节所用模型进行数值模拟，采样时间1 s，分析CCD、CFD、CD和ONAD四种方法的计算效率，各方法使用的数组大小、个数和算法时间复杂度等列于表 3.
|
|
表 3 不同方法计算效率比较 Table 3 Comparison of computational efficiency between different methods |
从表 3中可以看出：CCD和ONAD差分方法的网格半径最少，仅需要3个点，这对于边界处理最有利，但ONAD空间精度仅为四阶，要低于CCD的六阶精度.CFD和CD差分方法虽然空间精度能达到六阶，但网格半径也随之增加.CCD和ONAD均为低数值频散的差分方法，满足无数值频散要求的空间网格长度大于CFD和CD方法.但是ONAD方法需要计算位移场的梯度，这就增加了数组的个数和计算时间.同样的声波方程，CCD仅需要12个数组，而ONAD需要25个，在网格长度相近的情况下，CCD方法占用内存最小.从稳定性来看，要达到同样的时间四阶精度，CFD和CD方法的Courant数最大，CCD次之，ONAD最小.考虑空间步长的不同，CCD、CFD和CD能允许较大的时间步长(7 ms).从算法的时间复杂度来看，CCD和CD算法的复杂度要大于CFD和ONAD，其中n表示网格点数.值得注意的是，在模型大小相同的时候，由于网格长度的差异，会造成四种算法的空间网格节点数n不同，也就造成了实际所需运算时间的不同，表中运算时间是在同一台计算机和相同编程环境下模拟得到的.综合模拟精度、空间网格长度、时间步长和占用内存四个因素来看，CCD方法进行声波方程数值模拟时，它具有高精度(时间四阶、空间六阶)、高计算效率(较大Courant数、占用内存少)和低数值频散(粗网格)的优势.
需要说明的是，表 3中四种方法各参数的选取时，均不考虑计算精度，空间网格和时间步长仅满足无数值频散和稳定性的要求.但在实际计算时，为了精度考虑，二者均需要按要求适当减小，这势必会进一步增加计算量和存储量.
3 弹性波方程组合型紧致有限差分方法各向同性完全弹性介质中的二维波动方程可以表示为
|
(43) |
其中，u(x, z)和w(x, z)表示x和z方向的质点位移，VP(x, z)和VS(x, z)表示纵横波波速.同声波方程一样，利用公式(8)和(9)，并代入波动方程，将u对时间的导数转换为u和w对空间的导数，可以得到时间四阶精度差分格式：
|
(44) |
同理可以得到w的时间四阶精度差分格式为
|
(45) |
弹性波差分格式与声波方程不同，公式(44)和(45)中除了位移u的空间二阶导数∂2ui, j/∂x2和∂2ui, j/∂z2外，还有二阶混合偏导数∂2ui, j/∂x∂z和其他四阶偏导数.假设偏导数∂2/∂x2，∂2/∂z2，∂/∂x和∂/∂z的差分算子为A′xx，A′zz，A′x和A′z，则二阶混合偏导数∂2/∂x∂z的差分算子为A′xz=A′xA′z，四阶偏导数∂4/∂x4的差分算子表示为A′x4=A′xxA′xx.即对于二阶混合偏导数，可以在求出∂ui, j/∂x和∂ui, j/∂z之后，再次使用公式(19)，对另一个方向求导，从而求出空间二阶混合偏导数∂2ui, j/∂x∂z.同样，对差分公式中出现的u和w的其他高阶偏导数，也可以按照上述同样的方法进行迭代求取，这里不再赘述.
4 模型试算为了验证方法的数值性能，我们采用组合型紧致方法对公式(43)和(44)表示的二维各向同性完全弹性介质的弹性波方程进行数值模拟.
4.1 均匀介质模型模型长度和深度均为10 km，介质为泊松体，纵波速度3460 m·s-1，横波速度2000 m·s-1，在模型中心处激发10 Hz主频的雷克子波.由于最小速度为2000 m·s-1，按频散分析的结果，每个波长内需采样5.6个点，所以设置网格长度为35 m.同时，由于最大速度为3460 m·s-1，按稳定性条件，取时间步长为8 ms.长时程数值模拟是数值模拟中的一个重要内容(Chen, 2006; Li et al., 2012; Gao et al., 2016).对于大尺度模型和长时程传播，即便是十分微弱的误差也会随着传播时间的增加而不断累积，最终可能导致不容忽视的数值假象.为了分析CCD方法的长时程数值模拟结果，数值模拟时没有使用吸收边界条件，图 5(a—d)显示的是1 s、2 s、5 s和10 s时刻的位移水平分量波场快照，地震波场特征清晰，无可见的数值频散现象.为了比较，对该模型采用相同的参数，使用时间四阶空间六阶精度CFD格式进行数值模拟，图 5(e—h)为相同时刻的波场快照.可以看出，在满足频散要求和稳定性条件的情况下，CCD方法具有计算效率较高、长时程稳定的特点，能适应较长时间的数值模拟，计算结果较为理想.在相同参数下，CFD结果存在很明显的频散现象，并且随时间的增加，频散更加明显，严重地影响了数值模拟的效率和精度.
|
图 5 均匀介质模型弹性波数值模拟波场快照 (a) CCD(1 s); (b) CCD(2 s); (c) CCD(5 s); (d) CCD(10 s); (e) CFD(1 s); (f) CFD(2 s); (g) CFD(5 s); (h) CFD(10 s). Fig. 5 Wave field snapshots using elastic wave equation in homogeneous medium (a—d) are the snapshots at 1, 2, 5, 10 s for the displacement component in X direction using CCD; (e—h) are the snapshots at 1, 2, 5, 10 s for the displacement component in X direction using CFD. |
设置水平层状介质模型，共20层，每层厚度400 m，模型总厚度8 km，模型长度同样是8 km.第一层纵波速度为3500 m·s-1，往下各层纵波速度比上一层增加100 m·s-1.模型为泊松体，即纵横波波速比为1.732，所以模型中最小和最大波速为2020 m·s-1和5400 m·s-1.在地面中间处(X=4 km，Z=0 m)激发10 Hz雷克子波震源，人工边界采用PML边界条件(刘有山等，2012)处理.按前文分析结果，设置纵横向网格长度为36 m，时间步长为5 ms.图 6显示的是地面接收到的地震记录，记录长度8 s，由于直达波能量明显强于反射波，所以图中对地震记录进行了增益处理.
|
图 6 水平层状介质模型弹性波方程数值模拟地面地震记录 (a) X方向位移分量；(b) Z方向位移分量. Fig. 6 Seismic records of the model′s surface using elastic wave equation in horizontal layers media (a) Displacement component in X direction; (b) Displacement component in Z direction. |
多层模型中地面模拟接收到的地震记录非常清晰，没有数值频散和不稳定数值结果，地震记录中的直达波、反射波和转换波显示清楚，说明组合型紧致有限差分算法可以有效地模拟弹性波在多层各向同性介质模型中的传播.为了检验本文算法长时程数值模拟的结果，图 7给出在较粗网格条件下(网格长度36 m，时间步长5 ms)，CCD和CFD方法分别得到的检波点(X=4 km，Z=0 m)处的波形图，并与利用CFD方法在较细的网格(网格长度10 m，时间步长1 ms)情况下得到的近似理论解的参考波形进行对比.
|
图 7 接收点(4 km，0 km)处接收到的反射波水平位移波形图 (a)整道反射波记录对比；(b) 4000~4500 ms局部放大对比；(c) 6000~6500 ms局部放大对比.其中，红色加点实线是粗网格条件下的CCD数值结果，绿色实线是粗网格条件下的CFD数值结果，黑色实线是细网格条件下的近似理论解. Fig. 7 The horizontal displacement of reflected wave received at receiver (4 km, 0 km) (a) Comparison of the whole reflection wave record; (b) Comparison of the local amplification within 4000~4500 ms; (c) Comparison of the local amplification within 6000~6500 ms.The red dot solid line denotes the solution by CCD under coarse grid condition, and the green solid line denotes the solution by CFD under coarse grid condition, the black solid line denotes the approximate theoretical solution by CFD under tiny grid condition. |
图 7a显示三种情况下得到的反射波到达时间一致，从图 7(b，c)的局部放大结果来看，在粗网格和长时程模拟条件下，CCD方法得到的地震记录平滑无数值频散，与之前的分析结果一致，且CCD方法和近似解析解的波形曲线基本重合，这说明了CCD方法的正确性.而CFD方法在相同条件下存在较严重的数值频散，且模拟时间越长，模拟结果与近似解偏差越大.这些结果和对比说明，CCD方法能够长时间计算稳定性，而且也能很好地压制数值频散.
5 结论与建议有限差分法是地震波场数值模拟的重要方法之一，在目前常规有限差分数值模拟中为了获得更高的模拟精度，需要使用更多的节点，这样会增加计算量和所需的存储空间，还会增加人工边界处理的难度.
本文从组合型紧致差分格式出发，建立了时间四阶、空间六阶的声波方程的差分格式，对该差分格式进行了模拟精度、频散曲线和稳定性分析，并与常用的几种格式进行了综合比较，该格式可同时计算空间一阶和二阶导数，且具有使用节点少(3个)，低数值频散(每个波长采样5.6个点)，较高稳定性(Courant数为0.792)和模拟精度(空间6阶)的特点，适用于大尺度和较长时程地震波场数值模拟.文中还建立了各向同性介质弹性波方程的差分格式，并进行了数值模拟，模拟结果显示地震波场特征清晰，说明了该方法的适用性.该方法不仅可以用于弹性介质的波场模拟，还可以进一步推广到二维或三维的各向异性介质、黏弹介质和双相介质等复杂介质的数值模拟中，为今后研究地震波传播规律、逆时偏移、全波形反演等工作提供了一种有效的方法.组合型紧致差分格式仅需使用三个节点，就能使一阶和二阶空间偏导数的差分精度达到传统差分方法的六阶精度，因而能够大幅度节约内存需求和计算量，对于三维大尺度模型的正演和反演等问题均具有重要意义.
文中在对二阶混合偏导数的求取过程中，先求得位移x方向的一阶导数，然后对该一阶导数再求z方向的一阶偏导数，这样做可能会降低二阶混合偏导数的精度.在今后的研究中，可以考虑建立一种新的组合型紧致差分格式，一次性地求取二阶混合偏导数，以提高数值模拟中的精度.此外，可以通过有限差分系数优化的方法(Liu, 2013; Zhang and Yao, 2013a, b)进一步提升本文方法的紧致性.
致谢 感谢审稿专家对文章内容和格式提出的宝贵建议.
Adams N A, Shariff K. 1996. A high-resolution hybrid compact-ENO scheme for Shock-Turbulence interaction problems. J. Comput. Phys., 127(1): 27-51. DOI:10.1006/jcph.1996.0156 |
Alterman Z, Karal F C Jr. 1968. Propagation of elastic waves in layered media by finite difference methods. Bull. Seismol. Soc. Am., 58(1): 367-398. |
Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys., 114(2): 185-200. DOI:10.1006/jcph.1994.1159 |
Chen J B. 2006. Modeling the scalar wave equation with Nystr m methods. Geophysics, 71(5): T151-T158. DOI:10.1190/1.2335505 |
Cheng J B, Kang W, Wang T F. 2013. Description of qP-wave propagation in anisotropic media, Part Ⅰ:Pseudo-pure-mode wave equations. Chinese Journal of Geophysics (in Chinese), 56(10): 3474-3486. |
Chu C L, Phillips C, Stoffa P L. 2010. Frequency domain modeling using implicit spatial finite difference operators. SEG Technical Program Expanded Abstracts: 3076-3080. DOI:10.1190/SEGEAB.29 |
Chu P C, Fan C W. 1998. A three-point combined compact difference scheme. J. Comput. Phys., 140(2): 370-399. DOI:10.1006/jcph.1998.5899 |
Dennis S C R, Hudson J D. 1989. Compact h4 finite-difference approximations to operators of Navier-Stokes type. J. Comput. Phys., 85(2): 390-416. DOI:10.1016/0021-9991(89)90156-3 |
Dong L G, Ma Z T, Cao J Z, et al. 2000a. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese Journal of Geophysics (in Chinese), 43(3): 411-419. DOI:10.3321/j.issn:0001-5733.2000.03.015 |
Dong L G, Ma Z T, Cao J Z. 2000b. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation. Chinese Journal of Geophysics (in Chinese), 43(6): 856-864. DOI:10.3321/j.issn:0001-5733.2000.06.015 |
Du Q Z, Li B, Hou B. 2009. Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme. Appl. Geophys., 6(1): 42-49. DOI:10.1007/s11770-009-0008-z |
Du Q Z, Guo C F, Gong X F. 2015. Hybrid PS/FD numerical simulation and stability analysis of pure P-wave propagation in VTI media. Chinese Journal of Geophysics (in Chinese), 58(4): 1290-1304. DOI:10.6038/cjg20150417 |
Fletcher R, Du X, Fowler P J. 2008. A new pseudo-acoustic wave equation for TI media.//Proceedings of 2008 SEG Annual Meeting. Las Vegas, Nevada: SEG, 2082-2086. http://www.mendeley.com/catalog/new-pseudoacoustic-wave-equation-ti-media/
|
Gao Y J, Zhang J H, Yao Z X. 2016. Third-order symplectic integration method with inverse time dispersion transform for long-term simulation. J. Comput. Phys., 314: 436-449. DOI:10.1016/j.jcp.2016.03.031 |
Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bull. Seismol. Soc. Am., 86(4): 1091-1106. |
He X J, Yang D H, Wu H. 2015. A weighted Runge-Kutta discontinuous Galerkin method for wavefield modelling. Geophys. J. Int., 200(3): 1389-1410. DOI:10.1093/gji/ggu487 |
Huang C, Dong L G. 2009. Staggered-grid high-order finite-difference method in elastic wave simulation with variable grids and local time-steps. Chinese Journal of Geophysics (in Chinese), 52(11): 2870-2878. |
Jing H, Chen Y S, Yang D H, et al. 2017. Dispersion-relation preserving stereo-modeling method beyond Nyquist frequency for acoustic wave equation. Geophysics, 82(1): T1-T15. DOI:10.1190/geo2016-0104.1 |
Kosloff D, Pestana R C, Tal-Ezer H. 2010. Acoustic and elastic numerical wave simulations by recursive spatial derivative operators. Geophysics, 75(6): T167-T174. DOI:10.1190/1.3485217 |
Lele S K. 1992. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys., 103(1): 16-42. DOI:10.1016/0021-9991(92)90324-R |
Li X F, Wang W S, Lu M W, et al. 2012. Structure-preserving modelling of elastic waves:A symplectic discrete singular convolution differentiator method. Geophys. J. Int, 188(3): 1382-1392. DOI:10.1111/gji.2012.188.issue-3 |
Li Z C, Yang F S, Wang X D. 2016. High-precision numerical simulation of first-order qP-waves in the transversely isotropic (TI) medium optimized by the LS-RSGFD method. Chinese Journal of Geophysics (in Chinese), 59(4): 1477-1490. DOI:10.6038/cjg20160428 |
Liu S L, Yang D H, Lan C, et al. 2017. Modified symplectic schemes with nearly-analytic discrete operators for acoustic wave simulations. Comput. Phys. Commun., 213: 52-63. DOI:10.1016/j.cpc.2016.12.002 |
Liu Y, Sen M K. 2009a. A practical implicit finite-difference method:Examples from seismic modelling. J. Geophys. Eng., 6(3): 231-249. DOI:10.1088/1742-2132/6/3/003 |
Liu Y, Sen M K. 2009b. An implicit staggered-grid finite-difference method for seismic modelling. Geophys. J. Int., 179(1): 459-474. DOI:10.1111/gji.2009.179.issue-1 |
Liu Y. 2013. Globally optimal finite-difference schemes based on least squares. Geophysics, 78(4): T113-T132. DOI:10.1190/geo2012-0480.1 |
Liu Y. 2014. An optimal 5-point scheme for frequency-domain scalar wave equation. J. Appl. Geophys., 108: 19-24. DOI:10.1016/j.jappgeo.2014.06.006 |
Liu Y S, Liu S L, Zhang M G, et al. 2012. An improved perfectly matched layer absorbing boundary condition for second order elastic wave equation. Progress in Geophysics (in Chinese), 27(5): 2113-2122. |
Ma X, Yang D H. 2017. A phase-preserving and low-dispersive symplectic partitioned Runge-Kutta method for solving seismic wave equations. Geophys. J. Int., 209(3): 1534-1557. DOI:10.1093/gji/ggx097 |
Pan G D, Abubakar A, Habashy M T. 2012. An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme. Geophys. J. Int., 188(1): 211-222. DOI:10.1111/gji.2012.188.issue-1 |
Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics, 69(2): 583-591. DOI:10.1190/1.1707078 |
Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion, 31(1): 77-92. DOI:10.1016/S0165-2125(99)00023-2 |
Song G J. 2008. The improved nearly analytic discrete method for 3D elastic equations and its numerical simulation [Master's thesis] (in Chinese). Beijing: Tsinghua University.
|
Tang J, Wang F, Liu F L. 2016. 3D wave equation forward modeling based on three-level parallel acceleration. Oil Geophysical Prospecting (in Chinese), 51(5): 1049-1054. |
Tong P, Yang D H, Hua B L, et al. 2013. A high-order stereo-modeling method for solving wave equations. Bull. Seismol. Soc. Am., 103(2A): 811-833. DOI:10.1785/0120120144 |
Virieux J. 1986. P-SVwave propagation in heterogeneous media:Velocity-stress finite-difference method. Geophysics, 51(4): 889-901. DOI:10.1190/1.1442147 |
Wang L C, Chang X, Wang Y B. 2016. Forward modeling of pseudo P waves in TTI medium using staggered grid. Chinese Journal of Geophysics (in Chinese), 59(3): 1046-1058. DOI:10.6038/cjg20160325 |
Wang S Q, Yang D H, Yang K D. 2002. Compact finite difference scheme for elastic equations. Journal of Tsinghua University (Science and Technology) (in Chinese), 42(8): 1128-1131. |
Yang D H, Wang S Q, Zhang Z J, et al. 2003. n-times absorbing boundary conditions for compact finite-difference modeling of acoustic and elastic wave propagation in the 2D TI Medium. Bull. Seismol. Soc. Am., 93(6): 2389-2401. DOI:10.1785/0120020224 |
Yang D H, Wang N, Chen S, et al. 2009. An explicit method based on the implicit Runge-Kutta algorithm for solving wave equations. Bull. Seismol. Soc. Am., 99(6): 3340-3354. DOI:10.1785/0120080346 |
Yang D H, Wang L. 2010. A split-step algorithm for effectively suppressing the numerical dispersion for 3D seismic propagation modeling. Bull. Seismol. Soc. Am., 100(4): 1470-1484. DOI:10.1785/0120090200 |
Yang D H, Tong P, Deng X Y. 2012a. A central difference method with low numerical dispersion for solving the scalar wave equation. Geophysical Prospecting, 60(5): 885-905. DOI:10.1111/gpr.2012.60.issue-5 |
Yang D H, Wang N, Liu E. 2012b. A strong stability-preserving predictor-corrector method for the simulation of elastic wave propagation in anisotropic media. Communications in Computational Physics, 12(4): 1006-1032. DOI:10.4208/cicp.010111.230911a |
Yang D H, Wang M X, Ma X. 2014. Symplectic stereomodelling method for solving elastic wave equations in porous media. Geophys. J. Int., 196(1): 560-579. DOI:10.1093/gji/ggt393 |
Yang K D, Song G J, Li J S. 2011. FCT compact difference simulation of wave propagation based on the Biot and the squirt-flow coupling interaction. Chinese Journal of Geophysics (in Chinese), 54(5): 1348-1357. DOI:10.3969/j.issn.0001-5733.2011.05.024 |
Zhang J H, Yao Z X. 2013a. Optimized finite-difference operator for broadband seismic wave modeling. Geophysics, 78(1): A13-A18. DOI:10.1190/geo2012-0277.1 |
Zhang J H, Yao Z X. 2013b. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm. J. Comput. Phys., 250: 511-526. DOI:10.1016/j.jcp.2013.04.029 |
Zhang W S. 2006. Finite Difference Methods for Partial Differential Equations in Science Computation. Beijing: Higher Education Press.
|
Zhou B, Greenhalgh S A. 2003. Crosshole seismic inversion with normalized full-waveform amplitude data. Geophysics, 68(4): 1320-1330. DOI:10.1190/1.1598125 |
Zhou Y J, Yang D H, Ma X, et al. 2015. An effective method to suppress numerical dispersion in 2D acoustic and elastic modelling using a high-order Padé approximation. J. Geophys. Eng., 12(1): 114-129. DOI:10.1088/1742-2132/12/1/114 |
程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述Ⅰ:伪纯模式波动方程. 地球物理学报, 56(10): 3474-3486. DOI:10.6038/cjg20131022 |
董良国, 马在田, 曹景忠, 等. 2000a. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报, 43(3): 411-419. DOI:10.3321/j.issn:0001-5733.2000.03.015 |
董良国, 马在田, 曹景忠. 2000b. 一阶弹性波方程交错网格高阶差分解法稳定性研究. 地球物理学报, 43(6): 856-864. DOI:10.3321/j.issn:0001-5733.2000.06.015 |
杜启振, 郭成锋, 公绪飞. 2015. VTI介质纯P波混合法正演模拟及稳定性分析. 地球物理学报, 58(4): 1290-1304. DOI:10.6038/cjg20150417 |
黄超, 董良国. 2009. 可变网格与局部时间步长的交错网格高阶差分弹性波模拟. 地球物理学报, 52(11): 2870-2878. |
李振春, 杨富森, 王小丹. 2016. 基于LS-RSGFD方法优化的横向各向同性(TI)介质一阶qP波高精度数值模拟. 地球物理学报, 59(4): 1477-1490. DOI:10.6038/cjg20160428 |
刘有山, 刘少林, 张美根, 等. 2012. 一种改进的二阶弹性波动方程的最佳匹配层吸收边界条件. 地球物理学进展, 27(5): 2113-2122. |
宋国杰. 2008.三维弹性波方程的改进近似解析离散化方法及波场模拟[硕士论文].北京: 清华大学.
|
唐佳, 王凡, 刘福烈. 2016. 三维波动方程正演的三级并行加速. 石油地球物理勘探, 51(5): 1049-1054. |
王璐琛, 常旭, 王一博. 2016. TTI介质的交错网格伪P波正演方法. 地球物理学报, 59(3): 1046-1058. DOI:10.6038/cjg20160325 |
王书强, 杨顶辉, 杨宽德. 2002. 弹性波方程的紧致差分方法. 清华大学学报(自然科学版), 42(8): 1128-1131. DOI:10.3321/j.issn:1000-0054.2002.08.037 |
杨宽德, 宋国杰, 李静爽. 2011. Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟. 地球物理学报, 54(5): 1348-1357. DOI:10.3969/j.issn.0001-5733.2011.05.024 |
张文生. 2006. 科学计算中的偏微分方程有限差分法. 北京: 高等教育出版社.
|