0 引言

随着油气勘探技术的不断发展，勘探目标逐步由简单构造油气藏向隐蔽油气藏、复杂岩性油气藏和非常规油气藏等转变，主要特征表现在：目标地质体的精细化、非均质性变强、横向变速快等方面。叠前深度偏移算法是复杂构造和强横向变速区域成像的关键技术，主要包括：以几何射线理论为基础的偏移算法和以波动方程理论为基础的偏移算法。前者计算效率高，比较灵活，同时能面向目标成像；后者具有较高成像精度，但对速度场的精度要求较高，且计算量较大，尤其是三维数据。

Kirchhoff偏移是业界广泛使用的射线类深度域成像算法，在计算效率和适用性方面有着无可比拟的优势。然而由于渐近射线理论的局限性，不能处理多次波至，存在焦散区和阴影区等缺陷。高斯束偏移不仅能够对高陡构造、复杂断块进行成像，而且继承了射线类算法高效、灵活的特性^[1-4]。逆时偏移首先将地表接收到的地震记录进行逆时延拓，与震源正向传播波场进行互相关以获取成像结果，不受陡倾角的限制，具有很高的成像精度，不仅对速度场的精度要求很高，而且对计算机硬件水平要求较高，在实用性上需要进一步提升。

高斯束逆时偏移是一种通过高斯束的加权积分构建格林函数的逆时偏移算法，既保留了逆时偏移的高精度，又具有高斯束偏移的灵活性、高效性，能够面向目标成像^[5-9]。高斯束偏移算法的成像精度与初始参数的选取有很大关系：当选择较小的初始波束宽度时，浅层构造具有较高的成像精度，但其波束宽度会随着距离的增加而快速变大，影响中深部的成像效果；当选择较大的初始宽度时，波束宽度沿射线路径的变化比较缓慢，不仅降低了中心射线振幅及旅行时的计算精度，而且也会影响最终成像质量。针对上述问题，Nowack^[10]从高斯束的基本性质出发，将高斯束宽度最窄处(束腰)移至地下目标区域，提出了聚焦束的思想，并将其应用于高斯束偏移，提出了聚焦束偏移算法。Nowack^[11]进一步根据多次聚焦思想，将高斯波束沿射线路径约束为相同宽度，实现了动态聚焦束偏移方法，避免了高斯射线束随传播距离的增加而快速发散，提高了复杂构造区域的成像精度。Liu等^[12]通过中心射线求取的菲涅耳带对旅行时场计算范围进行约束。杨继东等^[13]将菲涅耳体的追踪理论应用于高斯束偏移，提出了菲涅耳束理论，实现了复杂地表条件下的菲涅尔束叠前深度偏移。随后，杨继东^[14]基于射线理论和高斯束理论，通过动态的束参数选取策略，将地震波场的能量约束在合理的范围内，发展了一种动态参数控制的波束偏移方法。

本文基于动态聚焦束思想，通过选取合适的束算子，实现格林函数的构建；然后通过格林函数实现正、反向波场的延拓；最后通过正反向波场的互相关求取成像结果，实现了各向异性介质中动态聚焦束逆时偏移。应用数值模型验证了算法的正确性和有效性。

1 基本原理 1.1 动态聚焦束算子

图 1为二维中心射线坐标系，n和t为两个基矢量。对于一条从初始点O(s₀, 0)传播到点R(s_R, 0)的高斯射线束，中心射线上任一点的传播矩阵可以表示为

$ \boldsymbol{\varPi}\left(s_{0}, s_{R}\right)=\left[\begin{array}{ll} Q_{1}\left(s_{R}\right) & Q_{2}\left(s_{R}\right) \\ P_{1}\left(s_{R}\right) & P_{2}\left(s_{R}\right) \end{array}\right] $

(1)

式中：(P₁, Q₁)表示动力学射线追踪方程的平面波解；(P₂, Q₂)表示动力学射线追踪方程的球面波解。

根据互易原理，可以得到点R到初始点O的传播矩阵

$ \boldsymbol{\varPi}\left(s_{R}, s_{0}\right)=\boldsymbol{\varPi}^{-1}\left(s_{0}, s_{R}\right)=\left[\begin{array}{cc} P_{2}\left(s_{R}\right) & -Q_{2}\left(s_{R}\right) \\ -P_{1}\left(s_{R}\right) & Q_{1}\left(s_{R}\right) \end{array}\right] $

(2)

假定该条高斯束在R点聚焦，且束腰处的宽度为L₀，可求得初始点处的动力学射线追踪参数

$ \left[\begin{array}{l} Q\left(s_{0}\right) \\ P\left(s_{0}\right) \end{array}\right]=\boldsymbol{\varPi}\left(s_{R}, s_{0}\right)\left[\begin{array}{c} -\mathrm{i} \omega L_{0}^{2} \\ 1 \end{array}\right] $

(3)

令

$ M\left(s_{0}\right)=\frac{P\left(s_{0}\right)}{Q\left(s_{0}\right)} $

(4)

则复值束参数可以表示为

$ \alpha=\frac{\operatorname{Re}\left[M\left(s_{0}\right)\right]}{\left|M\left(s_{0}\right)\right|^{2}}+\mathrm{i} \frac{\operatorname{Im}\left[M\left(s_{0}\right)\right]}{\left|M\left(s_{0}\right)\right|^{2}}=V_{O} S_{0}-\mathrm{i} V_{O} L_{0}^{2} $

(5)

式中：V_O为初始点速度；S₀为束腰位置到初始点的距离。该束参数的选取仅仅与射线束的聚焦位置有关，不会随射线路径的改变而变化，是一种静态选取方法。一般情况下，通过式(5)确定束参数的高斯束称为聚焦束，R点为聚焦点。

然而聚焦束只能在特定位置处聚焦，不能对整条射线束进行能量约束，因而在应用中存在一定的局限性。为此，Nowack^[11]对复值参数α(s)进行了修改，构建了动态聚焦型传播算子

$ \alpha(s)=\frac{-Q_{2}(s)-\mathrm{i} \omega_{\mathrm{ref}} L^{2}(s) P_{2}(s)}{Q_{1}(s)+\mathrm{i} \omega_{\mathrm{ref}} L^{2}(s) P_{1}(s)} $

(6)

式中：L(s)为动态聚焦束在s处宽度；ω_ref为参考圆频率。通过构建动态聚焦型传播算子，可将地震波束能量约束在有限范围内，进而提高成像质量。

1.2 动态聚焦束表征格林函数

在高斯束逆时偏移方法中，格林函数可以近似表示为

$ G(x)=\int_{0}^{2 {\rm{ \mathsf{ π} }}} \varPhi(\theta) U_{\theta}(s, n) \mathrm{d} \theta $

(7)

式中：Φ(θ)为出射角为θ的高斯束权系数；U_θ(s, n)为高斯束。具体可表示为

$ \left\{\begin{array}{l} U_{\theta}(s, n)=\sqrt{\frac{V(s)}{Q(s)}} \exp \left[-\mathrm{i} \omega \tau(s)+\frac{\mathrm{i} \omega}{2} \frac{P(s)}{Q(s)} n^{2}\right] \\ \varPhi(\theta)=\frac{\mathrm{i}}{4 {\rm{ \mathsf{ π} }}} \sqrt{\frac{\alpha(s)}{V(s)}} \end{array}\right. $

(8)

式中：τ为中心射线走时；V为地震波速度；复值动力学参数P和Q可以表示为

$ \left\{\begin{array}{l} P(s)=\alpha P_{1}(s)+P_{2}(s) \\ Q(s)=\alpha Q_{1}(s)+Q_{2}(s) \end{array}\right. $

(9)

对于动态聚焦束而言，为了达到聚焦的效果，通常按照式(6)的复值束参数进行高斯束数值求解，进而实现格林函数的构建。动态聚焦束本质上是先给出两组实初值求解动力学射线追踪方程组，然后将计算出的解通过复系数进行线性组合以获得最终的动力学射线追踪参数，是一种间接解法。在高斯束偏移中则仅需要对动力学射线追踪方程组进行一次求解，但由于是对复数进行计算，需要分为实部和虚部单独求解，计算量和聚焦束间接法相差不大。

在图 2a所示的速度场CDP=500处分别试射高斯束、聚焦束和动态聚焦束，其形态如图 2b~图 2d所示。可以看出高斯束的宽度随着传播距离的增加而快速发散；聚焦束宽度在束腰处最窄，随着远离束腰位置也呈现出快速发散趋势；而动态聚焦束宽度在整条路径上基本保持不变。

1.3 运动学和动力学射线追踪

为了克服传统基于弹性参数的各向异性射线追踪算法存在的困难，Zhu等^[15-16]基于相速度对各向异性射线追踪方程进行了重新定义，其运动学射线追踪方程为

$ \left\{\begin{array}{l} \frac{\mathrm{d} x_{m}}{\mathrm{~d} \tau}=V_{\mathrm{g}} \\ \frac{\mathrm{d} p_{m}}{\mathrm{~d} \tau}=-\frac{1}{V_{\mathrm{p}}} \frac{\partial V_{\mathrm{p}}}{\partial x_{m}} \end{array}\right. $

(10)

式中：m=1、2，分别代表x方向和z方向；p为慢度；V_p和V_g分别为相速度和群速度，二者在各向异性介质中关系为

$ V_{\mathrm{g}}=\frac{1}{V_{\mathrm{p}}} \frac{\partial V_{\mathrm{p}}}{\partial p_{m}}+V_{\mathrm{p}}^{2} p_{m} $

(11)

在式(10)所示的各向异性运动学射线追踪方程中，其右边函数更加简约，相比传统基于弹性参数的运动学射线追踪方程，更易于编程实现，同时具有更高的计算效率。

Zhu等基于相速度推导的各向异性介质中动力学射线方程为

$ \left\{\begin{array}{l} \frac{\mathrm{d} P_{s}}{\mathrm{~d} \tau}=-A_{s n} Q_{n}-B_{s n} P_{n} \\ \frac{\mathrm{d} Q_{n}}{\mathrm{~d} \tau}=C_{s n} Q_{n}+D_{s n} P_{n} \end{array}\right. $

(12)

式中A_sn、B_sn、C_sn、D_sn为相关系数，具体可表示为

$ \left\{\begin{array}{l} A_{s n}=\frac{1}{V_{\mathrm{p}}} \frac{\partial^{2} V_{\mathrm{p}}}{\partial y_{s} \partial y_{n}} \\ B_{s n}=\frac{\partial^{2} \ln V_{\mathrm{p}}}{\partial y_{s} \partial q_{n}} \\ C_{s n}=\frac{\partial^{2} \ln V_{\mathrm{p}}}{\partial y_{n} \partial q_{s}} \\ D_{s n}=\frac{\partial V_{\mathrm{g} n}}{\partial q_{s}} \end{array}\right. $

(13)

式中：y_s和y_n为中心射线坐标系下坐标；q=∂τ/∂y。

相比基于弹性参数的各向异性动力学射线追踪方程，基于相速度的动力学射线追踪方程更加简洁，仅需计算相速度和群速度的偏导数，避免了复杂的特征值求解，计效率更高。

1.4 正、反向波场构建

二维均匀介质中，格林函数满足

$ \left\{\begin{array}{l} G\left(\boldsymbol{x}, t ; \boldsymbol{x}_{0}, t_{0}\right)=\delta\left(t-t_{0}\right) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right) \\ \left.G\right|_{t<t_{0}} \equiv 0 \end{array}\right. $

(14)

通过Kirchhoff积分可以得到地下偏移域内点x₀在t₀时刻波场

$ \begin{array}{l} U\left(\boldsymbol{x}_{0}, t_{0}\right)=\int_{t_{0}}^{T} \mathrm{~d} t \iint_{\partial \varOmega} \mathrm{d} S \times \\ \ \ \ \ \ \ \ \ {\left[\mathrm{G}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{0}, t_{0}\right) \frac{\partial}{\partial r} U^{(0)}-U^{(0)} \frac{\partial}{\partial r} G\left(\boldsymbol{x}, t ; \boldsymbol{x}_{0}, t_{0}\right)\right]} \end{array} $

(15)

式中：U⁰(x，t)为接收到的地震波场；T为最大接收时间；∂Ω表示闭合空间Ω的边界； $\frac{\partial }{{\partial r}}$表示沿Ω外法线方向求导。

通过聚焦束渐进式G_GB(x, t; x₀, t₀)满足Kirchhoff近似边界条件G_GB(x, t; x₀, t₀)|_z=0=0替代式(15)中的格林函数，反向延拓波场可以表示为

$ \begin{aligned} U_{\mathrm{B}}\left(\boldsymbol{x}_{0}, t_{0}\right) \approx &-2 \int_{t_{0}}^{T} \mathrm{~d} t \iint_{z=0} \mathrm{~d} x \mathrm{d} y \times \\ & U^{(0)}(\boldsymbol{x}, t) \frac{\partial}{\partial z} G_{\mathrm{GB}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{0}, t_{0}\right) \end{aligned} $

(16)

正向传播地震波场的表达式可以通过聚焦束表征的格林函数渐进式获得

$ U_{\mathrm{F}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{s}\right)=\int_{0}^{\infty} \frac{1}{{\rm{ \mathsf{ π} }}} G_{\mathrm{GB}}\left(\boldsymbol{x}, \boldsymbol{x}_{\mathrm{s}} ; \omega\right) f_{\mathrm{F}}(\omega) \mathrm{e}^{-\mathrm{i} \omega t} \mathrm{~d} \omega $

(17)

式中f_F(ω)表示地震子波f(t)的傅里叶变换。

1.5 互相关成像公式

根据反射波成像准则，成像结果可以通过震源处正向传播波场和接收点处反向延拓波场的互相关求取。单炮成像可以表示为

$ I\left(\boldsymbol{x}_{0}, \boldsymbol{x}_{\mathrm{s}}\right)=\int U_{\mathrm{B}}\left(\boldsymbol{x}_{0}, t_{0}\right) U_{\mathrm{F}}\left(\boldsymbol{x}_{0}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \mathrm{d} t_{0} $

(18)

1.6 算法流程

各向异性动态聚焦束逆时偏移主要分为以下五步：

(1) 读入速度场、各向异性参数场和炮记录；

(2) 通过式(10)和式(12)所示的各向异性运动学和动力学射线追踪方程计算中心射线的旅行时、路径和振幅等信息；

(3) 根据式(6)构建动态聚焦束算子，进而实现正、反向波场的延拓(式(16)、式(17))；

(4) 求取正、反向延拓波场的互相关，得到单炮成像结果；

(5) 将所有炮成像结果叠加，得到最终的成像结果。

2 数值实验 2.1 VTI断块模型

首先采用VTI断块模型对算法进行测试。模型速度和各向异性参数场如图 3所示，横向和纵向网格点数分别为1000和550，网格间距分别10m和8m。采用VTI有限差分正演算法合成地震数据，震源是主频为24Hz的Ricker子波。地震记录共250炮，为全孔径接收，炮间距和道间距分别为40m和10m，时间采样点数为4001，间隔为0.8ms。

模型的试算结果如图 4所示。在各向同性动态聚焦束逆时偏移结果(图 4a)中，由于忽略了介质各向异性因素的影响，导致剖面中反射波不能回归于真实位置，而且存在大量的成像噪声，同相轴连续性和聚焦性较差。在各向异性高斯束逆时偏移(图 4b)和动态聚焦束逆时偏移(图 4c)结果中，反射波能够准确归位，绕射波收敛效果较好，对成像噪声也进行了很好的压制，同相轴的连续性和聚焦性也得到了加强，成像质量得到很大的提升。对比可以看出，在本文方法的成像结果中同相轴的连续性和聚焦性更好，能量分布更为均匀(红色箭头所示)，总体成像效果更好。

2.2 TTI复杂构造模型

采用TTI复杂构造模型验证本文算法的适用性。图 5为速度场和各向异性参数场，模型尺寸为15km×4km，网格间距均为10m。地震数据由TTI有限差分正演算法获得，震源是主频为25Hz的Ricker子波。正演地震记录共300炮，每炮1501道接收，炮间距和道间距分别为50m和10m。记录长度为3.0s，采样间隔为1ms。

图 6a为各向同性动态聚焦束逆时偏移成像结果，可以看出由于未考虑各向异性因素，剖面中存在大量的绕射波和成像噪声，同相轴不够连续和聚焦(红色箭头所示)。在VTI动态聚焦束逆时偏移成像结果(图 6b)中，成像质量得到了一定提升，但由于忽略了θ(对称轴倾角)参数的影响，剖面中依然存在成像噪声，未能准确刻画反射界面(红色箭头所示)。图 6c和图 6d分别为TTI高斯束逆时偏移和动态聚焦束逆时偏移的成像结果，可以看出：二者均对地下构造进行了准确的成像，但聚焦束逆时偏移结果层位更加清晰，同相轴更连续、更聚焦，能量分布更均匀。

3 结论与认识

对于高斯束偏移算法，选取合适的初始参数对成像精度有着决定性的影响。当选择较小的初始宽度时，波束宽度会随着距离的增加而快速变大，浅层构造具有较高的成像精度，但中深部的成像效果较差；当选择较大的初始宽度时，波束宽度沿射线路径的变化较缓慢，降低了中心射线振幅及旅行时的计算精度，导致成像质量下降。本文将聚焦束的思想与高斯束逆时偏移方法相结合，通过选取合适的束算子构建格林函数，再通过格林函数构建正、反向延拓波场，最后通过正、反向延拓波场的互相关求取成像结果，实现了各向异性介质中聚焦束逆时偏移。模型数据试算验证了本文方法正确性和适用性。

将本文方法推广到弹性介质是今后的研究方向。