0 引言

叠前时间偏移是地震数据处理过程中的重要环节，其对偏移速度场的精度要求不高、且可以满足大部分探区的成像精度需求.相比于应用条件苛刻的有限差分偏移(Claerbout, 1985；王华忠等, 2009)和频率-波数域偏移(Stolt et al., 1987; Ekren and Ursin, 1999; Chapman, 2004)，Kirchhoff叠前时间偏移(Schneider, 1978；Geiger, 2002；邹振等，2010；初海红，2018；谢万学等，2018；岳玉波等，2018)实现方式灵活，可以很容易的处理起伏地表采集和非规则采样的地震数据，并且具备较高的计算效率和灵活性，因而成为时间偏移处理的首选技术手段.Kirchhoff偏移的核心是地震数据到成像剖面的映射累加运算，其计算量正比于叠前地震数据的道数，对当前动辄几十、几百TB甚至PB级的海量地震数据进行偏移成像处理，需要付出高昂的计算成本.

由于高密度地震数据在接收线(或炮线)方向的采样间隔很小，因此完全可以利用地表某点到成像点走时的一阶泰勒展开来近似计算临近点到成像点的走时，由此发展而来的射线束深度偏移方法(Hill, 2001; Gray, 2005; Gray and Bleistein, 2009; 岳玉波等, 2012, 2019a, b)在业界得到了迅速发展和广泛应用.该类方法计算效率接近于常规Kirchhoff叠前深度偏移，且具备处理多波至的能力、成像精度更高.射线束偏移同样可以应用于时间域成像，蔡杰雄(2010)提出了一种射线束时间偏移方法，目的是提高偏移的计算效率，然而作者只是利用中心道预先计算的射线参数简化了相邻道的走时计算过程，仍然需要将每道数据映射累加运算，并没有显著降低偏移的计算量.

本文提出了一种快速射线束叠前时间偏移方法(FBPSTM)，该方法基于地震波走时的一阶泰勒展开，利用倾斜叠加将地震数据分解为局部平面波，然后利用射线束中心处的地震波走时来完成局部平面波的成像映射累加运算.FBPSTM的实现过程分为以下三步：首先，根据给定的射线束中心间隔将炮集数据划分为一系列子集；接下来，利用局部倾斜叠加将数据子集分解为不同方向的平面波；最后，根据射线束中心处求取的双程地震波走时和射线参数拾取相应的局部平面波采样振幅累加到成像点上.同常规Kirchhoff叠前时间偏移(KPSTM)相比，FBPSTM不但保持了一致的成像精度，且由于仅在稀疏的射线束中心位置进行成像映射累加运算(计算量仅与选择的射线束中心间隔有关，和地震道数无关)，计算效率得到了大幅度提升，十分适用于宽方位、高密度地震数据的叠前时间偏移处理.

1 方法原理

本节以二维共炮域偏移为例，介绍FBPSTM的基本原理和实现过程，并同KPSTM进行计算量的对比分析.

1.1 KPSTM

如图 1所示，假设震源和接收点的位置分别为 x_s=(x_s, 0)和 x_r=(x_r, 0)，对应的地震记录为U( x_r, x_s, t)，地下成像点的位置为 x =(x, t₀)，那么地下成像点处的成像值R(x, t₀)可以用Kirchhoff积分偏移公式来表示，即：

(1)

其中，W(x_s, x, x_r)为偏移加权函数，为半导数滤波因子，用来保证成像相位的正确性，T为地震波双程走时，可以用双平方根公式表示为

(2)

其中，t_s和t_r分别为震源和接收点到成像点的单程走时，V_rms为成像点处的均方根速度.由式(1)不难看出，KPSTM的计算量正比于地震数据的总道数.

1.2 FBPSTM

当前地震勘探广泛采用宽方位、高密度的采集方式，相邻的接收点(或炮点)间距很小，因此完全可以使用地表某个接收点走时的一阶泰勒展开来近似相邻接收点的走时.如图 1所示，若已知射线束中心L =(L, 0)到成像点的走时为t_L，那么便可以利用式(3)来近似计算L附近接收点x_r到成像点的走时，即：

(3)

其中：

(4)

为射线束中心处射线参数的水平分量.此时，公式(1)可以表示为

(5)

式(5)中第三个累加项实际上是射线束中心附近空间窗|x_r－L| < ΔL内地震数据的倾斜叠加(Hill, 2001; Gray, 2005; Yue et al., 2019).式(3)的计算误差随接收点到射线束中心距离逐步增大，为减小该误差，在式(4)中引入高斯窗(Hill, 2001; Gray and Bleistein, 2009; Yue et al., 2019)，得到最终的FBPSTM成像公式为

(6)

其中，C为常数因子，用保证高斯窗的累加结果等于1，D(p_x, τ)为射线束中心附近满足|x_r－L| < 2ΔL的地震道的加窗局部倾斜叠加，公式为

(7)

1.3 FBPSTM计算效率分析

由式(6)和式(7)可知，FBPSTM的计算量主要分为两部分：一部分是式(6)所表示的成像累加过程，其计算量正比于射线束中心数，假设射线束中心间隔为接收点间隔的N倍，那么式(6)的计算量约为式(1)表示的KPSTM的1/N；另一部分是式(7)所表示的局部倾斜叠加.对于一个输入地震道，倾斜叠加的计算复杂度约为O(np_x×nτ)，KPSTM的计算复杂度为O(nx×nt₀)，由于偏移孔径内的成像道数nx要远大于射线参数采样数np_x，因此倾斜叠加的计算量同样要远小于KPSTM.此外，若地震数据采样规则，还可以利用傅立叶变换在频率波数域对式(7)进行更为高效的计算.上述分析结果将在接下来的试算测试中进行验证.

射线束中心间隔ΔL是FBPSTM的重要参数，选择较大的ΔL会提高计算效率，但成像精度会有所损失，选择较小的ΔL虽然可以保证成像精度，但同时也会降低计算效率.ΔL的选取同地震波的平均速度V_avg和主频f₀有关，本文使用了下述公式来求取射线束中心间隔，可以保证计算效率和成像精度的良好平衡，ΔL表达式为

(8)

1.4 实现过程

FBPSTM的计算过程可以大致概括为：

(1) 读入一炮地震数据，根据射线束中心间隔将炮记录划分为不同的子集；

(2) 对于不同的数据子集，根据给定的射线参数范围和间隔，计算式(7)得到对应不同射线参数的局部平面波D(p_x, τ)；

(3) 在不同的射线束中心位置，利用式(2)计算对应每个成像点的地震波双程走时T=t_s+t_L和水平射线参数p_L，然后提取D(p_L, T)对应的振幅值并根据式(6)累加到成像结果R(x, t₀)；

(4) 重复上述计算过程，直到所有炮记录计算完成.

2 模型及实际资料试算

分别使用二维模型和实际地震数据进行FBPSTM的测试，并同KPSTM进行成像精度和计算效率的对比.

2.1 二维洼陷模型

首先使用洼陷模型数据进行测试.图 2a为模型的深度域层速度场，模型网格为640×750，水平采样间距为15 m，垂直采样间距为4 m.使用有限差分法正演了120炮单炮记录，每炮120道，炮间距和道间距均为30 m.接下来以切除直达波后的第61炮记录(图 2b)为例，展示FBPSTM的实现过程.首先设定射线束中心间隔为210 m，射线参数采样为30，根据射线束中心间隔可以将炮记录划分为21个数据子集，然后对于每个数据子集应用式(7)所示的局部倾斜叠加求得对应不同射线参数的局部平面波(图 3a)，最后应用式(6)即可求得对应第61炮记录的FBPSTM成像结果(图 3b).将所有炮记录的偏移结果进行叠加，即可得到图 4a所示的最终FBPSTM成像结果，可以看到其同图 4b所示的KPSTM结果基本一致.对比两者的计算时间，KPSTM总计算时间为276.8 s，FBPSTM倾斜叠加的(时空域算法)计算用时为20.3 s，叠加成像的计算用时为49.8 s，总计算时间为71.3 s，仅为KPSTM的1/4.

2.2 二维实际数据

接下来使用实际地震数据进行测试.该数据共有180炮，炮间隔为40 m，道间距为20 m，其中第5炮地震记录如图 5a所示.均方根速度场(图 5b)共有1201个CDP点，CDP间隔为10 m，时间采样为1800，采样间隔为2 ms.本测试选择的射线束中心间隔为290 m，射线参数采样为27.最终的FBPSTM和KPSTM成像结果分别如图 6a和图 6b所示，可以看到两者成像效果基本一致，仅在成像剖面的极浅层存在细微差异.同样对比两者的计算时间，KPSTM的总计算时间为4560.7 s，FBPSTM中倾斜叠加的(时空域算法)计算用时为157.3 s，叠加成像的计算用时为391.5 s，总计算时间为552.9 s，不足KPSTM的1/8.

3 结论

本文发展了一种快速射线束叠前时间偏移方法，相比于需要逐道运算的Kirchhoff叠前时间偏移，本方法仅需在稀疏的射线束中心位置进行地震数据的局部平面波分解以及后续的成像映射累加运算，可以在保证成像精度的同时，大幅度降低高密度地震资料叠前时间偏移处理的计算成本.此外，本方法还具备良好的灵活性，可以很容易的处理起伏地表和非规则的采集地震数据.