0 引言

作为一种沿中心射线传播的局部瞬态场，delta波包源自高斯波束，是高斯波束在时空域的对偶表示。2016年，笔者^[1-2]将其引入到偏移领域：利用delta波包叠加表达时间域格林函数，通过互相关成像条件，给出了深度偏移算法。基于delta波包叠加的深度偏移算法既继承了高斯波束偏移的诸多优点，如包含多波至走时、计算焦散区振幅、处理陡倾角等，又保留了射线类偏移方法的缺陷，如高速岩问题。

高速岩是地震勘探中经常遇见的棘手情况。随着陆地和浅海探明油气储量的减少，深海油气勘探正成为世界油气勘探的新热点。全球深水油气资源丰富，据2008年资料统计，40%的油气储量发现于浅水( < 500 m)，29%的储量发现于深水(>500 m)^[3-4]。深海中探明油气储量有相当比重分布在高速岩(如墨西哥湾盐岩、南海火成岩等)下^[4-5]。由于这类岩石的发育和构造样式复杂多变，与围岩间速度差异明显，当射线穿过高速岩体时，受折射效应影响，射线密度降低，进而导致delta波包的分布密度降低，使成像质量变差，甚至无法成像。为弥补这一缺陷，我们需要选择一种可以有效提高射线密度的射线追踪算法。

波前构建算法是一种计算地震波走时、射线路径以及射线振幅的快速算法。与传统射线追踪方法不同，波前构建算法不再局限于单条射线追踪计算，而是随波前推进同时计算多条射线。它最早由Vinje等^[6]提出，采用样条函数进行插值，适用于二维各向同性介质。随后Coman等^[7-8]引入了双三次样条函数差插值方法，韩复兴^[9]引入了二维三次卷积插值方法，这两种方法既满足了计算精度又提高了计算效率。Sun等^[10]和Vinje等^[11]成功实现了三维复杂介质的射线追踪数值计算，孙小东等^[12]利用四面体网格剖分，提高了三维介质射线追踪效率。Gibson等^[13]、白海军等^[14]将波前构建方法拓展到各向异性介质的数值模拟中。

作为一种可同时考虑到整个射线场的计算技巧，波前构建算法在波前推进过程中可随时检验射线场密度，当射线场密度不足时，可通过插入新的射线提高射线覆盖率。本文采用波前构建法的深度偏移算法，在处理具有复杂形状的盐丘模型(如Sigsbee 2A)时，由于射线密度的提升，可有效改善高速体下方区域的成像质量。对于射线插入的标准，在总结前人工作的基础上，采用Coman等^[8]和韩复兴^[9]的描述。在具体实现上，采用链表结构替代数组结构存储波前信息，以节省CPU耗时。

1 偏移方法原理 1.1 波场表达

考虑2维单炮地震记录，震源位置为x_s=(x_s, 0)，地表接收点位置为x_r=(x_r, 0)，地下成像点位置为x=(x, z)，根据Claerbout^[15]提出的互相关成像条件，有

式中：I(x, x_s)为共炮道集下成像点x处的成像结果；t为时间；U_D(x, t; x_s)表示由源点x_s传播到成像点x的直达波场；U_P(x, t; x_r)表示由接收点x_r反向延拓至成像点x的延拓波场。

根据第一种Rayleigh积分和高频渐近理论，延拓波场U_P(x, t; x_r)满足表达式^[16]

式中：v_r为接收点x_r处的地表速度；θ_r为接收点x_r处射线出射方向与z坐标轴正方向的夹角；U⁽⁰⁾(x, t)为x_r处接收到由震源x_s引起的地震反射波场；G(x, -t)表示由接收点x_r传播到成像点x的时间域格林函数，其中，-t表示反因果传播；*表示褶积。

假设直达波场U_D(x, t; x_s)由x_s处点源产生，其震源函数为f(t)，则U_D(x, t; x_s)满足表达式

式中，G(x, t)表示由源点x_s传播到成像点x的时间域格林函数。

1.2 delta波包

格林函数表达式不同，使用的偏移方法也不同。若频率域格林函数取A_ARTexp (iωτ_ART)的形式，可使用常规Kirchhoff偏移，角标ART表示渐近射线理论(asymptotic ray theory)，A_ART为振幅，τ_ART为实走时，ω为频率；若取高斯波束叠加形式，可利用高斯波束叠加法进行深度偏移^[16]；本文利用delta波包叠加表示时间域格林函数，

式中：w_DP(γ, t)被Červený^[17]称为delta波包函数，是单条高斯波束在时空域的对偶形式；角标DP表示delta波包(delta packet)；γ为射线坐标，用以标识中心射线，二维情况下可取极坐标或笛卡尔坐标；区域D为射线坐标γ的积分范围。式(4)表明，将区域D内delta波包叠加可得时间域格林函数G_DP(x, t)。

delta波包是高斯波束在时空域的对偶形式，满足高斯波束的傅里叶逆变换：

式中：权函数Φ、振幅A和走时τ皆为复数，满足高斯波束理论，可利用运动学和动力学射线追踪方程求取^[18]。设复走时τ=τ^R+iτ^I，其中τ^R和τ^I分别为复走时的实部和虚部，并且τ^I≥0。当振幅与频率无关时，波包函数w_DP(γ, t)在时间域存在解析式^[17]：

利用式(6)可直接在时间域计算波包函数w_DP(γ, t)。将公式(6)代入公式(4)，可得delta波包叠加形式的时间域格林函数表达式。

1.3 算法流程

结合之前的推导，利用delta波包叠加进行时间域深度偏移的计算流程为：1)输入光滑速度场和地震反射数据。因为动力学追踪方程含有速度场的二阶导数，所以需要对速度场的二阶导数进行光滑。2)数值求解运动学和动力学射线追踪方程。3)利用公式(3)、(4)和(6)计算直达波场U_D(x, t; x_s)，再利用公(2)、(4)和(6)计算延拓波场U_P(x, t; x_r)。4)将直达波场与延拓波场进行互相关成像，可得到单炮偏移记录。5)重复步骤3)和4)，将所有单炮偏移记录进行叠加。

2 射线追踪的波前构建法

为了将波前构建法与深度偏移算法相结合，需要讨论三个问题：射线插入标准，射线插入方法和波前数据存储。

2.1 插入射线

为了能在波前推进时插入新的射线，需要解决什么时候插入射线和如何插入射线两个问题。对于前一个问题，采用Coman^[8]和韩复兴^[9]的描述，将新射线插入的标准分为三部分：1)同一波前面上相邻射线点间距离大于距离预设值时，插入新射线；2)当不满足条件1)、但相邻射线曲率差大于曲率差预设值时，插入新射线；3)满足条件1)、但两条射线相交时，在交点插入新射线。以上三个条件同时使用时，既避免了小范围内速度异常引起的采样不足，又可控制计算网格节点处的走时误差。

当相邻两条射线满足插入标准时，使用虚拟震源方法进行新射线插入计算^[19]。如图 1所示，射线A1和射线A2为两条相邻射线，M_A1和M_A2分别为它们与波前面的交点，交点处的慢度方向为p_A1和p_A2。当它们满足插入新射线标准时，过点M_A1和M_A2分别沿慢度方向p_A1和p_A2做直线，其直线交点为虚拟震源，直线夹角为α。假设新插入射线由虚拟震源发射，为夹角α的角平分线，根据M_A1和M_A2到虚拟震源的距离和慢度方向，可确定新插入波前点M_new的位置。可以看出，利用虚拟震源方法可以快速方便地插入新射线。但需要注意的是，条件1)中判断射线是否插入的距离预设值不宜过大，特别是当介质复杂、速度变化较大时，以减少该方法的数值误差。

2.2 链表式存储

进行常规Kirchhoff偏移时，我们只需计算网格节点上的走时和振幅；此时波前构建程序可只保留当前计算的两层波前，而将已被计算过的波前释放掉。但利用delta波包叠加进行深度偏移时，需要使用完整射线路径，这要求波前构建程序必须保留所有波前信息，并从中恢复出射线路径。

为了适应波前构建法频繁进行新射线的插入操作，将使用链表结构对波前数据进行存储。链表是一种重要的计算机存储结构，用它表示线性表时，利用指针表示节点间的逻辑关系。因此，一个存储节点通常包含两个部分：数据(data)和指针(next)。数据域用于存储数据，指针域用于存放下一个节点的地址。用链表存储数据的优点在于可采用动态内存分配，不会造成内存浪费和溢出。链表执行插入和删除操作时，只需修改指针即可，无需移动大量元素^[20]。

本文用一条链表存储一条射线，插入新射线等于插入新的链表，同时建立一条射线索引链表用于储存每条射线链表的首节点地址。具体实现方法(图 2)如下：首先，构建一个射线索引链表(图 2中深色阴影部分)。假设有三条射线A1，A2和A3，对应射线索引链表上的3个节点，每个节点主要包含三个地址(以射线A1为例，图 2a)：下条射线(A2)的索引链表节点地址(next_A2)、该条射线(A1)的首个波前点地址(next_A11)，该条射线(A1)的当前计算波前点地址(next_A13)。其次，每条射线对应一条链表(图 2中浅色阴影部分，以射线A1为例)，射线链表的每个节点表示一个波前点，包含两部分(以节点A11为例)：下个节点(A12)的地址(next_A12)和当前节点(A11)的数据(data_A11)。数据data主要包含了当前波前点的走时、慢度、振幅等信息。最后，插入新的射线等于插入新的链表(以射线A2和A3间插入新射线Anew为例，图 2b)，当插入新链表时，只需将射线索引列表中A2节点中下个节点的地址改为next_Anew、Anew节点中下个节点的地址设为next_A3即可。

与常规的数组法相比，链表式存储的优势在于插入射线的效率。假设利用二维数组w[m][n]存储二维波前数据点，m表示波前，n表示射线。当需要在n=i和n=i+1间插入射线时，需要将数组w的一部分w[m][i+1]~w[m][n]整体向右移动到w[m][i+2]~w[m][n+1]，以空出原有的w[m][i+1]用来存储新插入的射线。若利用链表储存，则无需这种移动，只需改变节点地址即可。

此外，利用数组存储时，数组的大小在创建时就已固定。数组太大会造成存储空间的浪费；而数组太小，当内存外溢时，需要将原有数组移动到更大的数组，降低了计算效率。链表式存储的另一个优势就是需要多少内存就可以使用多少内存，充分利用现有内存空间。

3 数值算例

本文选用SMAART公司提供的Sigsbee 2A模型和模拟数据对结合了波前构建法的时间域深度偏移进行试算。图 3为Sigsbee 2A的速度模型。横向网格节点数为2 133，网格间隔为11.43 m；纵向网格节点数为1 201，网格间隔为7.62 m。模拟地震记录共500炮，最大道数为348道，炮间距和道间距均为22.86 m。

可以看出，Sigsbee 2A模型包含一个形状复杂的高速体(盐丘)，其对波场能量的影响非常复杂。受折射效应影响，盐下区域经常出现射线覆盖率不足的情况。这里选取了位于盐丘上的两个出射点(11 659，0)、(16 231，0)进行射线追踪计算(图 4)，射线出射角度范围为[-50°, 50°]。出射点(11 659，0)为第250炮的位置，出射点(16 231，0)为第350炮的位置。为了提高盐下射线覆盖率，射线间隔不采用角度，而采用慢度分量p_x，即任意相邻两条射线间的出射角角度差是不固定的，但慢度分量p_x差是固定的。同时选取两种方法进行射线追踪计算：一种是缩小射线间隔，直接在震源对射线条数进行加密(初始射线200条，如图 4a、d)；另一种是之前讨论的波前构建法(初始射线100条，插入100条，共200条射线，如图 4b、e；初始射线50条，插入150条，共200条射线，如图 4c、f)。通过对比可以看出，同样是200条射线，与加密射线方法相比，利用波前构建方法进行射线追踪可有效提高盐下区域的射线密度，提升射线覆盖率。但同样是波前构建方法，初始射线不同，射线加密区域却略有区别。这是因为，在我们的波前构建算法中，射线还无法实现精准插入，只能根据波前推进情况进行自动插入，所以初始射线的位置和条数将影响射线追踪的最终结果。虽然射线路径略有不同，但出射点下方对应盐下区域的射线覆盖率都得到了显著提升。

表 1给出了数组法和链表法的CPU耗时对比。可以看出：当不插入射线时(初始200条)，数组法要略优于链表法，这主要因为数组法寻址要快于链表法；而当需要插入射线时，链表法要明显优于数组法，CPU耗时至少缩短9%，插入新射线越多，CPU耗时越短。

图 5给出了两种不同方法试算出的偏移成像结果。对比图 5a、b、c可知：由于盐下射线覆盖率的提升，结合了波前构建法的偏移结果，其盐下区域的成像能量得到了显著提高；而其他射线覆盖率本就充足的地方，成像结果变化不大。这意味着当射线密度低时，提升射线密度可有效改善成像质量。但需要注意的是，这种成像能量的提升是无法区分有效信号和噪声的，即有效信号和噪声同时得到了加强。

对比图 5b和图 5c可知，虽然插入射线的条数增多，但成像能量并没有明显提高。这意味着当总射线条数固定时，一味提高插入射线数量，并不会无休止地改善成像质量；相反，由于使用虚拟震源法插入新射线，射线插入越多，其引起的数值误差也越大。

4 结论

1)本文将波前构建法与基于delta波包叠加的时间域深度偏移相结合，既保留了偏移算法的原有优点，如包含多波至走时、计算焦散区振幅、处理陡倾角等，同时又有效提升了射线覆盖率不足区域的射线密度，改善了成像质量。

2)由于深度偏移算法需要完整射线路径，采用链表结构进行波前数据存储，可简洁地进行射线插入操作，并从中轻易恢复出射线路径。

3) Sigsbee 2A模型的数值试算表明，改进后的偏移算法对于复杂盐丘模型具有良好的成像能力，并且利用链表存储波前信息要比利用数组至少节省9%的CPU耗时。但我们同时注意到，由于使用虚拟震源法插入射线，当插入射线过多时，其引起的数值误差可能降低图像的信噪比。如何更精确、有选择地插入射线，以提高波前构建法的计算精度和计算效率，还需我们做更进一步的研究。

致谢: 吉林大学移动平台探测研发中心对本文提供了计算设备支持，在此表示感谢。