1 引言

地震数据正演，逆时偏移及全波形反演等方法均需要以高精度的地震波场数值模拟算法为基础.有限差分数值方法求解波动方程易于实现，内存需求少，因此在波场模拟和逆时偏移中得到广泛的应用(Hemon，1978；Mcmechan，1983；Villarreal et al., 1997；Matthew and Clive, 2006).但是，有限差分法中，空间网格和时间采样步长的选取受稳定性条件限制，而且差分计算的误差会导致地震波场中出现数值频散干扰.提高有限差分阶数(刘红伟等，2010)或采用通量校正传输(Tong and Ken, 1995)方法能够一定程度上压制频散，但需要引入额外的计算量；Zhang和Yao(2013)采用模拟退火法优化有限差分系数，无需增加额外计算量，即可降低频散影响.上述方法均能压制频散干扰，但无法从根本上消除频散，而且为保持差分算法的稳定性，通常无法选择较大的时间步长参数(Carcione et al., 2002；Moczo et al., 2007)，因此对实际数据进行偏移成像时，往往需要预先对数据进行插值处理，增加了输入的数据量.

选择稳定和精确的波动方程求解方法是实现高精度地震波场正演模拟，提高逆时偏移成像质量的关键.Tal-Ezer等(1986, 1987)提出利用时间积分法实现声波方程波场模拟，能够获得比有限差分法波场模拟更高的计算精度和算法稳定性；Zhang和Zhang(2009)从单程波动方程的解析表达式出发，提出了一步波场延拓方法，该方法在双程波动方程中引入复波场，将其重构为类似单程波动方程的形式，进而得出波场延拓计算的解析递推公式，该方法可实现大时间步长的波场延拓，算法稳定性较强；Pestana和Stoffa(2010)提出快速展开法，仅保留了时间积分法(Tal-Ezer，1986；Tal-Ezer et al., 1987)递推公式的实数项，避免了复数波场运算，该类方法也被称为“两步法”(Robert and Zhang, 2008)；Formel等(2013)进一步采用低秩法近似波动方程的递推公式，有效降低了波场模拟的计算量；Daniel和Reynam(2016)基于矩阵形式推导了一步波场延拓法(后文简称一步法)，采用Chebyshev多项式近似e指数项，时间步长的选取无稳定性条件限制，在逼近Nyquist采样率的情况下仍可实现波场的稳定传播；Du等(2014)将上述波场模拟方法(Tal-Ezer et al., 1986, 1987；Zhang and Zhang, 2009；Pestana and Stoffa, 2010；Formel et al., 2013)归类命名为“递归时间积分延拓法”(RITE，Recursive integral time-extrapolation)，指出该类方法均具有较高的稳定性，能够实现大时间步长条件下的波场模拟.近年来，“递归时间积分延拓法”的应用不再局限于各向同性介质，在复杂各向异性介质(Pestana et al., 2011；Pestana et al., 2012；Zhan et al., 2012；黄金强和李振春，2017)以及粘性介质(Sun et al., 2015)的正演、偏移、最小二乘偏移及全波形反演等方法中均取得了成功的应用(Sun et al., 2014；Santos and Pestana, 2015).

为降低逆时偏移方法的波场存储量，诸多学者展开了细致的研究.Feng等(2012)提出保存边界处单层波场的策略，通过变阶数有限差分实现波场逆推重建，但波场重建的精度有所降低；与之类似的有效边界存储方法(王保利等，2012)虽然能改善精度问题，实现波场的无损重建，但该方法也局限于在有限差分波场模拟方法中使用，并不适用于基于“递归时间积分延拓法”的逆时偏移计算.Symes(2007)提出的优化checkpoint方法，通过优选checkpoint降低波场存储量，但会引入较多的额外计算量，且实现过程较为复杂.Robert(2009)提出了随机边界条件，虽然不受波场模拟算法限制，而且无需存储临时波场，通过方程反推即可实现波场重建，但随机噪声的引入会影响偏移效果.综上所述，基于“递归时间积分延拓法”的逆时偏移方法的波场重建策略仍有待进一步研究探索.

本文从二维常密度声波方程出发，引入复波场，利用Chebyshev多项式和Jacobi-Anger展开式，推导了一步法波场模拟方程，相比于有限差分法，一步法更能适应大时间步长波场模拟，无数值频散干扰，算法稳定性强；随后本文将一步法地震波场模拟引入逆时偏移成像，并提出了一种适用于一步法逆时偏移方法的双缓冲区存储策略，在不增加计算量的前提下，二级缓冲区的引入大幅提升了内存重复利用率，有效降低了波场存储量.采用双缓冲区存储策略时，参数的选取对波场的内存占用量影响明显，因此本文进一步推导出令波场存储量占用最小的参数选取方案.模型试算结果表明，文中方法波场模拟精度高，逆时偏移成像质量好，波场存储成本低.

2 方法原理 2.1 矩阵表示的一步法波场延拓方法

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

(1)

式中，P=P(x, z, t)是压力场，v=v(x, z, t)是速度场.

在(1)式中引入复波场(Zhang and Zhang, 2009)

(2)

式中，Q(x, z, t)是P(x, z, t)在时间方向上的希尔伯特变换，i是虚数单位.

定义伪差分算子则复波场对如下一阶方程成立(Zhang and Zhang, 2009)：

(3)

(3) 式等价于双程波方程，可获得和常规声波波动方程相同的走时和振幅信息.

若速度v为常数，(3)式的解可表示为

(4)

式中，F和F^－1分别为正、反傅里叶变换.

若速度v变化，(3)式的解可表示为(Zhang and Zhang, 2009)

(5)

为了简化复数运算，把(2)式代入(3)式，可得如下的耦合方程组：

(6)

令则(6)式可用矩阵简化表示为：

(7)

方程(7)的解可表示为

(8)

为计算指数算子e^AΔt，引入如下Jacobi-Anger展开式

(9)

式中，J_n为n阶一类贝塞尔函数.

若令

(10)

(11)

则有

(12)

式中，R为A的最大本征值.在二维情况下，可以选取其中，v_max为最大速度，Δx和Δz分别为水平和垂直方向的空间网格间距.综上可知，e^AΔt可以表示为

(13)

由第一类Chebyshev多项式可知，若定义x=cos(θ)，T_n(x)=cos(nθ)，T₀(x)=1，T₁(x)=x=cos(θ)，则有如下递推公式成立：

(14)

若令=其中i为虚数单位，不难得出递推公式

(15)

则(13)式可进一步化简为

(16)

为了确保(16)式收敛，需保证M＞RΔt(Tal-Ezer，1986).将(16)式代入(8)式中并展开可得

(17)

(17)式即为一步法波场延拓递推公式.

对(17)式中Γ_n项进一步展开，可得

(18)

(19)

(20)

由于矩阵A中的在时域内无法计算，因此需在波数域内进行计算(Zhang and Zhang, 2009)得

(21)

式中，k_x和k_z为空间波数.为降低计算成本，本文采用GPU加速技术(刘国峰等，2009；李博等，2010；Liu et al., 2012；石颖等，2010；2013)实现上述算法，结合基于CUDA平台的CUFFT进行傅里叶变换计算，可大幅提升计算效率(郭雪豹等，2016).

一步法的算法稳定性要明显优于有限差分法，在时间采样率和空间网格间距的选取上具有更高的宽容度.理论上，时间采样间隔的选取仅受限于Nyquist采样定理(Daniel and Reynam, 2016)，即

(22)

式中，Δt为时间采样间隔，f_max为震源子波的最大频率.给定震源子波最大频率f_max后，关系式(Robert and Zhang, 2008；Du et al., 2014)

(23)

在满足Nyquist定理的情况下，Δt的最大值为代入(23)式则有

(24)

假设Δx=Δz，有

(25)

由于地下介质为变速介质，为保证算法的全局稳定性，当v=v_max为地下介质最大速度时，将(25)式代入(24)式可得

(26)

综上可知，(26)式即为一步法的稳定性条件.

2.2 双缓冲区存储策略

为节约存储空间，一些学者在利用有限差分数值方法进行逆时偏移计算时，提出存储正向传播的局部震源波场，再反向重建历史震源波场的方法(Feng，2012；王保利等，2012)，但是此类方法不适用于本文研究的基于一步法波场延拓的逆时偏移方法.为此，本文提出了一种双缓冲区(一级缓冲区和二级缓冲区)存储策略.

该方法首先定义一级缓冲区(Buffer-L1)和二级缓冲区(Buffer-L2)，如图 1所示，再沿时间方向等间隔定义N个checkpoint，一级缓冲区用于存储所有checkpoint处的正传震源波场(含吸收边界)，二级缓冲区用于存储某两个checkpoint间的正传震源波场(不含吸收边界).震源波场正传时，保存所有checkpoint处波场至一级缓冲区；当震源波场反传时，以预先载入一级缓冲区内的某checkpoint时刻波场为初始条件，将波场正传至下一个checkpoint时刻，期间保存两个checkpoint之间所有时刻的波场至二级缓冲区.以第N-2个checkpoint为例，当从第N-2个checkpoint向N-1个checkpoint正传震源波场时，保存第N-2个至N-1个checkpoint之间所有时刻的波场(不含吸收边界)至二级缓冲区，此时第N-2至N-1个checkpoint之间的波场覆盖更新了二级缓冲区内的原有波场(第N-1个checkpoint至第N个checkpoint之间的波场).检波点波场的反传与上述过程同时进行，并对二级缓冲区中的震源波场与反传的检波点波场进行零延迟互相关成像，以此类推，直至完成全部时刻的成像计算.

假设一级缓冲区需保存N个全波场(含吸收边界)数据，二级缓冲区需保存n个有效波场(不含吸收边界)数据，则N和n的选取直接影响缓冲区空间的占用，因此本文推导了使缓冲区存储占用最小的参数选取方案.为方便推导，存储量占用以数据网格点数为单位.

设二维模型的水平和垂直网格点数分别为N_x，N_z，吸收边界层数为B，则一个时刻的包含吸收边界的全波场数据的网格点数为S₁=(N_x+2B)(N_z+2B)，一个时刻的不包含吸收边界的有效波场数据的网格点数为S₂=N_xN_z.为满足震源波场能够递推计算的初始条件，需在每个checkpoint处保存P(t)和Q(t)的全波场信息，则一级缓冲区所需保存的数据的网格点数为C₁=2S₁N；二级缓冲区用于临时保存两个一级缓冲区点之间的有效波场，用于互相关成像，仅需保存P(t)的有效波场信息，则二级缓冲区所需保存的数据的网格点数为C₂=S₂n.设时间方向的采样点数为T，理想情况下假设T可被N整除，则满足T=Nn.

据上可知，所需存储的所有数据的网格点数为

(27)

为使总波场存储量占用最小，即选取N和n的值，满足f(N, n)最小.设g(N, n)=T－Nn作为限定条件，由拉格朗日乘子法可构建目标方程为

(28)

(28)式对自变量求偏导可得方程组：

(29)

求解(29)式可得

(30)

在方法实现过程中，N和n需取整数，因此当时，波场的存储量占用最小，〈〉代表四舍五入取整操作.将N和n代入(27)式，即可推算出波场的存储占用情况.

基于双缓冲区存储策略的逆时偏移算法如图 2所示.图中，it表示波场递推过程中时间方向的计数变量，i为读取二级缓冲区内的波场时的计数变量，N，n和T的含义同上文，虚线框标识部分即为根据图 1所示采用双缓冲区存储方法的波场重建过程.

3 模型测试 3.1 波场模拟测试

采用有限差分法模拟地震波场时，空间网格间距和时间步长的大小是影响波场模拟精度的两个重要因素.为了对比分析，下文分别在不同空间网格间距的情况下，测试有限差分法和一步法模拟地震波场的精度.也进一步测试了在不同时间步长情况下，一步法波场模拟的精度.

3.1.1 不同空间网格间距波场模拟

本文设计如图 3所示的水平层状模型，测试不同网格间距对有限差分法和一步法波场模拟精度的影响.模型沿水平和垂直方向均为3 km，速度为1700 m·s^-1的水层和速度为3000 m·s^-1的岩层的分界面位于垂向1.5 km处.激发震源采用主频为20 Hz的雷克子波，位于水平方向1.5 km，垂直方向1.2 km处.选择空间网格间距分别10，15，20 m进行测试，因此三种情况下，模型沿水平方向和垂直方向的网格点数分别为301×301，201×201，151×151，时间采样步长恒定为1 ms，分别采用空间十阶有限差分法和一步法波场模拟，0.6 s时刻的波场快照如图 4和5所示.观察可知，随着空间网格间距的增加，有限差分法模拟的波场频散现象逐渐严重；一步法模拟的波场均较为稳定，误差干扰较小.

3.1.2 不同时间步长波场模拟

本文也设计如图 6a所示的复杂字母模型，用于测试时间步长对一步法波场模拟精度的影响.假定字母模型沿水平和垂直方向的网格点数均为201，空间网格间距为15 m，激发震源位于水平方向1.5 km，垂直方向1.5 km处，分别选择1 ms，2 ms和5 ms的时间步长进行一步法波场模拟测试，在0.4 s时刻的波场快照分别如图 6b，6c和6d所示，观察可知，随着时间步长的增加，采用一步法模拟的地震波场仍具有良好的稳定性.图 7显示了沿图 6b，6c和6d中虚线抽取的地震数据.对比分析进一步说明，对采用不同时间步长的三个模拟结果，基于一步法的地震波场模拟的精度高、稳定性好，抽取的单道数据吻合较好，地震波走时、振幅和相位都能保持一致.

3.1.3 波场模拟耗时对比

虽然一步法稳定性较好，但相比于有限差分法，其算法略加复杂，因此计算成本也有所增加.本文分别采用一步法和有限差分法进行地震记录正演模拟的耗时测试.正演模拟采用的速度模型为如图 8所示的复杂字母模型，水平和垂直方向的网格点数均为201，空间网格间距为20 m.震源位于地表水平方向1.5 km，震源子波选取主频15 Hz的雷克子波，记录时长为2 s.一步法中参数M=7，有限差分法采用时间二阶，空间十阶差分.一步法时间间隔选取5 ms，有限差分法选取时间间隔选取1 ms，即一步法正演的时间采样点数为500，有限差分法正演的时间采样点数为2000，正演所得共炮点道集数据如图 9所示，具体参数及耗时对比如表 1所示.

由图 9可知，两种方法均能准确有效的进行地震数据正演模拟，由于一步法能采用更大的时间步长，因此所需模拟的时间点数相对更少，但由于其计算方法更加复杂，因此计算效率仍要低于有限差分法.

3.2 逆时偏移测试与存储量分析 3.2.1 逆时偏移测试

本文采用Marmousi2的截断模型测试一步法逆时偏移成像，图 10a和10b分别为精确速度模型和平滑速度模型.模型沿水平方向和垂直方向的网格点数分别为400和300，空间网格间距为15 m，时间采样间隔为4 ms，采样点数为1500，设计主频15 Hz的雷克子波作为激发震源，20个震源的炮间距为300 m，沿水平方向0~6000 m范围内均匀分布，每炮199个检波器接收，均匀对称分布于震源两侧，检波器间距为15 m，最大偏移距为1485 m.分别采用一步法和有限差分法进行逆时偏移计算，两种情况均采用了e指数吸收边界条件和互相关成像条件.为了满足稳定性条件，利用有限差分法进行逆时偏移成像时，本文采用sinc函数插值方法对炮集数据进行时间方向插值，插值重建后的时间步长为1 ms，时间方向采样点数为6000.图 10c和10d分别为一步法逆时偏移和有限差分法逆时偏移成像结果.对比可知，一步法波场模拟不存在频散现象，波形传播稳定无畸变，聚焦性更强，逆时偏移成像的同相轴更清晰，剖面具有更高的信噪比.有限差分法受频散现象干扰，逆时偏移剖面分辨率相对较低，且具有较多成像噪声.

上述测试选取了100层网格的偏移孔径和50层网格的边界层，则计算全波场的网格点数为水平方向700，垂直方向400；其中有效波场网格点数为水平方向600，垂直方向300.根据2.2节内容可知，对应的S₁为280000，S₂为180000，T为1500，进而可计算出n和N分别为68和22，代入(17)式可得双缓冲区存储策略的波场存储量约为93.68 MB，与实测数据相符.在不增加额外计算量的情况下，以空间十阶差分为例，采用有效边界存储策略(王保利等，2012)计算逆时偏移所需波场存储约为205.08 MB，如表 2所示.

观察可知，相比于节约存储的有效边界存储方法，基于双缓冲区存储方法的逆时偏移成像所需的波场存储量更小，主要是由于一步法波场模拟可采用大时间步长，从而减小时间采样点数，进而减少了波场重建数量.另外，双缓冲区存储策略中二级缓冲区的重复利用以及缓冲区参数的最优化选取也进一步降低了波场的存储需求.

3.2.2 存储量分析

为了定量表示本文方法在不同模型尺寸下的波场存储量占用情况，也测试了N和n随模型尺寸的变化情况，假设水平和垂直方向空间网格点相同，吸收边界层数为50，时间采样点数1000，根据(30)式可计算得到N和n随模型尺寸的变化曲线如图 11所示.由于本文方法可采用大时间步长进行波场模拟，时间采样点为1000时即可进行时长约4 s的波场模拟.因此由图 12可知，在二维空间情况下(假设横向和纵向网格点数相同，且均小于1000；每个网格点数据占用4bytes内存空间)，本文方法能够将波场存储量控制在100 MB以下.较小的存储需求使本文方法对硬件要求较低，另根据前文所述可知，算法对一级缓冲区的访问为非连续访问，因此，应用GPU加速技术时，一级缓冲区可设置在主机内存端，能够进一步降低显存存储负担，提升本文方法的适用性.

4 结语

(1) 本文基于一步法波场延拓方法开展波场模拟研究，并将其引入逆时偏移成像，提出了降低波场存储量的双缓冲区存储方法.该方法通过设置二级缓冲区以及选取最优参数，可在不增加逆时偏移计算量的情况下，大幅降低逆时偏移的波场存储量.相比于有限差分方法，一步法波场模拟稳定性强，受模型网格间距和时间步长影响小，无频散假象，波场模拟精度更高.但由于一步法波场延拓算法略加复杂，需要多次执行傅里叶变换，因此也带来了计算耗时的增加.

(2) 多个模型验证了文中方法是一种具有广阔发展前景的地震波场模拟方法，可成功应用于叠前逆时偏移成像，具有较大的工业应用潜力.此外，文中方法也可推广应用于最小二乘逆时偏移及全波形反演等方法中.

致谢

感谢审稿专家和编辑部的大力支持.