地球物理学报  2021, Vol. 64 Issue (7): 2480-2493   PDF    
基于多倍角公式的一步波场外推法
柯璇1,2, 石颖1,2,3,4, 王银凤5,6     
1. 东北石油大学地球科学学院, 大庆 163318;
2. 黑龙江省油气藏形成机理与资源评价重点实验室, 大庆 163318;
3. 东北石油大学非常规油气研究院, 大庆 163318;
4. 东北石油大学陆相页岩油气成藏与高效开发教育部重点实验室, 大庆 163318;
5. 吉林大学地球探测科学与技术学院, 长春 130021;
6. 东北石油大学数学与统计学院, 大庆 163318
摘要:为了提高地震波场正演模拟的准确性和稳定性,针对一步波场外推法地震波场正演,本文提出了基于多倍角公式的耦合方程组解法.借助欧拉公式,将一步波场外推法的复数波场延拓方程转化为两个实数波场耦合的方程组,结合多倍角公式和泰勒展开式精确逼近包含拟微分算子的简谐函数算子,利用谱方法求解拟微分算子,进而推导了一种基于多倍角公式的一步波场外推法的耦合方程组.相比于常规一步波场外推法中复数方程的矩阵解法,本文方法能够显著减少傅里叶变换次数,降低计算成本.此外,本文推导了稳定性条件,为正确选取地震波场模拟参数提供了理论依据.基于二维匀速模型和复杂构造模型的数值测试表明,本文方法能够在大时间步长情况下保持外推波场稳定,计算效率较高.
关键词: 一步外推法      耦合方程组      多倍角公式     
One-step extrapolation method based on the multiple-angle formula
KE Xuan1,2, SHI Ying1,2,3,4, WANG YinFeng5,6     
1. School of Earth Sciences, Northeast Petroleum University, Daqing 163318, China;
2. Key Laboratory of Oil and Gas Reservoir Formation Mechanism and Resource Evaluation in Heilongjiang Province, Daqing 163318, China;
3. Institute of Unconventional Oil and Gas, Northeast Petroleum University, Daqing 163318, China;
4. Key Laboratory of Continental Shale Hydrocarbon Accumulation and Efficient Development, Ministry of Education, Northeast Petroleum University, Daqing 163318, China;
5. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130021, China;
6. School of Mathematics and Statistics, Northeast Petroleum University, Daqing 163318, China
Abstract: We propose a coupled-equations solution for the one-step extrapolation method to solve the acoustic wave equation. We use the Euler formula to transform the one-step extrapolation equation for a complex wavefield into coupled equations of two real wavefields. Then this multiple-angle formula is applied to more accurate and stable approximations of the harmonic function operators of pseudodifferential operators. The polynomial Taylor expansions are adopted in the coupled equations and then the Fourier method is used to solve this operator and we derive a kind of coupled equations of the one-step extrapolation method. The method workflow requires fewer iterations of the Fourier transform than the existing one-step wave extrapolation matrix method. Besides, we also derive the stability conditions, which could provide a quantitative basis for choosing the multiple-angle parameters and the orders of polynomials to ensure the proposed algorithm stable. Numerical experiments for a two-dimensional constant-velocity model and complex structural model show that our method can keep stable of the extrapolation of the wavefield with a large time step and is more effective than the one-step wave extrapolation matrix method.
Keywords: One-step extrapolation    Coupled equations    Multiple-angle formula    
0 引言

近年来,混合域空间波数算子被用于地震波场递推,Fowler等(2010)Du等(2014)将这些方法称为递归积分时间外推(Recursive Integral Time Extrapolation,RITE)法,此类方法能在大时间步长情况下进行地震波场模拟,且不受数值频散干扰.根据对波动方程求解策略的不同,RITE法主要分为两类:

一类方法是基于多项式展开,Tal-Ezer(1986)Tal-Ezer等(1987)提出利用Chebyshev多项式求解声波波动方程,在此基础上,Kosloff等(1989)提出了快速展开(Rapid-Expansion,RE)法,后被Pestana和Stoffa(2009, 2010)、Pestana和Revelo(2017)应用至逆时偏移方法中;Revelo和Pestana(2019)将RE法推广应用于地震波场的上下行波分离、各向异性介质的逆时偏移成像及噪声压制,实现了该方法在非均匀介质成像领域中的应用;Araujo和Pestana(2020)在RE法中引入完全匹配层(Perfectly Matched Layer,PML)吸收边界条件,显著减少了边界反射干扰;Araujo同年提出使用多倍角公式以及二、四阶泰勒展开式近似余弦算子,提升了RE法的计算效率和精度(Araujo and Pestana, 2019).理论上,RE法可以看做是Lax-Wendroff方法(Lax and Wendroff, 1964Chen,2011)的一种,即利用等效的空间导数代替时间导数,计算精度和稳定性均优于传统方法中对二阶时间偏导项的直接近似策略.基于多项式展开方法的另一分支,即Zhang和Zhang (2009)提出的一步波场延拓方法,简称一步法(One-Step Wave Extrapolation,OSE),通过构造解析波场,将二阶声波方程转化为单程波方程形式进行求解,进而通过多项式近似实现了地震波场的稳定外推,该方法同样不受数值频散和稳定性条件影响.在此基础上,Revelo和Pestana(2016)提出了OSE矩阵(OSE Matrix,OSEM)解法,通过引入Jacobi-Anger展开式和Chebyshev多项式近似指数算子,求解Zhang(2009)的OSE法方程;该方法后被应用至上下行波分离(宋利伟等,2018)、逆时偏移的存储策略(柯璇和石颖,2017)和成像条件(Revelo et al.,2016)等领域.

另一类方法则是基于低秩分解理论,通过使用低秩分解方法近似空间-波数域地震波场外推算子,此类方法优势在于能够灵活折衷于计算效率和计算精度,实现波场的高效高精度模拟(Fomel et al.,2013).Sun等(2016a)将该方法应用于OSE逆时偏移方法中,以期实现稳定的波场递推和成像算法的同时,降低计算成本.近年来,低秩分解法同样广泛应用于黏声介质、各向异性介质中的地震波场数值模拟和成像算法中(Sun et al.,20152016b黄金强和李振春, 2017, 2019).

本文从声波方程的OSE解法出发,利用欧拉公式将OSE法的复数方程转化为关于虚实波场的耦合方程组,结合多倍角公式将方程组中的简谐函数算子(正弦函数算子、余弦函数算子)转换为关于简谐函数算子的多项式形式;分别采用三、四阶泰勒展开对简谐函数算子多项式进行逐项近似.相比于直接对简谐函数算子应用泰勒展开,能够获得更高的计算精度.针对本文方法进行了稳定性分析,为保证地震波场稳定外推的参数选取提供参考.理论分析和模型测试表明,本文方法计算精度高,且计算效率优于OSEM方法,更适合扩展应用至地震数据正演(Liao et al.,2009段沛然等,2019崔晓娜等,2020刘立彬等,2020)、成像(Shi and Wang, 2016Chen and Sacchi, 2017Ke et al.,2018Zhang and Shi, 2019Shi et al.,2019Liao et al.,2017Li et al.,2020刘伟等,2020赵超峰等,2020Zhang et al.,2020a2021)、反演(郭雪豹等,2016Yang et al.,2015; Zhang et al.,2020b)等方法中.

1 方法原理 1.1 OSE法的耦合方程组解法

二维常密度声波波动方程可表示为:

(1)

式中,P(x, t)为压力波场,▽2代表拉普拉斯算子,v(x)为介质速度场,其中x=(x, z),xz为二维空间笛卡尔坐标,t为时间.

根据Zhang和Zhang(2009)提出的OSE法,构建复数波场

(2)

式中,Q(x, t)是波场P(x, t)在时间方向的希尔伯特变换,i是虚数单位.方程(1)可重新整理为:

(3)

式中,.对于变速介质,方程(3)的解可表示为(Zhang and Zhang, 2009):

(4)

引入欧拉公式,式(4)可整理为:

(5)

分离(5)式中的实部虚部,可以得到两个耦合方程,为便于表述,省略波场中的空间坐标并将其整理为方程组形式:

(6)

然而,式(6)中的简谐函数算子无法直接进行数值求解,因此,后文将阐述如何利用多倍角公式和泰勒展开式对其进行近似求解.

1.2 简谐函数的多倍角公式和泰勒展开

Bromwich和MacRobert(1991)给出了正弦函数多倍角表达式:

(7)

为简化推导,本文只考虑n为奇数的情况,以保证(7)式中仅包含sin函数.因此,令n=2r+1,其中,r=1, 2, 3…,同时引入Ψr, j,且令:

(8)

经计算,系数Ψr, j如附录表A1所示,式(7)可整理为:

(9)

余弦函数多倍角表达式为(Bromwich and, 1991Araujo and Pestana, 2019b):

(10)

式中,为向下取整运算,系数.如果引入Γn, j,且令:

(11)

经计算,系数Γn, j如附录表A2所示,式(10)可整理为:

(12)

式中,"\"为取余计算.

将公式(9)中的x替换为x/(2r+1),公式(12)中的x替换为x/n,整理可得:

(13)

(14)

在阶数固定的情况下,泰勒展开式对小角度的简谐函数的近似效果更好.如图 1所示,当弧度x>1.5时,采用三、五阶泰勒展开式近似sin(x)函数(分别见图 1a中圆形、正三角形符号线),采用二、四阶泰勒展开式近似cos(x)函数(分别见图 1c中圆形、正三角形符号线)所取得的近似精度显著下降;而采用多倍角公式将sin(x)和cos(x)函数分别展开为关于sin[x/(2r+1)]和cos(x/n)的多项式形式后,再进行泰勒展开式近似则能取得较好效果,如图 1ac中标记为"T(r=1~5)"和"T(n=3~9)"的其他图形符号曲线所示.根据图 1bd所示的误差曲线也可看出,"多倍角+泰勒展开式"的近似方案能够取得更好的近似效果,且多倍角参数越大,近似效果越好.

图 1 多倍角公式和泰勒展开对简谐函数的近似结果 (a) 应用泰勒展开式和公式(13)近似的正弦函数;(b) 正弦函数近似的误差绝对值;(c) 应用泰勒展开式和公式(14)近似的余弦函数;(d) 余弦函数近似的误差绝对值. Fig. 1 Approximation results of harmonic function obtained by combining multiple-angle formula and Taylor expansion (a) Analytical sine function and approximations using the Taylor series and polynomial expression for equation (13); (b) Absolute error for approximations of the sine function; (c) Analytical cosine function and approximations using the Taylor series and the polynomial expression for equation (14); (d) Absolute error for approximations of the cosine function.

因此,利用式(13)、(14)即可将式(6)中简谐函数算子展开,进而对展开所得多项式再逐项进行泰勒展开,能够获得更高的近似精度,下节将对此进行详细论述.

1.3 OSE法的耦合方程组解法

为方便表示,引入算子Ir=sin[LΔt/(2r+1)]和On=cos(LΔt/n),并将式(13)和(14)代入式(6),可得:

(15)

方程组(15)是本文方法的基本框架,其中,IrOn则可应用诸如泰勒展开式、Chebyshev多项式、Hermite展开式或Legendre展开式等多项式形式近似求得.本文分别采用三、四阶泰勒展开式近似IrOn,则有:

(16)

(17)

式中的拟微分算子L,很难在空间域直接求解,需在波数域进行计算,即:

(18)

式中,FTIFT分别代表正、反傅里叶变换.将式(16)和(17)代入方程组(15)中,并将其写为波数域形式,可得:

(19)

式中,是波数矢量k=(kx, kz)的模,代表波场P(t)和Q(t)的波数域形式.根据式(19)可知,随着方程组中参数j的每一次递增,泰勒展开项中每个关于的指数项的计算均需一次正、反傅里叶变换.根据式(18)不难得知,对于不同幂次数的指数项的计算,主要区别在于所乘的波数补偿项次数的不同,再经反傅里叶变换即可获得最终结果.因此,理论上只要将方程中的泰勒展开项整理为适当的形式,即可避免冗余的正傅里叶变换(Zhang and Zhang, 2009).

为此,本文将方程组(19) 中的泰勒展开项进行展开并重新整理为关于的显式多项式的形式,经整理可得:

(20)

式中,SC分别为整理后的多项式系数.上标代表泰勒展开的阶数,下标rn分别表示正、余弦算子的多倍角参数,m为多项式的项数序号,经计算,本文给出了关于Sr, m3Cn, m4的系数表,详见附录表A3、A4.

理论上,方程组(20)完全等效于方程组(19).然而,由于方程组(20)中,经泰勒展开的正弦、余弦算子项被重新整理为关于的多项式形式,在程序实现时,针对的指数项的计算则只需要进行一次正傅里叶变换和若干次反傅里叶变换即可(具体次数等于项中的多项式数量),该方法能够大幅度降低计算量,提升计算效率.

1.4 稳定性分析及算法实现流程

本节分析了所提出方法的稳定性条件,应用平面波方程,并将其代入方程组(20),得到如式(21)所示的稳定性条件:

(21)

具体推导过程详见附录B.随着vkΔt的增大,算法稳定的条件是式(21)中不等号左侧多项式的值应分布于区间[-1, 1]内,根据式(21)所绘曲线如图 2所示.

图 2 据式(21)所绘稳定性曲线 (a) 采用不同r值和相应系数Sr, m3所绘关于的曲线;(b) 采用不同n值和相应系数Cn, m4所绘关于的曲线. Fig. 2 Stability curves drawn from equation (21) (a) Curves calculated using with different values of r and the corresponding coefficient Sr, m3; (b) Curves calculated using with different values of n and the corresponding coefficient Cn, m4.

显然,图 2中,不同的多倍角参数rn会导致曲线在区间[-1, 1]内延伸的长度不同,且延伸长度与rn呈正相关,表明rn的值越大,算法的稳定性越强.但需要注意的是,采用较大的多倍角参数能够提升算法的稳定性,但会随之带来更高的计算成本.根据式(21),可计算得出选取不同的多倍角参数rn后,算法能容许的最大vkΔt值,如表 1表 2所示.

表 1 方程组(20)中正弦函数算子的稳定性条件 Table 1 Stability conditions for the coupled-equations solution with one-step extrapolation for the sine function operator in equations (20)
表 2 方程组(20)中余弦函数算子的稳定性条件 Table 2 Stability conditions for the coupled-equations solution with one-step extrapolation for the cosine function operator in equations (20)

实际应用中,在确定了模型参数的情况下,即可计算出最大速度vmax和最大波数kNyquist=π,进而根据表 1表 2选取多倍角参数rn即可保证算法的稳定.

基于方程组(20)的波场延拓方法实现如流程1所示.

表 流程1 基于方程(20)进行波场延拓的实现流程
2 测试分析 2.1 计算量分析

由于RITE波场模拟方法中的拟微分算子需在波数域进行求解,因此需要进行多次傅里叶变换计算.本文方法中:①采用欧拉公式对方程(4)中的指数项进行展开,提出了基于多倍角公式OSE法的方程组解法,可对不同传播时间的实波场P(t)和虚波场Q(t)及其拟微分算子实现解耦计算,而这在OSEM法中是无法实现的(Revelo and Pestana, 2016);②在此基础上,如前节所述,通过将方程组(19)转化为方程组(20)的形式,使得本文方法在进行编程实现时,能够大幅减少正傅里叶变换次数,提升算法的执行效率.

为便于分析本文方法对一步法计算效率的提升,在满足稳定性条件的情况下,将本文方法与OSEM法(Revelo and Pestana, 2016)进行对比.通过比较二者所需的傅里叶变换次数来评估算法的计算成本,结果如表 3所示.

表 3 本文方法和OSEM法所需的傅里叶变换次数 Table 3 Numbers of Fourier transforms needed by the proposed method and the OSEM method

可以看出,随着vmaxkNyquistΔt的增加,两种算法均需引入更多次的傅里叶变换以维持算法稳定,当vmaxkNyquistΔt较小时,二者所需的傅里叶变换次数相当,随着vmaxkNyquistΔt的增大,本文提出的算法的计算量要明显小于OSEM法.

2.2 匀速介质模型测试

采用3000 m·s-1的匀速介质进行测试,模型网格尺寸为201×201.水平和垂向空间网格间距均为20 m,采用主频8 Hz的雷克子波作为震源子波,布置于(xz)=(2 km,2 km)处.截取400 ms时波场快照,如图 3所示,当多倍角参数为r=4、n=5,时间步长分别选取1 ms、2 ms、4 ms和8 ms时,根据所得波场快照可知,即使给定较大时间步长,本文方法仍能保持稳定的波场外推,且所得结果与应用小时间步长参数时无明显差异.

图 3 基于方程组(20),多倍角参数选取为r=4, n=5,计算所得波场快照(a) Δt=1 ms;(b) Δt=2 ms;(c) Δt=4 ms;(d) Δt=8 ms. Fig. 3 Snapshots of calculation based on equation (20) with r=4, n=5 and different time steps

为进一步分析本文方法的计算精度和稳定性,采用相同的模型参数,并在(x, z)=(2 km,1.6 km) 位置处设定检波器用于记录采用不同参数时所得的单道地震记录.由于时间步长越小,波场模拟的数值结果越接近解析解,因此在本次测试中,选择时间步长为1 ms时,OSEM法所得结果作为参考解,如图 4ad中的方形符号曲线所示.

图 4 基于方程组(20),采用不同时间步长、多倍角参数组合所得地震记录对比(a) Δt=1 ms;(b) Δt=2 ms;(c) Δt=4 ms;(d) Δt=8 ms. Fig. 4 Comparison of seismic records in the constant-velocity model based on equation (20) using different time steps

时间步长分别选取为1 ms、2 ms、4 ms和8 ms,多倍角参数设置分别为r=1,2,3,4和n=2,3,4,5的组合时,所得结果如图 4所示,采用不同参数的情况下,本文方法均能实现波场的稳定模拟,且各组合方案所得结果均与参考解无明显差异.

为进一步验证本文方法的适用性,模型参数保持不变,分别采用一、二阶泰勒展开近似公式(15)中的正、余弦算子,具体推导过程不再赘述,多倍角参数选取为r=4和n=8,分别选取了1 ms、2 ms、4 ms和8 ms的时间步长参数进行波场模拟测试.所得结果如图 5ad所示,根据波场快照可以发现,即使一定程度上降低泰勒展开阶数,采用不同时间步长参数的情况下,本文方法均能保持稳定的波场外推,所得结果间并无明显差异.

图 5 采用一、二阶泰勒展开近似,多倍角参数选取为r=4, n=8,计算所得波场快照(a) Δt=1 ms;(b) Δt=2 ms;(c) Δt=4 ms;(d) Δt=8 ms. Fig. 5 Snapshots from calculation using first- and second-order Taylor expansion approximations with r=4, n=8 and different time steps

在(x, z)=(2 km,1.6 km)位置处设定检波器记录采用不同参数时所得的单道地震记录,测试结果如图 6ad所示,当时间步长为1 ms、2 ms和4 ms时,本文选取的几组多倍角参数组合方法均能实现波场的稳定模拟,且各组合方案所得结果均与参考解无明显差异,如图 6ac所示.当时间步长取8 ms时,较小的多倍角参数无法满足稳定性条件,因此仅展示r=3、n=6和r=4、n=8时所得单道地震记录,如图 6d所示,所得结果与参考解整体匹配较好,局部放大后,虽然波场模拟结果和参考解存在稍许误差,但仍在可接受范围内.

图 6 采用一、二阶泰勒展开近似,不同时间步长、多倍角参数组合所得地震记录对比(a) Δt=1 ms;(b) Δt=2 ms;(c) Δt=4 ms;(d) Δt=8 ms. Fig. 6 Comparison of seismic records in the constant-velocity model using first- and second-order Taylor expansion approximations and different time steps
2.3 复杂构造模型测试

选用Marmousi模型,如图 7a所示,对本文方法和OSEM法进行波场模拟的对比测试,以比较分析两种方法的计算量.速度模型网格尺寸为水平方向737×垂直方向240,水平和垂向空间网格间距均为20 m.在(x, z)=(7 km,2 km)位置处载入主频为10 Hz的雷克子波作为震源.时间步长为2 ms,计算得出vmaxkNyquistΔt=2.44.因此,为满足稳定性条件,本文方法的参数选取为r=0、n=1;OSEM法的参数选取为M=3(Revelo and Pestana, 2016).两种方法的结果几近相同,均能稳定精确地模拟地震波场,如图 7b—c所示.在(x, z)=(7 km,1.6 km)位置处设定检波器记录所得的单道地震记录,如图 7d所示,本文方法所得结果与参考解匹配一致.然而,通过对本次测试所需计算量和实测计算耗时对比,如表 4所示,可知在每个时间步长内的计算中,本文方法需要更少的傅里叶变换计算,具有更高的计算效率;通过归一化整体计算耗时进行对比,本文方法的计算耗时也明显更小,且耗时比例与各方法所需傅里叶变换次数比例相近,也侧面印证了傅里叶变换计算贡献了算法的主要耗时.

图 7 Marmousi模型波场模拟测试 Marmousi速度模型;(b) 基于本文方法,多倍角参数选取为r=0, n=1所得波场快照;(c) 基于OSEM法,参数M=1所得波场快照;(d) 地震记录对比. Fig. 7 Wavefield simulation tests on Marmousi model (a) The Marmousi velocity model; (b) Based on proposed method with r=0, n=1; (c) Based on OSEM method with M=3; (d) Comparison of seismic.
表 4 基于Marmousi模型的本文方法和OSEM法的算法参数、傅里叶变换次数、计算耗时对比 Table 4 Test parameters, numbers of Fourier transforms needed, and normalized time for proposed method, and the OSEM method for the wavefield extrapolation test with the Marmousi model
3 结论与探讨

(1) 本文首先以OSE法为理论基础,提出了一种方程组解法并采用多倍角公式和泰勒展开进行求解.相比于OSEM法,本文方法能够显著减少算法所需傅里叶变换次数,提升算法执行效率;推导了所提方法的稳定性条件,以量化本文所提出的波场外推算法的稳定执行标准.为提升本文方法的可复制性,给出了伪代码流程;多个模型测试表明,本文方法能够适应大时间步长条件下的波场稳定模拟,对简单、复杂构造介质均具有良好的适应性.

(2) 本文所提方程组解法和多倍角公式的联合可作为一种理论框架,框架内可应用其他多项式进行简谐函数算子的近似求解,而不仅限于泰勒展开式;另外,针对方程组(20)中多项式系数可进行优化,以提升算法的稳定性,但相关内容本文尚未涉及,有待进一步研究.

附录A
表 A1 式(9)中系数Ψr, j Table A1 Coefficients Ψr, j in equation (9)
表 A2 式(12)中系数Γn, j Table A2 Coefficients Γn, j in equation (12)
表 A3 方程组(20)中的系数Sr, m3m=0~13 Table A3 Coefficients Sr, m3, m=0~13 in equations (20)
表 A4 方程组(20)中的系数Cn, m4m=0~10 Table A4 Coefficients Cn, m4, m=0~10 in equations (20)
附录B稳定性条件推导

均匀介质情况下,将方程组(6)转换至波数时间域形式:

(B1)

式中,分别为波数域的地震波场及其时间方向的希尔伯特变换形式(省略波数、时间坐标),k=是波数矢量k=(kx, kz)的范数,v为速度值.

分别将的平面波形式:

(B2)

(B3)

代入方程组(B1),其中i是虚数单位,能够得到:

(B4)

化简方程(B4),能够获得:

(B5)

采用欧拉公式展开(B5)中的e指数项,可得:

(B6)

分别令(B6)中等号左右两侧实部虚部对应相等,可得:

(B7)

对于任意正弦、余弦函数,其函数值应在闭区间[-1, 1]内,即:

(B8)

应用三阶、四阶泰勒展开近似sin(vkΔt)和cos(vkΔt)项,结合多倍角公式,经过推导可获得稳定性条件为:

(B9)

式中,Sr, m3Cn, m4详见附录A.针对不同的模型参数vkΔt,为满足稳定性条件,通过数值计算,能够确定不等式组(B9)中多项式项数rn,详见表B1和B2.

表 B1 方程组(20)满足稳定性条件的最大r Table B1 Maximum value of r for equations (20) meeting stability conditions
表 B2 方程组(20)满足稳定性条件的最大n Table B2 Maximum value of n for equations (20) meeting stability conditions
References
Arau jo E S, Pestana R C. 2019. Time-stepping wave-equation solution for seismic modeling using a multiple-angle formula and the Taylor expansion. Geophysics, 84(4): T299-T311. DOI:10.1190/GEO2018-0463.1
Araujo E S, Pestana R C. 2020. Perfectly matched layer boundary conditions for the second-order acoustic wave equation solved by the rapid expansion method. Geophysical Prospecting, 68(2): 572-590. DOI:10.1111/1365-2478.12868
Bromwich T J I, MacRobert T M. 1991. An Introduction to the Theory of Infinite Series. 3rd ed. New York: Chelsea: 202-207.
Chen J B. 2011. A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation. Geophysics, 76(2): T37-T42. DOI:10.1190/1.3554626
Chen K, Sacchi M D. 2017. Elastic least-squares reverse time migration via linearized elastic full waveform inversion with pseudo-Hessian preconditioning. //87th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts: 4364-4369.
Cui X N, Liu H. 2020. Fast multipole expansion of Helmholtz equation integral solution. Progress in Geophysics (in Chinese), 35(5): 1751-1758. DOI:10.6038/pg2020DD0214
Du X, Fowler P J, Fletcher R P. 2014. Recursive integral time-extrapolation methods for waves: A comparative review. Geophysics, 79(1): T9-T26. DOI:10.1190/GEO2013-0115.1
Duan P R, Li Q Y, Zhao Z Q, et al. 2019. High-order finite-difference method based on equivalent staggered grid scheme for scalar wavefield simulation. Progress in Geophysics (in Chinese), 34(3): 1032-1040. DOI:10.6038/pg2019CC0064
Fomel S, Ying L X, Song X L. 2013. Seismic wave extrapolation using lowrank symbol approximatio. Geophysical Prospecting, 61(3): 526-536. DOI:10.1111/j.1365-2478.2012.01064.x
Fowler P J, Du X, Fletcher R P. 2010. Recursive integral time extrapolation methods for scalar waves. //80th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts: 3210-3215.
Guo X B, Liu H, Shi Y. 2016. Time domain full waveform inversion based on frequency attenuation. Chinese Journal of Geophysics (in Chinese), 59(10): 3777-3787. DOI:10.6038/cjg20161022
Huang J Q, Li Z C. 2017. Modeling and reverse time migration of pure quasi-P-waves in complex TI media with a low-rank decomposition. Chinese Journal of Geophysics (in Chinese), 60(2): 704-721. DOI:10.6038/cjg20170223
Huang J Q, Li Z C. 2019. Least-squares reverse time migration of pure qP-wave pre-stack plane-waves based on low-rank finite difference for anisotropic media. Chinese Journal of Geophysics (in Chinese), 62(8): 3106-3129. DOI:10.6038/cjg2019M0188
Ke X, Shi Y. 2017. Forward simulation and reverse time migration imaging based on one step wavefield extrapolation. Chinese Journal of Geophysics (in Chinese), 60(11): 4468-4479. DOI:10.6038/cjg20171131
Ke X, Shi Y, Wang W H. 2018. An efficient wavefield simulation and reconstruction method for least-squares reverse time migration. Journal of Seismic Exploration, 27(2): 183-200.
Kosloff D, Filho A Q, Tessmer E, et al. 1989. Numerical solution of the acoustic and elastic wave equation by new rapid expansion method. Geophysical Prospecting, 37(4): 383-394. DOI:10.1111/j.1365-2478.1989.tb02212.x
Lax P D, Wendroff B. 1964. Difference schemes for hyperbolic equations with high order of accuracy. Communications on Pure and Applied Mathematics, 17(3): 381-398. DOI:10.1002/cpa.3160170311
Li C, Liu J X, Liao J P, et al. 2020. 2D high-resolution crosswell seismic traveltime tomography. Journal of Environmental and Engineering Geophysics, 25(1): 47-53. DOI:10.2113/JEEG19-003
Liao J P, Guo Z W, Liu H X, et al. 2017. Application of frequency-dependent traveltime tomography to 2D crosswell seismic field Data. Journal of Environmental and Engineering Geophysics, 22(4): 421-426. DOI:10.2113/JEEG22.4.421
Liao J P, Wang H Z, Ma Z T. 2009. 2-D elastic wave modeling with frequency-space 25-point finite-difference operators. Applied Geophysics, 6(3): 259-266. DOI:10.1007/s11770-009-0029-7
Liu L B, Duan P R, Zhang Y Y, et al. 2020. Overview of mesh-free method of seismic forward numerical simulation. Progress in Geophysics (in Chinese), 35(5): 1815-1825. DOI:10.6038/pg2020DD0117
Liu W, Li H, Zhang W. 2020. Comparative study on elastic reverse time migration based on two methods of wavesfield separation. Progress in Geophysics (in Chinese), 35(3): 1000-1009. DOI:10.6038/pg2020DD0137
Pestana R C, Stoffa P L. 2009. Rapid expansion method (s) for time-stepping in reverse time migration (RTM). //79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts: 2819-2823. DOI:10.1190/1.32255434
Pestana R C, Stoffa P L. 2010. Time evolution of the wave equation using rapid expansion method. Geophysics, 75(4): T121-T131. DOI:10.1190/1.3449091
Pestana R P, Revelo D. 2017. An improved method to calculate the analytical wavefield for causal imaging condition. //89th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts: 9-17.
Revelo D, Pestana R. 2019. De-primary TTI-RTM using the P-pure analytical wavefield. //89th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts: 1-5. DOI:10.3997/2214-4609.201900838
Revelo D, Pestana R C, Gomez L. 2016. Reverse Time Migration (RTM) using analytical wavefield and causal imaging condition. 86th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts: 1-5. DOI:10.3997/2214-4609.201601372
Revelo D E, Pestana R C. 2016. One-step wave extrapolation matrix method for reverse time migration. Geophysics, 81(5): S359-S366. DOI:10.1190/GEO2015-0555.1
Shi Y, Wang Y H. 2016. Reverse time migration of 3D vertical seismic profile data. Geophysics, 81(1): S31-S38. DOI:10.1190/GEO2015-0277.1
Shi Y, Zhang W, Wang Y H. 2019. Seismic elastic RTM with vector-wavefield decomposition. Journal of Geophysics and Engineering, 16(3): 509-524. DOI:10.1093/jge/gxz023
Song L W, Shi Y, Ke X. 2018. Reverse time migration imaging based on analytic wavefield separation. Progress in Geophysics (in Chinese), 33(4): 1573-1578. DOI:10.6038/pg2018BB0445
Sun J Z, Fomel S, Ying L X. 2016a. Low-rank one-step wave extrapolation for reverse time migration. Geophysics, 81(1): S39-S54. DOI:10.1190/GEO2015-0183.1
Sun J Z, Fomel S, Zhu T Y, et al. 2016b. Q-compensated least-squares reverse time migration using low-rank one-step wave extrapolation. Geophysics, 81(4): S271-S279. DOI:10.1190/GEO2015-0520.1
Sun J Z, Zhu T Y, Fomel S. 2015. Viscoacoustic modeling and imaging using low-rank approximation. Geophysics, 80(5): A103-A108. DOI:10.1190/geo2015-0083.1
Tal-Ezer H. 1986. Spectral methods in time for hyperbolic equations. Siam Journal on Numerical Analysis, 23(1): 11-26. DOI:10.1137/0723002
Tal-Ezer H, Kosloff D, Koren Z. 1987. An accurate scheme for seismic forward modelling. Geophysical Prospecting, 35(5): 479-490. DOI:10.1111/j.1365-2478.1987.tb00830.x
Yang P L, Gao J H, Wang B L. 2015. A graphics processing unit implementation of time-domain full-waveform inversion. Geophysics, 80(3): F31-F39. DOI:10.1190/geo2014-0283.1
Zhang W, Gao J H, Gao Z Q, et al. 2020a. 2D and 3D amplitude-preserving elastic reverse time migration based on the vector-decomposed P- and S-wave records. Geophysical Prospecting, 68(9): 2712-2737. DOI:10.1111/1365-2478.13023
Zhang W, Gao J H, Gao Z Q, et al. 2020b. Adjoint-driven deep-learning seismic full-waveform inversion. IEEE Transactions on Geoscience and Remote Sensing: 1-20. DOI:10.1109/TGRS.2020.3044065
Zhang W, Gao J H, Li F P, et al. 2021. Correlative least-squares reverse time migration in viscoelastic media. Journal of Applied Geophysics, 185: 104256. DOI:10.1016/j.jappgeo.2021.104256
Zhang W, Shi Y. 2019. Imaging conditions for elastic reverse time migration. Geophysics, 84(2): S95-S111. DOI:10.1190/geo2018-0197.1
Zhang Y, Zhang G Q. 2009. One-step extrapolation method for reverse time migration. Geophysics, 74(4): A29-A33. DOI:10.1190/1.3123476
Zhao C F, Tian J T, Zhang W, et al. 2020. Application research of dual-well Walkaway VSP technology in Shuguang high-steep structure of Liaohe depression. Progress in Geophysics (in Chinese), 35(1): 265-271. DOI:10.6038/pg2020DD0002
崔晓娜, 刘洪. 2020. Helmholtz方程积分解的快速多极子展开. 地球物理学进展, 35(5): 1751-1758.
段沛然, 李青阳, 赵志强, 等. 2019. 等效交错网格高阶有限差分法标量波波场模拟. 地球物理学进展, 34(3): 1032-1040.
郭雪豹, 刘洪, 石颖. 2016. 基于频域衰减的时域全波形反演. 地球物理学报, 59(10): 3777-3787. DOI:10.6038/cjg20161022
黄金强, 李振春. 2017. 基于Low-rank分解的复杂TI介质纯qP波正演模拟与逆时偏移. 地球物理学报, 60(2): 704-721. DOI:10.6038/cjg20170223
黄金强, 李振春. 2019. 各向异性介质Low-rank有限差分法纯qP波叠前平面波最小二乘逆时偏移. 地球物理学报, 62(8): 3106-3129. DOI:10.6038/cjg2019M0188
柯璇, 石颖. 2017. 基于一步法波场延拓的正演模拟和逆时偏移成像. 地球物理学报, 60(11): 4468-4479. DOI:10.6038/cjg20171131
刘立彬, 段沛然, 张云银, 等. 2020. 基于无网格的地震波场数值模拟方法综述. 地球物理学进展, 25(5): 1815-1825.
刘伟, 李慧, 张伟. 2020. 基于两种波场分离方法的弹性波逆时偏移对比研究. 地球物理学进展, 35(3): 1000-1009.
宋利伟, 石颖, 柯璇. 2018. 基于解析法波场分离的逆时偏移成像. 地球物理学进展, 33(4): 1573-1578. DOI:10.6038/pg2018BB0445
赵超峰, 田建涛, 张伟, 等. 2020. 双井Walkaway VSP技术在辽河坳陷曙光高陡构造地区的应用研究. 地球物理学进展, 35(1): 265-271. DOI:10.6038/pg2020DD0002