地球物理学报  2016, Vol. 59 Issue (7): 2650-2662   PDF    
引用本文
郑婉秋, 孟小红, 刘建红 等. 2016. 基于余弦调制Chebyshev窗的弹性波高精度正演. 地球物理学报, 59(7): 2650-2662, doi: 10.6038/cjg20160728.  
ZHENG W Q, MENG X H, LIU J H, et al. 2016. High precision elastic wave equation forward modeling based on cosine modulated Chebyshev window function. Chinese J. Geophys. (in Chinese), 59(7): 2650-2662, doi: 10.6038/cjg20160728.  
基于余弦调制Chebyshev窗的弹性波高精度正演
郑婉秋1 , 孟小红1 , 刘建红2 , 王建1     
1. 中国地质大学(北京)地球物理与信息技术学院, 北京 100083;
2. 中国石油集团东方地球物理勘探有限责任公司物探技术研究中心, 河北 涿州 072751
收稿日期 2015-11-30, 2016-06-11 收修定稿
基金项目: 国家重大科研装备研制项目(ZDYZ2012-1-02-04)和国家自然科学基金(41474106)联合资助.
作者简介: 郑婉秋,在读硕士,主要从事复杂介质中地震波模拟及偏移成像方面的研究.E-mail:453294543@qq.com
通讯作者: 孟小红,教授,主要从事计算与综合地球物理研究.E-mail:mxh@cugb.edu.cn
摘要: 有限差分时间域正演是弹性波逆时偏移和全波形反演的基础,正演的计算精度也控制着偏移结果的准确性,若精度不高,则在偏移、反演后会带来假象.为了有效提高正演精度,本文结合窗函数优化方法,在窗函数截断伪谱法空间褶积序列以逼近有限差分算子的基础上,提出了一种基于Chebyshev窗的余弦调制模型,在原始Chebyshev窗的基础上引入了调制次数和调制范围,通过调节这两个参数可以人工可视化的调节截断误差,新的窗函数继承了Chebyshev窗的特点,在不明显降低截断谱范围的基础上明显降低了截断误差.本文针对不同正演阶数N,给出了一组经验调制系数,并通过数值模拟方法,对比了新方法、改进二项式窗和基于最小二乘优化方法的正演效果.结果表明,基于余弦调制的Chebyshev窗控制数值频散的能力更强,在大网格下可以得到更精确的正演结果.从经济角度分析,该方法减小了计算花费,提高了计算效率.
关键词: 有限差分      数值频散      窗函数      弹性波      余弦调制     
High precision elastic wave equation forward modeling based on cosine modulated Chebyshev window function
ZHENG Wan-Qiu1, MENG Xiao-Hong1, LIU Jian-Hong2, WANG Jian1     
1. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China;
2. BGP Research and Development Center, CNPC, Hebei Zhuozhou 072751, China
Abstract: The finite difference forward modeling is the basis of elastic wave reverse-time migration and full waveform inversion in the time domain. The accuracy of forward modeling also controls the accuracy of seismic imaging and inversion. The migration or inversion will bring illusion if the accuracy is not high. We can get optimized explicit finite difference operators by using the window function to truncate spatial convolution counterpart of the pseudo-spectral method. Based on this, a cosine modulated Chebyshev window is designed. On the basis of the original Chebyshev window, the modulation times and modulation domain are introduced, and we can adjust truncation error visually by controlling these two parameters. As the new window function inherits the character of Chebyshev window, we observe that the spectral range using the modulated window function for truncation is significantly broader than using the conventional window function with stable error. For different forward modeling orders N, we give a set of empirical modulation factors and compare the forward modeling effect of the new method and improved binomial window by the numerical simulation method. The results demonstrate that the operators based on the cosine modulated Chebyshev window can efficiently suppress the numerical dispersion and get more accurate forward modeling results on the large grid. From economic perspective, this method reduces the computational cost and improves efficiency..
Key words: Finite difference      Numerical dispersion      Window function      Elastic wave      Cosine modulation     
1 引言

弹性波逆时偏移是依据弹性波动理论，以弹性波方程为基础的逆时偏移方法，弹性矢量波场逆时偏移技术采用全波波动方程，可获得更加丰富的细节信息，不论在各向同性介质还是各向异性介质，弹性波对地下介质的成像效果要优于单纯声波的成像 效果(Chang and McMechan，1987; Sun and McMechan，2001). 而弹性波正演是弹性波逆时偏移的基础，准确的正演可以避免成像假象，提高成像精度.常用的数值模拟分为三类：几何射线法、积分方程法和波动方程法，波动方程法的模拟结果包含了地震波的运动学和动力学特征，是最常用的(程冰洁等，2008).

较多使用的时域模拟方法有伪谱法(Gazdag，1981;Kosloff and Baysal，1982;Carcione et al.，1992)和有限差分方法(Alterman and Karal，1968;Kelly et al.，1976;Dablain，1986;Igel et al.，1995).

伪谱法的计算精度较高，在时间轴上其采用有限差分方法进行近似，而在空间方向采用傅里叶变换求解，理论上其利用了计算区域的所有点，这样的算法避免了网格频散的现象，但是大大增加了正演 耗时(Gazdag，1981；Kosloff and Baysal，1982； Carcione et al.，1992).有限差分方法即在时间方向和空间方向皆采用差分近似微分的方法，自20世纪60年代起就被地球物理学家广泛使用，逐步发展出常规网格有限差分、普通交错网格有限差分、旋转交错网格有限差分.Alterman和Karal(1968)首次利用有限差分法进行弹性波正演.Madariaga(1976)首次实现了基于交错网格的弹性波数值模拟.Saenger等(2002，2004)提出了适应于各向异性介质的旋转交错网格有限差分方法.刘洋等(1998)实现了任意偶数阶弹性波的正演模拟.董良国等(2000)实现了一阶弹性波的交错网格有限差分.但是对于有限差分数值模拟来说，数值频散是一个不可回避的问题.随着震源子波主频的升高，空间采样间隔增大，数值频散会愈发严重.根本上来讲，利用有限点的差分近似代替连续的微分算子，在截断后就产生了数值频散.降低数值频散的简单办法就是采用更细致的网格并降低子波主频，但是这会带来海量的计算.前人在如何降低数值频散方面做了诸多工作，其一是通量校正技术(FCT)，自Boris和Book(1973)首先提出通量校正技术(FCT)后，Fei和Larner(1995)将该技术推广到弹性波正演，随后杨顶辉和腾吉文(1997)将FCT技术运用到各向异性介质的模拟中.其二是通过最优化方法，如模拟退火法(Zhang and Yao，2013)、最小二乘法(Yang，2014)等.再有就是本文所利用的窗函数法，据Fornberg(1987)的研究，伪谱法可以表示为有限差分法的高阶近似，换而言之，即利用时间域窗函数可以截断伪谱法序列以得到有限差分系数.前人利用了不同的窗函数以得到优先 差分系数，如加权Hanning窗(Zhou and Greenhalgh，1992)、高斯窗(Igel，1995)、改进二项式窗(Chu and Stoffa，2012)、基于自褶积的Chebyshev窗(王之洋等，2015).

在前人的研究基础上，本文提出了一种基于余弦调制的Chebyshev窗，即将Chebyshev窗与一定范围的余弦函数相乘数次.这里的范围即调制范围，次数即调制次数.研究表明，经过余弦调制后的Chebyshev窗提高了原有窗函数的截断效果，在不明显降低截断谱范围的基础上，大大提高了截断误差的稳定性.通过实验分析，本文给出了针对不同阶数的调制参数，并通过数值模拟证明，该方法可适用于较大网格的正演模拟，且在压制数值频散方面较其他窗函数和优化方法有较大的优势，从经济上讲，我们可以利用改进后的低阶模拟代替原始高阶计算，大大降低了计算花费.

2 窗函数逼近法 2.1 差分算子及截断误差

据前人研究表明，有限差分算子可以通过sinc函数插值理论推导得到，sinc函数可以重构均匀采样后的带限信号(Diniz et al.，2012)，sinc函数可以表示为

(1)

其中，Δx为空间采样间隔，代表截止波数.Chu和Stoffa(2012)利用窗函数截断公式(1)的一阶导数和二阶导数，获取了有限差分算子.对公式(1)求导后可表示为

(2)

(3)

截断后，可以表示为

(4)

(5)

因n=0为公式(2)和公式(3)的一个奇异点，根据Lee和Seo(2002)，公式(2)和公式(3)可以表示为

(6)

(7)

其中dn2和dn1为fn的系数，这里ξ代表黎曼ξ函数.加窗截断后可以表示为

(8)

(9)

这里，

将公式(8)和公式(9)变换到频率域后为

(10)

(11)

对于一阶导数，可以将误差函数表示成

(12)

对于二阶导数，可以将误差函数表示成

(13)

2.2 传统窗函数下的逼近

如今信号领域广泛使用的窗函数包括：Hanning窗、 Kaiser窗、Chebyshev窗等，而在利用窗函数截断逼近有限差分算子方面，也有不同学者使用了不同的窗函数，如：广义加权的Hanning窗(Zhou and Greenhalgh，1992)、高斯窗(Igel，1995)、基于自褶积的Chebyshev窗(王之洋，2015)和本文将重点对比的改进二项式窗(Chu and Stoffa，2012)，二项式窗的窗函数形如公式(14).

(14)

Chu和Stoffa认为，当N=N+M时，窗函数的形式会发生改变，其截断效果更佳，可以更好地控制数值频散，而当M=0时，截断参数为传统的有限差分算子系数，从这个角度来说，改进二项式窗是有限差分法和伪谱法连接的桥梁.

在空间域，将伪谱法的空间褶积序列截断就得到了有限差分法，如果从信号处理层面来讲，截断相当于加了一个窗函数(王之洋等，2015).这就是说，我们可以将时间域的窗函数看成一个有限长的滤波器，滤波器主瓣宽度和旁瓣衰减也就决定了其性质.Diniz等(2012)的分析结论表明，主瓣宽度与截断的谱范围呈负相关关系，即主瓣越宽，则截断的谱范围越低.而旁瓣衰减与截断误差的稳定性呈正相关关系，即旁瓣衰减越大，误差稳定性越高.

图 1是不同窗函数在N=24时的时间域展布，图 2为各窗函数的幅值响应.从图中可以发现，当时间域窗函数系数变化率较高时，其幅值响应有着更宽的主瓣，而窗函数尾端曲率越大，则旁瓣衰减越明显.在图 1所示各种窗函数中，改进二项式窗系数有着最快的下降速率，也对应了最宽的主瓣范围，但也拥有最大的旁瓣衰减.对于Chebyshev窗，当参数γ增大时，其主瓣将变宽缓，旁瓣衰减变大.Kaiser窗与Chebyshev窗有着类似的性质.

图 1 窗函数对比 Fig. 1 Comparison of window functions
图 2 各窗函数的幅值响应 Fig. 2 Amplitude responses of window functions

综合以上分析可以得出结论：若想有效地控制数值频散，应该寻求主瓣较窄、旁瓣衰减较大的窗函数截断褶积序列以得到有限差分系数.图 3显示了图 1中所示各种传统窗函数在二阶导数下的截断误差，可以看出，Kaiser窗有着最大的截断谱范围，但从图 3b中看出，其截断误差的稳定性是几种窗函数中最差的，这一点也印证了我们之前对主瓣范围和旁瓣衰减在误差控制方面的分析.在另外三种窗函数中，改进二项式窗的误差稳定性最好，但是谱范围也是最小的，这就意味着其对高波数部分的逼近效果较差，会产生较强的数值频散.Hanning窗与Chebyshev窗的谱范围大致相同，其在低波数范围的误差较小，但是其误差随波数的增加而迅速增加，在高波数部分表现出了强烈的动荡，从这个角度来讲，虽然Chebyshev窗在整个区域内都存在着误差，但其较为稳定，不存在误差的急剧变化，在这几种传统窗函数中为较优的选择.

图 3 传统窗函数二阶导数下的截断误差(a)截断误差对比；(b)将(a)纵轴放大1000倍结果. Fig. 3 Truncation errors for second derivative of traditional window function(a)Comparison of truncation errors；(b)The longitudinal axis of(a)amplified 1000 times.
2.3 余弦调制Chebyshev窗下的逼近

所谓余弦调制，即将Chebyshev窗与一定范围的余弦函数在时间域相乘数次，可以表示为

(15)

其中，wdcs代表余弦调制后的窗函数，wdche表示Chebyshev窗函数，e表示调制次数，p表示调制范围，N+1表示窗函数长度.可以看出，我们在余弦调制的过程中定义了两个参数，即调制范围和调制次数.

图 4显示了不同调制参数下窗函数的幅值响应，可以看出，当p值增大，e值不变时，窗函数的主瓣宽度逐渐减小，旁瓣衰减变化不大.而当e值增大，p值不变时，主瓣宽度有所增加，但旁瓣衰减也逐步增大.也就是说，可以通过对这两个参数的动态选择，以控制所调制窗函数的主瓣宽度和旁瓣衰减，最终达到优化原有Chebyshev窗、优化有限差分系数的目的.以24阶有限差分算子为例，计算了不同调制范围和调制次数下的二阶导数截断误差，如图 5.可以看出，调制范围相同，调制次数增加时，截断谱范围降低、误差稳定性增强.而当调制次数不变，调制范围增大时，误差稳定性降低、谱范围增大.所以调制次数和调制范围为两个相互制约且起反作用的参数.为平衡两者，本文采用的做法是对不同阶数进行多次误差分析和正演模拟试验，以选择较为合适的经验调制参数，表 1为通过实验分析得到的不同阶数下截断效果较好的调制范围参数和调制次数参数.

图 4 不同调制参数下的窗函数幅值响应 Fig. 4 Amplitude responses for different modulation parameters of window function
图 5 不同调制参数下的截断误差 Fig. 5 Truncation errors for different modulation parameters
表 1 各阶数调制参数分布表 Table 1 Distribution of modulating parameters using different orders

通过实验分析，在N+1=25时确定p=3、e=3有着较优的截断效果，而N+1=9时确定p=3.5、e=1有较优的截断效果.图 5为不同调制参数下一阶导数的截断误差.

从图 6中可以看出，对于高阶数N=24时，虽然调制后的窗函数牺牲了一定的谱范围，但是余弦调制Chebyshev窗在控制误差稳定性上有了很大的提升，通过和改进二项式窗(M=16)的对比可见，新改造的窗函数对误差稳定性的控制已不亚于改进二项式窗，而且有更宽的谱范围.图 7显示了对于低阶N=8时，几种窗函数的截断效果，由图中可以看出，改进二项式窗(M=16)和Chebyshev窗出现了不稳定的情况，而此时的余弦调制窗函数则依旧保持了误差的稳定性，其截断的谱范围也没有明显降 低.综合图 3、图 5及图 6分析并结合正演阶数、截断谱范围和误差稳定性三方面综合考虑，基于余弦调制Chebyshev窗是最优的选择.表 2及表 3分别是一阶导数、二阶导数下的有限差分系数，图 8中整合了不同阶数下优化有限差分算子和基于自褶积Chebyshev窗(王之洋等，2015)截断的有限差分算子的二阶导数截断误差，从图中不难发现，余弦调制 后的Chebyshev窗有更好的误差稳定性.

图 6 调制窗函数截断误差效果(N=24)(a)截断误差对比；(b)将(a)纵轴放大1000倍结果. Fig. 6 Truncation errors of modulating window function(N=24)(a)Comparison of truncation errors；(b)The longitudinal axis of(a)amplified 1000 times.
图 7 调制窗函数截断误差效果(N=8) (a)截断误差对比；(b)将(a)纵轴放大1000倍结果. Fig. 7 Truncation errors of modulating window function(N=8) (a)Comparison of truncation errors；(b)The longitudinal axis of(a)amplified 1000 times.
表 2 一阶导数下的有限差分系数 Table 2 The finite difference coefficients for the first order derivative
表 3 二阶导数下的有限差分系数 Table 3 The finite difference coefficients for the second order derivative
图 8 余弦调制Chebyshev窗和自褶积Chebyshev窗二阶导数截断误差对比 Fig. 8 Comparison of truncation errors between modulated Chebyshev window and auto-convolution Chebyshev window
3 模型测试

对于有限差分法，正演时要对时间轴和空间轴分别进行差分表示，再按照时间轴方向递推，据此本文进行了基于应变-位移方程的各向同性和各向异性弹性波模拟，其中时间轴二阶导数可以表示为

(16)

(17)

对于空间轴，本文采用常用的2N阶离散格式：

(18)

(19)

在公式(16)—(19)中，ux和uz代表离散的时间域波场，am为有限差分系数，Δx和Δz代表空间步长.

可以说，弹性波动理论的基础是本构方程、运动微分方程和几何方程，想得到最终的离散正演公式，只需将公式(16)—(19)带入这三个方程即可，以x方向为例，最终的离散公式可以表示为

(20)

3.1 简单模型

对于均匀模型的弹性波数值模拟是验证窗函数截断所得到有限差分系数质量的有力手段，我们以40 Hz Ricker子波为震源，将震源置于模型中间，模型纵 波速度设定为2000 m·s-1，横波速度为1300 m·s-1，模型的空间采样间隔设定为10 m，网格大小为600×600，并采用海绵吸收边界改善模型边界反射.本文测试了改进二项式窗8阶、余弦调制窗8阶、改进二项式窗12阶算子在各向同性(ISO)和各向异性(VTI)介质模型下的正演效果，如图 9，数值验证表明：余弦调制窗8阶的ISO介质正演效果远远超过了改进二项式窗8阶的正演效果，其控制数值频散的能力甚至超过了改进二项式窗12阶算子，而新方法8阶算子的VTI正演效果明显优于改进二项式窗12阶算子，如图 10.数值模拟的结果与误差谱范围的覆盖及误差稳定性分析结果有着统一的结论，在误差允许范围内，谱范围越大，误差稳定性越强，控制数值频散的能力越好.这就意味着较低阶的有限差分算子可以代替高阶有限差分算子，在控制误差的同时，大大节约了计算成本.

图 9 均匀ISO介质测试 (a)8阶改进二项式窗波场快照X分量;(b)8阶改进二项式窗波场快照Z分量;(c)8阶余弦调制窗波场快照X分量;(d)8阶余弦调制窗波场快照Z分量;(e)12阶改进二项式窗波场快照X分量;(f)12阶改进二项式窗波场快照Z分量. Fig. 9 Tests of homogeneous isotropic media (a)Using the 8th improved binomial window function(X component);(b)Using the 8th improved binomial window function(Z component);(c)Using the 8th modulated window function(X component);(d)Using the 8th modulated window function(Z component);(e)Using the 12th improved binomial window function(X component);(f)Using the 12th improved binomial window function(Z component).
图 10 均匀VTI介质测试 (a)12阶改进二项式窗波场快照X分量;(b)12阶改进二项式窗波场快照Z分量;(c)8阶余弦调制窗波场快照X分量;(d)8阶余弦调制窗波场快照Z分量. Fig. 10 Tests of homogeneous anisotropic media (a)Using the 12th improved binomial window function(X component);(b)Using the 12th improved binomial window function(Z component);(c)Using the 8th modulated window function(X component);(d)Using the 8th modulated window function(Z component).
3.2 BP模型

为了进一步测试基于余弦调制有限差分算子的稳定性，我们在复杂的BP模型上进行了弹性波数值模拟.对比了常规有限差分算子、基于最小二乘优化方法(Yang，2014)和本文余弦调制方法在正演时的稳定性.模型规模为3792×1896，空间采样间隔为6.25 m，时间采样间隔为0.5 ms，时间采样点数为14000.震源位于表层中间位置，Ricker子波主频为30 Hz，图 11和图 12分别为BP纵波速度模型和密度模型.

图 11 BP纵波速度模型 Fig. 11 BP longitudinal wave velocity model
图 12 BP密度模型 Fig. 12 BP density model

(21)

其中VS为横波速度，VP为纵波速度.

我们将正演算法并行化，采用的GPU型号为GTX960，该图形处理器拥有4 GB显存，192个计算核心，每个流处理器可最多容纳2048个线程.CPU为Inter(R)i7-4790k@4GHZ，计算机内存为8 G.表 4显示了相同炮点不同阶数的常规有限差分算子和优化有限差分算子在BP模型上的计算耗时，可以看出，优化后的有限差分算子并没有增加过多的计算花费.

表 4 各阶算子计算效率 Table 4 Computational efficiency using different orders

从正演道记录的局部细节可以清晰地看出，基于余弦调制的Chebyshev窗在抑制数值频散方面有突出的优势.由图 13d，13e，13f可见，在直达波下，传统方法出现了较明显的“拖尾”现象，而基于余弦调制的Chebyshev窗很好地控制了直达波下的数值频散，且控制频散的能力强于最小二乘优化方法.对比图 13g，13h，13i中有明显反射的局部细节，可以发现余弦调制方法可以使有效反射更加清晰，其抑制频散的能力也是这三种方法中最好的.

图 13 BP模型炮记录 (a)常规8阶算子；(b)最小二乘8阶算子；(c)调制窗8阶算子；(d)(e)(f)分别为(a)(b)(c)左侧红框相同位置的细节放大；(g)(h)(i)分别为(a)(b)(c)右侧红框相同位置的细节放大. Fig. 13 Shot records for the BP model (a)Conventional 8th operator;(b)Least squares 8th operator;(c)Modulated 8th operator;(d)(e)(f)Zoomed in left red box of(a)(b)(c);(g)(h)(i)Zoomed in right red box of(a)(b)(c).
4 结论

本文在窗函数截断伪谱法空间褶积序列的基础上，分析了现有的多种三角类窗函数和改进二项式窗函数，对比分析不同窗函数的截断效果后，提出了一种基于余弦调制的Chebyshev窗.从时间域来看，窗函数可以看成一个有限长的滤波器，滤波器主瓣宽度和旁瓣衰减决定了其性质.具体来讲，主瓣宽度控制了截断的频谱覆盖范围，主瓣越窄，截断的谱范围越宽，即负相关关系.而旁瓣衰减则控制了误差的稳定性，即偏差程度，旁瓣衰减越大，误差稳定性越好.通过对几种传统窗函数的截断误差分析我们发现，无论从截断的谱范围还是误差稳定性来看，Chebyshev窗都是一个较优的选择，但是实际上，其误差稳定性仍较改进二项式窗差.

据此，本文提出的基于余弦调制的Chebyshev窗在保证原始截断谱范围的基础上，大大增强了Chebyshev窗截断误差的稳定性，使其提升了两个数量级.在余弦调制的过程中，我们引入了两个参数，分别为调制范围和调制次数，调制范围主要控制了窗函数的主瓣范围，而调制次数影响了旁瓣衰减.经过对不同调制参数进行误差分析和数值模拟，本文总结了一套经验系数.

通过对均匀模型的试算，验证了在理论误差分析中得出的结论：余弦调制窗8阶的均匀介质正演效果远远超过了改进二项式窗8阶的正演效果，其控制数值频散的能力甚至超过了改进二项式窗12阶算子，而新方法8阶算子的VTI介质正演效果明显优于改进二项式窗12阶算子，对于BP模型的模拟也验证了本文方法的有效性.在保证误差稳定的前提下，可以以低阶数的运算代替较高阶数的运算，这减小了计算耗时，提高了经济效率.

参考文献
Alterman Z, Karal F C Jr. 1968. Propagation of elastic waves in layered media by finite difference methods. Bull. Seism. Soc. Am. , 58(1): 367–398.
Boris J P, Book D L. 1973. Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works. Journal of Computational Physics , 11(1): 38–69.
Carcione J M, Kosloff D, Behle A, et al. 1992. A spectral scheme for wave propagation simulation in 3-D elastic-anisotropic media. Geophysics , 57(12): 1593–1607. doi: 10.1190/1.1443227.
Chang W F, McMechan G A. 1987. Elastic reverse time migration. Geophysics , 52: 1367–1375.
Cheng B J, Li X F, Long G H. 2008. Seismic waves modeling by convolutional Forsyte polynomial differentiator method. Chinese J. Geophys. (in Chinese) , 51(2): 531–537.
Chu C L, Stoffa P L. 2012. Determination of finite-difference weights using scaled binomial windows. Geophysics , 77(3): W17–W26. doi: 10.1190/GEO2011-0336.1.
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.
Diniz P S R, da Silva E A B, Netto S L. Digital Signal Processing System Analysis and Design. Beijing: China Machine Press, 2012 .
Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese J. Geophys. (in Chinese) , 43(3): 411–419.
Fei T, Larner K. 1995. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport. Geophysics , 60(6): 1830–1842.
Fornberg B. 1987. The pseudospectral method: Comparisons with finite differences for the elastic wave equation. Geophysics , 52(4): 483–501. doi: 10.1190/1.1442319.
Gazdag J. 1981. Modeling of the acoustic wave equation with transform methods. Geophysics , 46(6): 854–859. doi: 10.1190/1.1441223.
Igel H, Mora P, Riollet B. 1995. Anisotropic wave propagation through finite-difference grids. Geophysics , 60(4): 1203–1216. doi: 10.1190/1.1443849.
Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms: A finite-difference approach. Geophysics , 41(1): 2–27. doi: 10.1190/1.1440605.
Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method. Geophysics , 47(10): 1402–1412. doi: 10.1190/1.1441288.
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 Y, Li C C, Mou Y G. 1998. Finite-difference numerical modeling of any even-order accuracy. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 33(1): 1–10.
Madariaga R. 1976. Dynamics of an expanding circular fault. Bull. Seism. Soc. Am. , 65: 163–188.
Saenger E H, Shapiro S. 2002. An effective velocities in fractured media: A numerical study using the rotated staggered finite-difference grid. Geophysical Prospecting , 50: 183–194.
Saenger E H, Thomas B. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rots-ted staggered grid. Geophysics , 69(2): 583–591.
Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics , 66(5): 1519–1527.
Wang Z Y, Liu H, Tang X D, et al. 2015. Optimized finite-difference operators based on Chebyshev auto-convolution combined window function. Chinese J. Geophys. (in Chinese) , 58(2): 628–642. doi: 10.6038/cjg20150224.
Yang D H, Teng J W. 1997. FCT finite difference modeling of three-component seismic records in anisotropic medium. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 32(2): 181–190.
Yang L, Yan H Y, Liu H. 2014. Least squares staggered-grid finite-difference for elastic wave modelling. Exploration Geophysics , 45: 255–260.
Zhang J H, Yao Z X. 2013. Optimized finite-difference operator for broadband seismic wave modeling. Geophysics , 78(1): A13–A18. doi: 10.1190/GEO2012-0277.1.
Zhou B, Greenhalgh S A. 1992. Seismic scalar wave equation modeling by a convolutional differentiator. Bull. Seism. Soc. Am. , 82(1): 289–303.
程冰洁, 李小凡, 龙桂华. 2008. 基于广义正交多项式褶积微分算子的地震波场数值模拟方法. 地球物理学报 , 51(2): 531–537.
董良国, 马在田, 曹景忠, 等. 2000. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报 , 43(3): 411–419.
刘洋, 李承楚, 牟永光. 1998. 任意偶数阶精度有限差分法数值模拟. 石油地球物理勘探 , 33(1): 1–10.
王之洋, 刘洪, 唐祥德, 等. 2015. 基于Chebyshev自褶积组合窗的有限差分算子优化方法. 地球物理学报 , 58(2): 628–642.
杨顶辉, 腾吉文. 1997. 各向异性介质中三分量地震记录的FCT有限差分模拟. 石油地球物理勘探 , 32(2): 181–190.