0 引言

随着采集技术的不断进步，越来越多的宽方位和大炮检距数据被采集到，速度各向异性对地震偏移成像的影响变得不可忽略。对于构造成像，目前主要关注的是沉积岩石中的各向异性，其总体可归结为横向各向同性(TI)介质，根据对称轴方向的不同，主要分为VTI介质和TTI介质。

叠前深度偏移是处理强横向变速和复杂构造地区成像的关键技术，依据理论的不同主要分为射线类偏移和波动方程类偏移。Kirchhoff偏移是工业界普遍使用的一种射线类偏移方法，不仅具很高的计算效率，而且能够面向目标成像，然而存在不能处理多次波至、焦散区和阴影区等固有缺陷。高斯束偏移作为一种改进的射线类方法，不仅继承了Kirchhoff偏移灵活、高效的特点，而且克服了其固有缺陷，同时具有接近波动方程的成像精度。高斯束方法起源于20世纪70年代，Červený等^[1]首先将高斯束方法引入到地震波场正演模拟。随着束叠加方法^[2]的提出，高斯束方法才开始应用于偏移成像。Hill^[3]将地震数据进行平面波分解，进而通过高斯束实现平面波的反向延拓成像，实现了高斯束方法在叠后偏移中的应用；Alkhalifah^[4]通过引入各向异性射线追踪方程，将叠后高斯束偏移推广到各向异性介质。随后Hill^[5]采用最速下降法对高斯束偏移方程的多重积分进行优化，并结合Claerbout成像准则^[6]，实现了叠前高斯束深度偏移；为了进一步提升高斯束方法对观测系统的适应性，Nowack等^[7]和Gray^[8]研究了共炮域的叠前高斯束深度偏移算法；Zhu等^[9-10]基于相速度重新定义了各向异性射线追踪方程体系，将叠前高斯束偏移应用于各向异性介质；Gray等^[11]在广义真振幅偏移理论^[12]的基础上，提出了共炮域真振幅叠前高斯束偏移；岳玉波等^[13]、李振春等^[14]对复杂介质高斯束方法进行了研究，并给出了基于高斯束传播角度信息的角度域共成像点道集(ADCIG)直接提取方法；段鹏飞等^[15]通过近似各向异性射线追踪方程组相关系数，并结合角度域成像准则，实现了TI介质局部角度域的叠前高斯束偏移；蔡杰雄等^[16]研究了高斯束偏移和层析速度反演建模；吴娟等^[17]利用高斯束偏移压制鬼波，增强了同相轴连续性，提高了分辨率。

波动方程偏移基于波动方程的数值解法，主要包括单程波偏移和双程波偏移。相比于射线类方法，单程波偏移具有较高的成像精度，但是无法对陡倾地层成像。双程波偏移(即逆时偏移)不受倾角的限制，然而不仅对速度模型的要求比较高，而且计算量巨大。Popov等^[18]提出了采用高斯束构建经典波动方程中格林函数的逆时偏移算法。黄建平等^[19]在Kirchhoff积分的基础上，利用高斯束叠加积分表征格林函数，实现了格林函数高斯束逆时偏移。张凯等^[20]通过引入基于弹性参数的各向异性射线追踪方程，将高斯束逆时偏移推广到VTI介质中，然而算法比较复杂，而且在每步射线追踪时需要对特征值进行求解，计算效率相对较低。

本文将基于相速度的各向异性射线追踪体系引入到高斯束逆时偏移中，并结合高斯束计算过程中的角度信息，实现了一种适用于广义TI介质的角度域高斯束逆时偏移方法。最后通过模型试算对本文方法的正确性和适用性进行了验证。结果表明：在保证成像精度的前提下，本文方法在计算效率方面优势更加明显；提取的ADCIG不仅可以为后续各向异性偏移速度分析提供支撑，而且可以进行分角度叠加成像，以压制成像噪声、提高成像质量。

1 基本原理

各向异性高斯束逆时偏移与各向同性高斯束逆时偏移原理类似：首先通过高斯束叠加积分构建格林函数，进而利用格林函数实现正、反向波场的延拓，最后通过正、反向延拓波场的互相关求取成像值。

1.1 正、反向波场延拓

二维均匀各向同性介质中，空间x处波场U(x，t)满足波动方程

$ LU = \left[ {\Delta - \frac{1}{{{C^2}(\mathit{\boldsymbol{x}})}}\frac{{{\partial ^2}}}{{\partial {t^2}}}} \right]U = 0 $

(1)

式中：Δ为拉普拉斯算子；C(x)为地震波速度。利用Kirchhoff积分可以得到炮检距内x₀点t₀时刻的反向延拓波U(x₀, t₀)。设格林函数G满足

$ \left\{ {\begin{array}{*{20}{l}} {LG\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_0},{t_0}} \right) = {\rm{ \mathsf{ δ} }}\left( {t - {t_0}} \right){\rm{ \mathsf{ δ} }}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_0}} \right)}\\ {{{\left. G \right|}_{t < {t_0}}} = 0} \end{array}} \right. $

(2)

记在0≤t≤T时间间隔内接收的波场为U⁰(x, t)，此时反向延拓波场为

$ \begin{gathered} U\left( {{\mathit{\boldsymbol{x}}_0},{t_0}} \right) = \int_{{t_0}}^T {\text{d}} t\iint_{\partial \Omega } {\left[ {G\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_0},{t_0}} \right)\frac{\partial }{{\partial {n_x}}}{U^0} - } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{U^0}\frac{\partial }{{\partial {n_x}}}G\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_0},{t_0}} \right)} \right]{\text{d}}{S_x} \hfill \\ \end{gathered} $

(3)

式中：∂Ω为闭合区间Ω的边界；$ \frac{\partial }{{\partial {n_x}}}$为沿Ω外法线方向的导数；G(x, t; x₀, t₀)为高斯束表征的格林函数渐近式，在Kirchhoff近似下满足边界条件

$ {\left. {G\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_0},{t_0}} \right)} \right|_{z = 0}} = 0 $

(4)

此时反向延拓波场可以表示为

$ \begin{gathered} {U_{\text{B}}}\left( {{\mathit{\boldsymbol{x}}_0},{t_0}} \right) = - 2\int_{{t_0}}^T {\text{d}} t\iint_{z = 0} {{U^0}}(\mathit{\boldsymbol{x}},t) \times \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{\partial }{{\partial z}}G\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_0},{t_0}} \right){\text{d}}x{\text{d}}y \hfill \\ \end{gathered} $

(5)

利用高斯束表征的格林函数渐近式，可以求得正向传播波场

$ {U_{\text{F}}}\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_{\text{s}}}} \right) = {\rm{Re}} \int_0^\infty {\frac{1}{{\rm{\pi }}}} {G_{{\text{GB}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\text{s}}};\omega } \right){f_{\text{F}}}(\omega ){{\text{e}}^{ - {\text{i}}\omega t}}{\text{d}}\omega $

(6)

式中：x_s为震源坐标；G_GB(x, x_s; ω)为高斯束表征的格林函数渐近式；f_F(ω)为地震子波的傅里叶变换。

1.2 格林函数渐近式

Popov^[21]给出了高斯束表征的格林函数渐近式为

$ G\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_0},{t_0}} \right) = {\rm{Re}} \int_0^\infty {\frac{1}{{\rm{\pi }}}} {G_{{\text{GB}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_0};\omega } \right){{\text{e}}^{ - {\text{i}}\omega \left( {t - {t_0}} \right)}}{\text{d}}\omega $

(7)

式中G_GB(x, x₀; ω)为频率域格林函数。在高斯束逆时偏移中，格林函数是由一系列从x₀出射、具有不同出射角的高斯束的叠加积分表征的，具体表达式为

$ {G_{{\text{GB}}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_0};\omega } \right) = \mathit{\Phi }\int u \left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{p}},\omega } \right){\text{d}}\theta $

(8)

式中：p为中心射线慢度；θ为出射角；u(x, x₀, p, ω)为高斯束；Φ为初始振幅系数。在二维情况下，有

$ \left\{ {\begin{array}{*{20}{l}} {u\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_0},\mathit{\boldsymbol{p}},\omega } \right) = \sqrt {\frac{V}{Q}} \exp \left[ { - {\rm{i}}\omega \tau (s) + \frac{{{\rm{i}}\omega }}{2}\frac{P}{Q}{n^2}} \right]}\\ {\mathit{\Phi } = \frac{{ - {\rm{i}}}}{{4{\rm{ \mathsf{ π} }}V\left( {{\mathit{\boldsymbol{x}}_0}} \right)}}} \end{array}} \right. $

(9)

式中：P和Q为动力学射线追踪参数；V和τ分别为中心射线的速度和旅行时。

1.3 互相关成像

根据反射波成像准则，通过正向传播波场和反向延拓波场的互相关可以得到最终成像

$ I\left( {{\mathit{\boldsymbol{x}}_0},{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \int {{U_{\rm{B}}}} \left( {{\mathit{\boldsymbol{x}}_0},{t_0}} \right){U_{\rm{F}}}\left( {{\mathit{\boldsymbol{x}}_0},t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}{t_0} $

(10)

式中I(x₀, x_s)为单炮成像值，将所有炮成像值叠加即可得到最终的成像结果。

2 各向异性射线追踪

各向异性和各向同性高斯束逆时偏移的本质区别在于射线追踪。运动学射线追踪主要用于计算中心射线的旅行时和路径信息，而动力学射线追踪则是用来求取动力学射线追踪参数。

2.1 各向异性运动学射线追踪

Červený^[22]基于各向异性弹性波波动方程推导了各向异性运动学射线追踪方程为

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}{x_i}}}{{{\rm{d}}\tau }} = {a_{ijkl}}{p_l}{g_j}{g_k}}\\ {\frac{{{\rm{d}}{p_i}}}{{{\rm{d}}\tau }} = - \frac{1}{2}\frac{{\partial {a_{mjkl}}}}{{\partial {x_i}}}{p_m}{p_l}{g_j}{g_k}} \end{array}} \right. $

(11)

式中：p和g分别为慢度矢量和极化矢量；a_ijkl为弹性参数(i, j, k, l, m=1, 2, 3, 分别代表x、y、z三个方向)；τ为中心射线的旅行时。可以看出上式右边函数比较复杂，求解相对困难，计算效率低。为了解决该问题，Zhu等^[9-10]基于相速度对各向异性运动学射线追踪方程进行了重新定义

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}{x_i}}}{{{\rm{d}}\tau }} = V}\\ {\frac{{{\rm{d}}{p_i}}}{{{\rm{d}}\tau }} = - \frac{1}{v}\frac{{\partial v}}{{\partial {x_i}}}} \end{array}} \right. $

(12)

式中：V为群速度；v为相速度。在TI介质中二者满足

$ V = \frac{1}{v}\frac{{\partial v}}{{\partial {p_i}}} + {v^2}{p_i} $

(13)

上述基于相速度的各向异性运动学射线方程组更为简洁，易于实现，具有更高的计算效率。

2.2 各向异性动力学射线追踪

在各向异性介质中，射线中心坐标系不再正交，此时需要引入一个沿射线路径的权值变量消除这种非正交性引起的误差。Hanyga^[23]给出了基于弹性参数的各向异性动力学射线追踪方程组

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}Q}}{{{\rm{d}}\tau }} = MP + NQ}\\ {\frac{{{\rm{d}}P}}{{{\rm{d}}\tau }} = - NP - HQ} \end{array}} \right. $

(14)

式中M、N、H为程函方程特征值对射线中心坐标系中坐标n和慢度p_n的偏导数，定义为

$ \left\{ {\begin{array}{*{20}{l}} {M = \frac{1}{2}\frac{{{\partial ^2}Z}}{{\partial {n^2}}} - \frac{1}{4}{{\left( {\frac{{\partial Z}}{{\partial n}}} \right)}^2}}\\ {N = \frac{1}{2}\frac{{{\partial ^2}Z}}{{\partial {p_n}\partial n}} - \frac{1}{4}\frac{{\partial Z}}{{\partial {p_n}}}\frac{{\partial Z}}{{\partial n}}}\\ {H = \frac{1}{2}\frac{{{\partial ^2}Z}}{{\partial p_n^2}} - \frac{1}{4}{{\left( {\frac{{\partial Z}}{{\partial {p_n}}}} \right)}^2}} \end{array}} \right. $

(15)

式中Z为Christoffel方程的特征值。式(14)所示的方程组不仅在每一步射线追踪时都需要对特征值进行求解，而且需要确定的弱各向异性参数，在适用性上有一定的局限性。Zhu等^[9-10]给出了基于相速度动力学射线追踪方程组

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}{P_J}}}{{{\rm{d}}\tau }} = - {A_{JK}}{Q_K} - {B_{JK}}{P_K}}\\ {\frac{{{\rm{d}}{Q_J}}}{{{\rm{d}}\tau }} = {C_{JK}}{Q_K} + {D_{JK}}{P_K}} \end{array}} \right. $

(16)

式中：J、K为射线中心坐标系下的坐标方向；P_J、P_K和Q_J、Q_K分别为动力学射线追踪参数沿J、K方向的分量；A_JK、B_JK、C_JK、D_JK为相关系数，具体表达式为

$ \left\{ {\begin{array}{*{20}{l}} {{A_{JK}} = \frac{1}{v}\frac{{{\partial ^2}v}}{{\partial {y_J}\partial {y_K}}}}\\ {{B_{JK}} = \frac{{{\partial ^2}\ln v}}{{\partial {y_J}\partial {q_K}}}}\\ {{C_{JK}} = \frac{{{\partial ^2}\ln v}}{{\partial {y_K}\partial {q_J}}}}\\ {{D_{JK}} = \frac{{\partial V}}{{\partial {q_J}}}} \end{array}} \right. $

(17)

其中：y_J和y_K为射线中心坐标系(y_J, y_K, τ)中的坐标；$ {q_J} = \frac{{\partial \tau }}{{\partial {y_J}}}$。

上述基于相速度的动力学射线追踪方程组在求解时仅需要简单计算相速度和群速度的偏导数，不仅避免了每步射线追踪时特征值的求解，而且能够避免求解地下介质弱各向异性时的不确定性，具有更高的计算效率和更广的适用性。

3 ADCIG提取

根据Hill^[3]给出的高斯束初始值，可以确保高斯束的振幅和走时在一个波长范围内平缓的变化，因此可以利用高斯束的实值走时信息直接提取ADCIG。

在二维射线中心坐标系(图 1)下，中心射线附近点R处的实值走时可以表示为

$ T(R) = T(X) + \frac{1}{2} {\rm{Re}} \left[ {\frac{{P(X)}}{{Q(X)}}} \right]{n^2} $

(18)

式中：T(X)分别为各向异性射线追踪求取的旅行时；n是射线附近R点到X点的距离，可以表示为

$ n = \left( {x - {x_X}} \right){t_z} - \left( {z - {z_X}} \right){t_x} $

(19)

其中t_x和t_z为单位切向量t的x分量和z分量。将式(18)两边分别对x求偏导数，可得

$ \frac{{\partial T(R)}}{{\partial x}} = \frac{{\partial T(X)}}{{\partial x}} + n {\rm{Re}} \left[ {\frac{{P(X)}}{{Q(X)}}} \right]\frac{{\partial n}}{{\partial x}} $

(20)

将式(19)代入式(20)，并记

$ \left\{ {\begin{array}{*{20}{l}} {{p_x}(X) = \frac{{\partial T(X)}}{{\partial x}}}\\ {{p_x}(R) = \frac{{\partial T(R)}}{{\partial x}}} \end{array}} \right. $

(21)

可得

$ \begin{array}{l} {p_x}(R) = {p_x}(X) + \left[ {\left( {x - {x_X}} \right)t_z^2 - \left( {z - {z_X}} \right){t_x}{t_z}} \right] \times \\ \;\;\;\;\;\;\;\;\;\;\;{\rm{Re}} \left[ {\frac{{P(X)}}{{Q(X)}}} \right] \end{array} $

(22)

同理可得

$ \begin{array}{l} {p_z}(R) = {p_z}(X) - \left[ {\left( {x - {x_X}} \right){t_x}{t_z} - \left( {z - {z_X}} \right)t_x^2} \right] \times \\ \;\;\;\;\;\;\;\;\;\;\;\;{\rm{Re}} \left[ {\frac{{P(X)}}{{Q(X)}}} \right] \end{array} $

(23)

慢度矢量和传播角度ϕ满足

$ \phi = \left\{ {\begin{array}{*{20}{l}} {\arctan \frac{{{p_x}}}{{{p_z}}} - \pi }&{{p_x} < 0,{p_z} < 0}\\ {\arctan \frac{{{p_x}}}{{{p_z}}} + \pi }&{{p_x} > 0,{p_z} < 0}\\ {\arctan \frac{{{p_x}}}{{{p_z}}} - \pi }&{其他} \end{array}} \right. $

(24)

利用上式计算出炮点和检波点处的传播角度之后，可以得到偏移张角，进而可以实现ADCIG的提取。

4 计算流程

各向异性高斯束逆时偏移具体实现流程可以分为以下五步：

(1) 读入速度场和炮记录，以及相关各向异性参数场；

(2) 通过各向异性运动学(式(12))和动力学(式(16))射线追踪求取中心射线的运动学和动力学信息；

(3) 在震源处利用式(6)计算正向传播波场，同时在接收点处利用式(5)构建反向延拓波场；

(4) 利用式(10)的成像公式，计算震源处正向传播波场和接收点处反向延拓波场的互相关，得到单炮成像值；

(5) 所有炮计算完后，将单炮成像值进行叠加，得到最终的成像结果和ADCIG。

5 模型试算 5.1 VTI介质ZY模型

首先通过VTI介质ZY模型验证角度域高斯束逆时偏移的正确性。模型速度场和各向异性参数场如图 2所示，模型中包含断层、洼陷等。模型纵横向网格点数为550×1000，网格距和CDP间距均为10m。合成地震数据共251炮，时间采样点数为4001，采样间隔为0.8ms。

图 3a为各向同性高斯束逆时偏移结果，可以看出由于忽略了各向异性因素的影响，剖面中存在大量成像噪声，反射界面不能准确归位，同相轴的连续性和聚焦性也比较差。而VTI高斯束逆时偏移(基于弹性参数)(图 3b)和本文方法(图 3c)成像结果的断层和洼陷聚焦性均较好，各向异性层位也得到了清晰的刻画，总体成像质量更高。对比图 3b和图 3c可以见，本文方法结果中同相轴更加聚焦和连续，总体成像剖面更加清晰(图中方框所示)。

为了进一步验证各向异性对成像结果的影响，提取了不同方法CDP500处的ADCIG(图 4)。由图可以看出，各向同性方法ADCIG中同相轴未能准确归位，在大角度处出现了明显的上翘；VTI高斯束逆时偏移(基于弹性参数)和本文方法提取的ADCIG中同相轴不仅能够准确归位，并且整体处于水平状态，而且本文方法提取的ADCIG同相轴更为连续，能量分布更为均匀，信噪比更高。

图 5是本文方法不同角度ADCIG的叠加成像的结果，可以看出，当选用小角度ADCIG叠加成像，陡倾构造的成像结果清晰准确，偏移噪声得到较好压制^[24]。而随着叠加角度的增加，陡倾构造的同相轴能量逐渐减弱，直至消失，但平层构造同相轴构造未被破坏，同相轴能量相对得到增强。

在不考虑回转波的情况下，地层倾角与入射角之和需小于90°，在小角度叠加成像结果中包含比较多的陡倾地层信息，而大角度叠加结果中包含的平层信息较多。因此在叠加成像时，可以根据地下介质的实际构造情况，选择合适的叠加角度进行叠加成像，进而压制成像噪声，提高成像质量。

5.2 TTI洼陷模型

采用洼陷模型测试本文算法在TTI介质中的适用性。图 6为洼陷模型的速度场和各向异性参数场，模型由三个平层和一个洼陷组成，纵横向网格点数为301×1801，网格距和CDP间距均为10m。正演炮记录一共400炮，炮间隔为40m，采用全接收方式，时间记录长度为4.0s，记录间隔为1ms。

图 7a是各向同性高斯束逆时偏移的结果，可以看出成像剖面中普遍存在噪声，底部水平反射同相轴出现上翘，而且同相轴的连续性和聚焦性较差。图 7b是VTI高斯束逆时的偏移结果，虽然成像质量得到了一定的提升，但是由于忽略了倾斜角的影响，反射界面并未正确归位，底层的水平反射同相轴也有轻微的上翘。图 7c和7d分别是基于弹性参数的TTI高斯束逆时偏移结果和本文方法偏移结果，可见反射界面能够正确归位，绕射波也得到了较好的收敛，成像剖面清晰，信噪比高。由于模型比较简单，本文方法和基于弹性参数的TTI高斯束逆时偏移的结果相差不大，只在同相轴能量分布和聚焦性上存在一定差异。

表 1是四种方法的计算效率对比，可以看出各向异性算法需要进行复杂的各向异性参数计算，因此计算时间均高于各向同性算法。但同传统基于弹性参数的各向异性算法相比，在保证成像精度的前提下，本文方法计算效率具有明显的优势。

6 结论

本文将基于相速度的各向异性射线追踪体系应用于高斯束逆时偏移，实现了一种更为高效的TI介质高斯束逆时偏移成像方法，并结合高斯束的传播角度信息，直接提取了ADCIG。模型试算结果表明，与传统基于弹性参数的各向异性算法相比，本文算法在计算效率方面具有明显优势，同时提取的ADCIG不仅能够为后续的各向异性偏移速度分析提供有力支撑，而且可以用于分角度叠加成像，压制成像噪声、提高成像质量。