频率域弹性波场模拟中高阶有限差分法的精度对比与改进策略

2014-07-05沈金松李曦宁

中国石油大学学报（自然科学版） 2014年5期

马超,沈金松,2,3,李曦宁

(1.中国石油大学地球物理与信息工程学院,北京 102249;2.中国石油大学油气资源与探测国家重点实验室,北京 102249; 3.中国石油大学中国石油集团公司物探重点实验室,北京 102249)

马超1,沈金松1,2,3,李曦宁1

推导基于二维交错网格的二阶和四阶频率域弹性波有限差分算子,并结合最优化差分系数和质量加权平均方法压制数值各向异性。基于均匀各向同性介质比较不同差分算子的模拟精度和频散关系,并详细比较二阶算子和四阶算子的精度和误差量级。结果表明,当采用交错网格最优化四阶差分算子,单位横波波长内含3个及以上网格点时,可将模拟误差控制在2%以内。同时,相对于常规网格,基于交错网格的差分算子可用于对含流体介质模型的模拟。

频率域;交错网格;有限差分算子;数值各向异性

频域有限差分(FDFD)模拟方法由于其格式简单、对计算资源的需求相对较小等特点,已被广泛应用于地震模拟。Pratt和Worthington[1]与Pratt[2]最早用FDFD法进行井间地震层析成像研究。Pan等[3]将频率四阶精度交错网格[3-4]FD方法应用于声波方程的模拟,分析了盐丘、Marmousi等典型复杂介质模型的响应特征。为了取得与时间域有限差分(FDTD)相当的模拟精度[5],单位波长内需要更多的离散网格点,这使得其在实际应用中受到很大限制。但同时频域波场模拟也具有多方面的优势。例如,对黏性介质及含衰减介质的波场模拟更具灵活性,以及易于并行化等优点。为了克服其局限性, Jo等[6]用优化9点FD算子来近似表示拉普拉斯算子,同时借用了有限元方法中的质量加权平均法; Shin和Sohn[7]将加权平均的范围扩大到25个点,使每个波长内的离散点减至3个以下。在弹性波模拟中,Štekl和Pratt[8]在前人方法的基础上,采用旋转网格给出了混合9点FD格式;Min等[9]构造了25点加权平均FD算子,使每个横波波长内的离散点降至3个,且维持最低程度的频散。总体上,交错网格FDFD弹性波模拟的应用还相对较少,主要由于低阶FD格式的精度低,而高阶格式的计算效率差,难以适应大模型的模拟计算。另外,在实际应用中,当处理含流固耦合界面情形时,基于常规网格的优化策略与泊松比相关,难以精确模拟泊松比变化剧烈的介质。笔者首先导出交错网格二阶和四阶精度的FD格式,并比较二阶和四阶精度FD格式的数值误差范围,同时对高阶FD格式采用优化差分和质量加权平均的方法压制数值各向异性的影响,最后应用提出的FDFD算子计算含流体介质的弹性波响应,以考察频域高阶FD格式的有效性和精度。

1 弹性波动方程的FD算子

1.1 基于常规网格的FD算子

在2D直角坐标系中,以x方向水平向右,z方向垂直向下,均匀各向同性介质中的弹性波二阶双曲型位移方程在频率域可以表示为

式中,u和v分别为空间位移场的水平和垂直分量; ω为圆频率;ρ为介质的密度;λ和μ为Lame常数,在频率域模拟时可设置为复数,以用于黏弹性介质的模拟。

对式(1)的FD求解算子主要有两种:一种为9点FD算子[2][8];另一种为25点FD算子[9]。对波动方程的FD求解实际上就是构造如下3个FD算子:Dxx=∂2/∂x2,Dzz=∂2/∂z2和Dxz=∂2/∂x∂z。二阶位移方程的FD求解通常使用常规网格,如图1所示,其中位移分量处于网格顶点的位置。

最早被提出来的频域FD算子为Pratt[2]基于二阶FD离散推导得到的9点FD算子;Štekl和Pratt[8]在经典二阶FD算子的基础上综合了来自旋转坐标系下的9点FD算子,从而形成了由两套FD算子加权平均得到的新的9点FD算子,有效提高了差分模拟的精度。以上两种FD算子都是9点FD算子,即关于1个点的微分表达式展开后共耦合了9个点的信息,Min等[9]提出的25点加权平均FD算子在空间上耦合了25个点的信息,采用了更为复杂的加权平均方法,并对频散关系进行最优化,得到可进一步压制频散的差分格式,理论上可给出均匀介质中的精确FD模拟结果。

图1 xoz平面内弹性FD的网格配置Fig.1 Grid stencil for elastic finite-difference in xoz-plane

1.2 基于交错网格的FD算子

在时间域,利用广义胡克定律,可以将二阶位移方程改写成一阶双曲型速度－应力方程[10],均匀各向同性介质的情况下,在频率域可以表示为

对式(2)和(3)的FD求解可以在经典的交错网格中进行,如图1所示,主应力分配在网格顶点中,剪应力分配在网格中心,而速度场则分配在网格边界的中点。交错网格FD法对介质泊松比变化剧烈的地层有很好的适应性,因此是时间域进行复杂介质地震模拟以及井中问题的重要方法。但在频率域的弹性波问题应用上还不是很广泛,注意到式(2)和(3)中既有速度场也有应力场的信息,这样将导致自由度的成倍增长,如果以这样的形式直接求解会引起不必要的巨大浪费。因此,需要将上式改造成只包含速度场分量的表达式进行求解。

基于二阶差分近似推导非均匀介质情况下弹性波方程的FDFD格式。求解过程中采用PML吸收边界条件[3-4,11-12]吸收来自人工边界的虚假反射。在二维情况下,引入复拉伸坐标系:

其中,n表示法向方向单位矢量;l表示位置;dn为与方向、吸收层厚度和所处位置有关的函数。

对式(4)～(9)进行层层代入,即先将式(8)、(9)代入到式(6)、(7)中,写出[σxx]i+1,k、[σxx]i,k、[τxz]i+1/2,k+1/2和[τxz]i+1/2,k－1/2的表达式;再将这些表达式代入到式(5)中;最后将结果代入到式(4)中即可消除所有的应力项,形成的FD格式中只含有速度项[13]。用层层代入方式形成的FD格式是基于一阶导数合成的,所以对PML的定义明确。同时,FD格式中包含了非均质的信息,这点在四阶精度FD格式中更为显著,因此对模拟复杂介质更具优势。

通过展开FDFD算子可以构造成一个大型稀疏矩阵[14-15]S(ω),实现FD离散近似的单频弹性波场数值解可用如下线性方程系统表示:

其中,v(ω)为频率为f=ω/(2π)时的单频速度场;F为源项。

以一维速度vx为例,四阶FD格式可以写成如下的形式[16-18]:

为了使显示更为简洁,令ρx=ρi+1/2,k,ρz=ρi,k+1/2,并用X代替vx、Z代替vz,ξ代替Sx、ζ代替Sz,则最终得到的交错网格四阶FD格式可以表示为

式中,c1和c2为FD系数,容易看出,当c2=0时,该式退化为二阶差分算子。为了确定c1和c2,分别将[vx]i+1、[vx]i、[vx]i+2和[vx]i－1在[vx]i+1/2处用泰勒展开至3阶,再将这4个表达式代入到式(11)中,整理得到二元一次方程组,求解得出FD系数的初值为c1=9/8,c2=－1/24。

2 数值各向异性的压制

2.1 质量加权平均方法

对波动方程的FD离散将不可避免地产生一定的网格频散,为了最小化数值各向异性,除了提高FD精度、采用加权平均FD算子,还可以引用有限元方法中的集总表示法[19]。在均匀各向同性介质的情况下,采用质量加权平均的方法,与常规网格的质量加权方式[8-9]类似,以交错网格系统中的vx为例,可以将vx(i+1/2,k)近似表示为(i+1/2,k)附近节点场值的加权平均,其中距离相同的节点具有相同的加权系数。当使用四阶FD算子时可以采用

式中,当a3～a4为零时可用于二阶算子的方案。交错网格四阶FD算子实际上是一个27点FD算子,从式(13)可以看出,参与质量加权平均的点共有13个,同样可以参照25点FD算子[9],将式(13)扩展为更多点(29点)的质量加权平均,此时能得到更加精确的结果,在模拟大尺度模型时具有明显优势,但同时也会使得计算代价急剧上升。

2.2 频散关系与校正

Štekl和Pratt[8]和Min等[9]详细介绍了二阶位移方程的数值相速度和群速度计算方法,以一阶速度－应力方程为例,简单推导四阶精度FD算子在均匀各向同性介质中的频散关系,其中二阶算子可以看做是四阶算子的特例。设在全空间均匀介质中的平面谐波为vω=Vωeik·r,其中V=(Vx,Vz)为水平分量和垂直分量的幅度,k=(kx,kz)表示波数向量,r= (x,z)表示位置向量。将vω代入四阶精度差分表达式中可以得到

式中,M为质量加权平均算子;α和β为纵、横波速度;FXX、FXZ、FZX、FZZ、以及FX和FZ分别表示空间导数的FD算子,在均匀各向同性介质情况下可表示为

式中,当a5～a7为零时对应式(13);kx=2πKcos θ, kz=2πKsin θ为波数分量,为分析方便,采用正方形网格剖分。

为了使方程(14)有解,则必须满足左边的矩阵行列式为零,对纵、横波,通过相速度vph=ω/k,群速度vgr=∂ω/∂k,得到纵、横波的频散条件:

利用梯度法[6]或者高斯牛顿法[9]求极小值,即可得到一系列最优化系数。

对一般情况下使用i=1～4即式(13)即可得到精确的模拟结果。

2.3 频散分析

根据相速度和群速度的计算关系,分别考察数值相速度和群速度随介质泊松比、传播角度以及1/ Gs的变化关系,其中Gs为单位横波波长内所占的网格点数。图2～5分别显示了不同差分算子的纵、横波相速度和群速度的频散关系:Vpph/Vp、Vpgr/Vp、Vsph/Vs和Vsgr/Vs。

图2 常规二阶和交错网格二阶FD算子采用质量加权平均方法的频散关系Fig.2 Dispersion curves of conventional&staggered grid 2nd-order FD operators with lumped mass scheme

常规二阶与交错网格二阶FD算子的计算代价和内存需求相当,为了使频散关系达到最好,使用质量加权平均方法,其中常规网格的差分及质量加权系数可参考文献[8]给出的差分格式。图2给出了当介质泊松比σ为0.25时,常规二阶和交错网格二阶FD算子的频散关系。可以看出,常规二阶算子在压制纵波频散问题上优势明显,但对横波的频散压制效果很差,当Gs小于20个的时候横波群速度的误差达到5%以上,而在相同误差标准时,交错网格二阶算子能将Gs减少到6.25个。事实上,对实际模拟的介质,纵波速度总在横波速度的1.732倍以上,因此横波的频散问题更为关键,所以交错网格二阶算子的方案比常规网格二阶算子的方案稍佳。

图3 混合二阶FD算子采用不同优化方案时的频散关系Fig.3 Dispersion curves of mixed 2nd-order FD operator with different optimized scheme

图4 优化25点FD算子的频散关系Fig.4 Dispersion curves of optimized 25-point FD operator

图5 交错网格四阶FD算子使用质量加权平均方法前后的频散关系对比Fig.5 Dispersion curves of staggered grid 4th-order FD operator with different lumped mass schemes

对混合二阶FD算子,需要对Fp和Fs或者max

(Fp,Fs)求极小化得出最优化差分系数[8],其中

图3分别给出了泊松比为0.25时,基于纵波压制、横波压制以及最大误差压制得到的最终频散关系。相对于Štekl和Pratt给出的基于最大误差压制的结果(图3(c)),基于横波的压制效果最佳,因为纵波波长较大,其频散条件较容易满足。

图4给出了泊松比分别为0.25和0.48时25点FD算子的频散关系(最优化系数见文献[9]),由图中可见其理论精度很高,当Gs不小于3.3个时可将误差控制在1%以内。但实际上采用多角度加权FD算子会引起源点波场失真,后面的模拟中将会提到;且基于均匀介质推导的FD系数在强非均质地层中的应用效果也将变差。

图5给出了交错网格四阶FD算子采用不同质量加权方案时的频散关系,图5(a)对应原始四阶差分格式的频散关系,该方法可有效模拟复杂介质波场传播。同时可见当经过频散压制后(图5 (b))可以实现更高精度的模拟,当Gs不小于2.9个的时候可将误差控制在2%以内,而且当采用29点质量加权平均时(图5(c))可以达到更高的精度。

交错网格FD算子还具有两大优势:首先,最优化差分系数的选取只与Gs有关,与介质泊松比无关,事实上随着泊松比的增大,相对纵波速度也将增大。常规网格FD算子与介质泊松比相关,当泊松比较大时,这些差分系数也会有所不同,当泊松比接近极限时无法克服网格频散,如图4(b)所示;另外,基于交错网格的FD算子可以自动适应含流固耦合界面的波场模拟,无需加入更多边界条件。

3 精度对比

用均匀各向同性介质模型验证FD算子的精度,密度和纵、横波速度分别为2 000 kg/m3、3 000 m/s、1800 m/s。图6(a)为10 Hz单频波场(实部)解析解(使用水平力源),图6(b)和6(c)分别为二阶和四阶FD数值解与解析解的残差波场。FD算子的吸收边界作用显著,另外四阶算子明显比二阶算子的精度要高的多,将误差从10－1降低到10－3量级。图7显示了图6(a)中红线所示的过源点沿垂直方向上解析解、四阶和二阶FD结果的对比。在水平方向上纵波占优,两者的误差水平接近;但垂直方向上横波占优,二阶差分的误差达到10－1量级,而四阶差分的误差依然在10－3量级。

图6 单频10 Hz波场水平分量的实部以及解析解与数值解之间的波场残差Fig.6 Real part of 10Hz monochromatic vxwave field and misfits between analytic and numerical results

图7 过源点垂直方向上波场信号比较Fig.7 Comparison of signals across source point on vertical direction

为了比较不同FD算子的模拟精度,在水平、垂直和45°方向各设置接收器A、B和C,源和接收器组合形成一个600 m边长的正方形。使用水平集中力源,将最小Gs控制为3个,模拟0～75 Hz的单频响应,得到主频为25 Hz的Ricker子波信号。图8 (a)显示了A上接收到的常规二阶和交错网格二阶FD解,可见交错网格二阶算子比常规网格的二阶算子模拟精度稍高。图8(b)显示了接收器B上,分别对纵、横波和综合压制时计算得到的信号,可见这三支信号相差很大,其中基于横波的压制方案效果最好。图8(c)给出了接收器C上用25点FD算子和最优化四阶FD算子模拟结果与解析解的对比图。两者的模拟精度都较高,但由于25点FD算子在脉冲源点处有振幅失真,因而使得高频信号段产生畸变,最后在接收信号上产生虚假的抖动。

图8 接收器A、B和C上接收到的由不同FD方法计算得到的信号和解析解Fig.8 Analytic signals and results based on different FD schemes from receiver A,B and C

4 流-固耦合介质的模拟

为进一步考察四阶精度FD格式对流固耦合界面[20]模拟问题的适应性,设计图9所示的模型,固体背景介质的纵、横波速度分别为3 000 m/s、1 800 m/s,流体块的纵、横波速度分别为1 500 m/s和0,固体和流体介质的密度分别为2 000 kg/m3、1 000 kg/m3。

图10显示了采用爆炸震源激发得到的弹性波场(实部)。图10(a)～10(b)为5 Hz情况下的水平波场;图10(c)～10(d)对应30 Hz时的垂直波场。其中图10(a)和10(c)使用常规网格算子;图10(b)和10(d)使用交错网格算子。由于常规网格对流固耦合界面的适应性差,即使在低频5 Hz情况下,在流体块内部也出现明显的频散;而交错网格模拟结果中基本没有表现出频散特征。对30 Hz激发时的模拟结果,用常规网格模拟的结果出现严重的混乱,而交错网格的模拟结果则呈现出稳定的波场,在满足频散条件内给出了正确的模拟结果。

图9 含流体介质的纵波速度模型Fig.9 P-wave velocity model of medium containing fluid

图10 用基于不同网格的FD算子计算含流体介质模型的5和30 Hz波场Fig.10 Wave field of 5&30 Hz in medium containing fluid by FD schemes based on different grids

5 离散线性稀疏方程组的求解

用FD方法离散得到的大型稀疏线性方程组

(10)可通过迭代法[14-15]和直接法[8-9]求解。相比而言,Krylov子空间迭代解法所需内存较少,适用于大尺度模型的求解。但由于采用PML等吸收边界条件,方程(10)左端的算子矩阵S(ω)是非Hermitian非对称复数矩阵,为了适应非Hermitian矩阵的特点,采用一种基于Lanczos分解原理的双共轭梯度法——广义乘积型双共轭梯度法(GPBiCG)求解。

以图2和图5所示的频散曲线为参考,对密度和纵、横波速度分别为2000 kg/m3、3000 m/s、1800 m/s的模型建立间距为10 m的100×100均匀网格,水平集中力源置于网格中心。对该模型用交错网格四阶FD算子准确求解的最低频率可达0.01 Hz,其中最高频率由频散条件决定。

图11 不同频率下和不同差分格式下迭代算法的残差与迭代次数之间的关系Fig.11 Relation between residual error and iteration times in different frequencies and FD schemes

图11(a)显示了四阶FD算子对应求解0.5、25和50 Hz单频波场的收敛情况,其中这三种频率下求解的波场与解析解之间的误差都在1%以内,残差范数为

从图中可以看出,求解25 Hz波场时最先收敛, 50 Hz高频时迭代稍慢,在0.5 Hz低频时迭代最慢。

在实际应用中,根据模型尺寸和所需频率之间的关系,可以灵活选择不同的FD算子进行正演模拟。对一个较大尺度的模型,模拟相对低频时使用二阶算子即可得到准确结果,而模拟相对高频时则适宜使用高阶算子。图11(b)、(c)分别显示了对不含流体介质的情况下,用四种FD算子迭代求解25 Hz和50 Hz波场时的收敛情况。在25 Hz情况下,用这四种算子均能得到较精确的模拟结果,从图11(b)中可以看出,常规网格二阶算子最先收敛,而交错网格算子需经过较多步迭代才收敛;在50 Hz情况下,根据频散条件,二阶算子已经无法准确求解,此时25点算子和交错网格四阶算子能实现高精度模拟,从图11(c)中可见,交错网格四阶算子比25点FD算子先达到收敛条件,这意味着对同样的频率,就迭代解法的速度而言,交错网格四阶算子比25点FD算子更适用于求解相对大尺度的问题。

6 结论

(1)基于交错网格的二阶FD算子在精度上要稍高于基于常规网格的二阶FD算子。为了进一步提高模拟精度:在常规网格的情况下,可以将9点差分扩展为25点差分,并采用质量加权平均法来克服频散现象;在交错网格的情况下,可以将二阶精度算子提高到四阶精度,从而得到标准四阶FD算子和最优化四阶FD算子。

(2)在均匀各向同性介质中,仅利用每个横波波长内3个网格点,即可获得横波群速度的模拟误差小于2%的精度模拟结果,且基于交错网格的FD算子的最优化系数不受介质泊松比的影响,因此在模拟复杂介质模型时具有优势。同时,交错网格FD算子能自动适应流固界面切变模量突变为零的情形,从而可以直接用于模拟含流体介质或者处理流体环境中的弹性波问题。

(3)当使用迭代方法进行波场求解时,为了提高效率,在满足频散关系的情况下,模拟相对低频宜采用低阶精度FD算子,而模拟相对高频宜采用高阶FD算子;同时,就迭代解法的速度而言,交错网格四阶算子比25点FD算子更适用于求解相对大尺度的问题。

[1] PRATT R G,WORTHINGTON M H.Inverse theory applied to multisource cross-hole tomography,part 1:acoustic wave-equation method[J].Geophys Prospecting, 1990,38:287-310.

[2] PRATT R G.Frequency-domain elastic wave modeling by finite differences:a tool for crosshole seismic imaging [J].Geophysics,1990,55(5):626-632.

[3] PAN G D,ABUBAKAR A,HABASHY T M.An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme[J].Geophys J Int,2012,188:211-222.

[4] HUSTEDT B,OPERTO S,VIRIEUX J.Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling[J].Geophys J Int,2004, 157:1269-1296.

[5] 方刚,FOMEL S,杜启振.交错网格Lowrank有限差分及其在逆时偏移中的应用[J].中国石油大学学报:自然科学版,2014,38(2):44-51.

FANG Gang,FOMEL S,DU Qizhen.Lowrank finite difference on a staggered grid and its application on reverse time migration[J].Journal of China University of Petroleum(Edition of Natural Science),2014,38(2): 44-51.

[6] JO C H,SHIN C,SUH J H.An optimal 9 point finite difference frequency-space 2-D wave extrapolator[J].Geophysics,1996,61(2):529-537.

[7] SHIN C,SOHN J H.A frequency-space 2-d scalar wave extrapolator using extended 25-point finite-difference operator[J].Geophysics,1998,63(1):289-296.

[8] ŠTEKL I,PRATT R G.Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators[J].Geophysics,1998,63(5):1779-1794.

[9] MIN D J,SHIN C,KWON B D,et al.Improved frequency-domain elastic wave modeling using weighted-averaging difference operators[J].Geophysics,2000,65 (3):884-895.

[10] VIRIEUX J.P-SV wave propagation in heterogeneous media,velocity stress finite difference method[J].Geophysics,1986,51(4):889-901.

[11] BÉRENGER J P.A perfectly matched layer for the absorption of electromagnetic waves[J].J Comput Phys, 1994,114:185-200.

[12] COLLINO F,TSOGKA C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics,2001,66(1):294-307.

[13] LUO Y,SCHUSTER G.Parsimonious staggered grid finite differencing of the wave equation[J].Geophys Res Lett,1990,17(2):155-158.

[14] VAN D V.ICCG and related methods for 3D problems on Vector Computers[J].Comput Phys Comm,1989, 53:223-235.

[15] VAN D V.Bi-CGSTAB:A fast and smoothly converging variant of Bi-CG for the solution of non-symmetric linear systems[J].Siam J Sci Stat Comput,1992,13:631-644.

[16] DABLAIN M A.The application of higher order differencing to the scalar wave equation[J].Geophysics, 1986,51(1):54-66.

[18] 裴正林.层状各向异性介质中横波分裂和再分裂数值模拟[J].石油地球物理勘探,2006,41(1):17-25.

PEI Zhenglin.Numeric simulation of S-wave splitting and second splitting in layered anisotropic media[J].OGP,2006,41(1):17-25.

[19] MARFURT K J.Accuracy of finite-difference and finiteelement modeling of the scalar and elastic wave-equation [J].Geophysics,1984,49(5):533-549.

[20] MA C,SHEN J S.2D elastic finite-difference modeling of seismic response in fractured anisotropic media[C].75th EAGE Conference,Expanded Abstract,2013.

(编辑修荣荣)

Accuracy comparison and improvement strategy in numerical modeling of elastic wave in frequency-domain by high-order finite-difference scheme

MA Chao1,SHEN Jinsong1,2,3,LI Xining1
(1.College of Geophysics and Information Engineering in China University of Petroleum,Beijing 102249,China; 2.State Key Laboratory of Petroleum Resources and Prospecting in China University of Petroleum,Beijing 102249,China; 3.CNPC Key Lab of Geophysical Exploration in China University of Petroleum,Beijing 102249,China)

2D finite-difference operators in second-order as well as fourth-order accuracy based on staggered grid in frequency-domain are derived.By combining optimal difference coefficients with lumped mass methods with a weighted average,numerical anisotropy can be suppressed,which enables accurate dispersion comparisons among these difference operators.U-sing a homogeneous isotropic medium,simulation precisions of different operators are compared to contrast the accuracy and error magnitude between second-order and fourth-order operators.It appears that the resulting error is smaller than 2%if the number of grid points required per smallest shear wavelength increases to 3.Meanwhile,in contrast to those based on conventional grid,difference operators based on staggered grid can be conducted in model simulations in fluids.

frequency-domain;staggered-grid;finite-difference operator;numerical anisotropy

P 631.4

1673-5005(2014)05-0048-11

10.3969/j.issn.1673-5005.2014.05.007

2014－04－28

国家“973”重点基础研究发展计划(2013CB228605);国家自然科学基金项目(41374141)

马超(1989－),男,博士研究生,主要从事地震和电磁模拟及反演研究。E－mail:oilmachao@gmail.com。

马超,沈金松,李曦宁.频率域弹性波场模拟中高阶有限差分法的精度对比与改进策略[J].中国石油大学学报:自然科学版,2014,38(5):48-58.

MA Chao,SHEN Jinsong,LI Xining.Accuracy comparison and improvement strategy in numerical modeling of elastic wave in frequency-domain by high-order finite-difference scheme[J].Journal of China University of Petroleum(Edition of Natural Science),2014,38(5):48-58.