0 引言

随着采集方位角分布范围的变宽和观测排列长度的增加，各向异性对地震波传播的影响越来越不容忽略。因此，需根据观测系统特性，在速度建模与偏移成像过程中充分考虑各向异性。对各向异性构造成像而言，现今主要研究目标是沉积岩石形成的各向异性，其总体可归结为TI介质(横向各向异性)模型，依据对称轴不同主要分为VTI(垂直对称轴)和TTI介质(倾斜对称轴)。

作为处理复杂构造和强横向变速地区成像的关键技术，叠前深度偏移在近年取得快速发展。高斯束偏移是一种改进的射线类成像方法，不仅具有很高的成像精度，且还能处理多次波至和阴影区等问题，因此更具研究价值和应用前景。

高斯束方法起源于20世纪70年代，最早应用于波场正演^[1-2]。Hill^[3]通过倾斜叠加对地震数据进行局部平面波分解，再利用高斯束对平面波反向延拓、成像，实现了叠后高斯束偏移。Hale^[4]提出最近点搜索法及粗网格递归算法，提高了高斯束偏移方法计算效率。Alkhalifah^[5]通过引入基于弹性参数的各向异性射线追踪算法，将叠后高斯束偏移应用于各向异性介质。Hill^[6]采用最速下降法简化了高斯束偏移方程中的多重积分，提出适用于共方位角、共炮检距道集的叠前高斯束偏移。Nowack等^[7]和Gray^[8]研究了适用于构造成像的共炮域叠前高斯束偏移，提升了高斯束方法对观测系统的适应性。Zhu等^[9-10]基于相速度重新定义了各向异性射线追踪方程组，将叠前高斯束偏移方法推广到各向异性介质。Gray等^[11]将保幅共炮域偏移理论^[12]与高斯束方法相结合，完成了共炮域的保幅高斯束偏移。李振春等^[13]实现了基于互相关成像条件的保幅高斯束偏移，并结合高斯束传播角度信息直接提取角度域共成像点道集。岳玉波等^[14]研讨了复杂介质条件下的高斯束方法。段鹏飞等^[15]讨论了TI介质局部角度域高斯束成像方法，给出一种更高效的高斯束合成方案。高成等^[16]通过调整积分顺序，将频率域的积分变换到时间域，实现了时间域高斯束偏移，并通过模型试算验证了该方法计算效率和成像精度。

随后杨晶等^[17]在局部倾斜叠加时引入汉宁窗进行滤波，提升了时间域高斯束方法对低信噪比数据的适用性。吴娟等^[18]利用高斯束压制鬼波，达到进一步提高同相轴的连续性和分辨率的目的。石星辰等^[19]通过引入品质因子Q进行衰减补偿，实现了矢量黏弹性衰减补偿高斯束偏移。吕考考等^[20]研究了数据驱动的控制束偏移方法，有效地压制了偏移假象和偏移噪声；张瑞等^[21]基于干扰信号与有效信号在τ-p域的差异，发展了一种相似系数阈值滤波的数据驱动控制束偏移算法。Yue等^[22]研究了最小二乘弹性波高斯束偏移方法。

本文在前人研究的基础上，将基于相速度的各向异性射线追踪算法与时间域高斯束偏移相结合，同时在地震波传播方向上考虑各向异性的影响，优化了动力学射线追踪方程相关系数，提出了一种更高效的各向异性声波时间域高斯束偏移方法。通过模型试算对该方法的正确性和适用性进行了验证。与其他各向异性算法相比，在保证成像精度的前提下，本文方法在计算效率和适用性方面具有很大优势。

1 基本原理 1.1 时间域高斯束成像

对于时间域高斯束偏移成像，各向异性与各向同性的基本原理类似，本质区别在于射线追踪。时间域高斯束偏移首先通过高斯束的叠加积分构建格林函数，进而利用格林函数实现正反向波场的延拓，最后求取正反向延拓波场的互相关以获取成像值。

频率域互相关成像条件^[23]可表示为

$ I(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}) = - {\rm{i}}\int {{U_{\rm{r}}}} (\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega )U_{\rm{s}}^*(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ){\rm{sgn}}(\omega ){\rm{d}}\omega $

(1)

式中：“*”为复共轭; sgn(·)为符号函数；U_r(x, x_s; ω)和U_s^*(x, x_s; ω)分别为检波点和震源点波场，基于格林函数的渐近式如下^[11]

$ \begin{array}{l} \begin{array}{*{20}{l}} {{U_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ) = 2\int {\sqrt { - {\rm{i}}\omega } } \frac{{{\rm{cos}}{\theta _{\rm{r}}}}}{{{V_{\rm{r}}}}}{U_{\rm{r}}}({\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ){A_{\rm{r}}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}( - {\rm{i}}\omega {\tau _{\rm{r}}}){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}} \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = - 2\int {\rm{d}} {\mathit{\boldsymbol{x}}_{\rm{r}}}{U_{\rm{r}}}({\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega )\frac{{{\rm{cos}}{\theta _{\rm{r}}}}}{{{V_{\rm{r}}}}}{\rm{i}}\omega {G^*}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ) \end{array} $

(2)

$ \begin{array}{*{20}{l}} {U_{\rm{s}}^*(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ) = - 2\frac{{{\rm{cos}}{\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}\sqrt { - {\rm{i}}\omega } {A_{\rm{s}}}{\rm{exp}}( - {\rm{i}}\omega {\tau _{\rm{s}}})}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = 2\frac{{{\rm{cos}}{\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}{\rm{i}}\omega {G^*}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega )} \end{array} $

(3)

式中：V_s、V_r分别为震源点、检波点处的速度；θ_s、θ_r分别为震源点、检波点处的出射角；A_s、A_r和τ_s、τ_r分别为通过渐近射线理论求得的震源点、检波点处的实值振幅和旅行时；G^*(x, x_s; ω)和G^*(x, x_r; ω)分别为震源点和检波点在成像点x处格林函数的复共轭，具有如下形式^[11]

$ {{G^*}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ) = \frac{{\rm{i}}}{{4\pi }}\int {A_{\rm{s}}^*} {\rm{exp}}( - {\rm{i}}\omega T_{\rm{s}}^*)\frac{{{\rm{d}}{p_x}}}{{{p_z}}}} $

(4)

$ {{G^*}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ) = \frac{{\rm{i}}}{{4\pi }}\int {A_{\rm{r}}^*} {\rm{exp}}( - {\rm{i}}\omega T_{\rm{r}}^*)\frac{{{\rm{d}}{p_x}}}{{{p_z}}}} $

(5)

其中T_s和T_r分别为震源点和检波点到成像点x处的旅行时；p_x、p_z为成像点处的水平、垂直慢度。

由于式(5)需对每个x_r点处的格林函数求解，计算量巨大。根据高斯束波前在初始位置曲率为零的特性，通过引入高斯窗对地震波场进行局部倾斜叠加，再采用高斯束进行延拓，可显著地提高计算效率。

高斯窗函数具有如下性质^[14]

$ \frac{{\sqrt 3 }}{{4\pi }}\left| {\frac{\omega }{{{\omega _{\rm{r}}}}}} \right|{\left( {\frac{{\Delta L}}{{{w_0}}}} \right)^2}\sum\limits_\mathit{\boldsymbol{L}} {{\rm{exp}}} \left( { - \left| {\frac{\omega }{{{\omega _{\rm{r}}}}}} \right|\frac{{|{\mathit{\boldsymbol{x}}_{\rm{r}}} - \mathit{\boldsymbol{L}}|}}{{2w_0^2}}} \right) \approx 1 $

(6)

式中：L为束中心位置(以下下标“L”对应束中心)；ΔL为束中心间隔；ω_r为参考频率；w₀为初始宽度。

局部倾斜叠加公式为^[6]

$ \begin{array}{l} {D_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}};\omega ) = {\left| {\frac{\omega }{{{\omega _{\rm{r}}}}}} \right|^{1/2}}\int {\rm{d}} {\mathit{\boldsymbol{x}}_{\rm{r}}}U({\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ) \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}\left[ {{\rm{i}}\omega {p_{{\rm{L}}x}}({\mathit{\boldsymbol{x}}_{\rm{r}}} - \mathit{\boldsymbol{L}}) - \left| {\frac{\omega }{{{\omega _{\rm{r}}}}}} \right|\frac{{{{({\mathit{\boldsymbol{x}}_{\rm{r}}} - \mathit{\boldsymbol{L}})}^2}}}{{2\omega _0^2}}} \right] \end{array} $

(7)

将震源点和检波点波场代入式(1)，可得

$ \begin{array}{*{20}{c}} {I(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}) = - 4\int {\frac{{{\rm{cos}}{\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}} {G^*}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega ){\rm{i}}\omega {\rm{sgn}}(\omega ){\rm{d}}\omega \times }\\ {\int {{U_{\rm{r}}}} ({\mathit{\boldsymbol{x}}_{\rm{r}}},{\mathit{\boldsymbol{x}}_{\rm{s}}};\omega )\frac{{{\rm{cos}}{\theta _{\rm{r}}}}}{{{V_{\rm{r}}}}}{\rm{i}}\omega {G^*}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}} \end{array} $

(8)

将式(4)~式(7)代入式(8)，可得

$ \begin{array}{*{20}{l}} {I(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}) = \frac{{\Delta L{\omega _{\rm{r}}}}}{{4\sqrt 2 {\pi ^{5/2}}{\omega _0}}}\sum\limits_\mathit{\boldsymbol{L}} {\int {\frac{{{\rm{cos}}{\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}} } {\rm{i}}\omega {\rm{d}}\omega \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int {A_{\rm{s}}^*} {\rm{exp}}( - {\rm{i}}\omega T_{\rm{s}}^*)\frac{{{\rm{d}}{p_{{\rm{s}}x}}}}{{{p_{{\rm{s}}z}}}}\int {\frac{{{\rm{cos}}{\theta _{\rm{L}}}}}{{{V_{\rm{L}}}}}} {\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int {A_{\rm{L}}^*} {\rm{exp}}( - {\rm{i}}\omega T_{\rm{L}}^*)\frac{{{\rm{d}}{p_{{\rm{L}}x}}}}{{{p_{{\rm{L}}z}}}}{D_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},\omega )} \end{array} $

(9)

式中：p_sx和p_sz分别为震源点处的水平慢度和垂直慢度；p_Lx和p_Lz分别为束中心处水平慢度和垂直慢度；A_L、V_L、T_L分别为束中心处振幅、速度和旅行时。

通过改变积分顺序，并令：A=A_sA_L，T=T_s+T_L；A_R=Re(A)，A_I=Im(A)；T_R=Re(T)，T_I=Im(T)。成像公式^[16]则可写为

$ \begin{array}{*{20}{l}} {I(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}) = \frac{{\Delta L{\omega _{\rm{r}}}}}{{4\sqrt 2 {\pi ^{5/2}}{\omega _0}}}\sum\limits_\mathit{\boldsymbol{L}} {\int {\frac{{{\rm{d}}{p_{{\rm{s}}x}}}}{{{p_{{\rm{s}}z}}}}} } \frac{{{\rm{d}}{p_{{\rm{L}}x}}}}{{{p_{{\rm{L}}z}}}}\frac{{{\rm{cos}}{\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}\frac{{{\rm{cos}}{\theta _{\rm{L}}}}}{{{V_{\rm{L}}}}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {I_1}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{r}}})} \end{array} $

(10)

式中

$ \begin{array}{l} \begin{array}{*{20}{l}} {{I_1}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_\text{s}}) = \int {\rm{i}} \omega {\rm{d}}\omega {A^*}{\rm{exp}}( - {\rm{i}}\omega {T_{\rm{R}}}) \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}( - \omega {T_{\rm{I}}}){D_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},\omega )}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \int {\rm{i}} \omega {\rm{d}}\omega {D_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},\omega ){\rm{exp}}( - \omega {T_{\rm{I}}}) \times } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{\rm{exp}}( - {\rm{i}}\omega {T_{\rm{R}}}){A_{\rm{R}}} + }\\ {\int {\rm{i}} \omega {\rm{d}}\omega {D_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},\omega ){\rm{exp}}( - \omega {T_{\rm{I}}}) \times }\\ {{\rm{exp}}( - {\rm{i}}\omega {T_{\rm{R}}})( - {\rm{i}}{A_{\rm{I}}})} \end{array} \end{array} $

(11)

傅里叶逆变换定义为

$ f(t) = \frac{1}{{2\pi }}\int {F(\omega )} {\rm{exp}}( - {\rm{i}}\omega t){\rm{d}}\omega $

(12)

此时有

$ \begin{array}{l} 2\pi {A_{\rm{R}}}{{\tilde D}_{\rm{s}}}(L,{p_{{\rm{L}}x}},{T_{\rm{R}}},{T_{\rm{I}}}) = \int {\rm{i}} \omega {\rm{d}}\omega {D_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},\omega ) \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}( - \omega {T_{\rm{I}}}){\rm{exp}}( - {\rm{i}}\omega {T_{\rm{R}}}){A_{\rm{R}}} \end{array} $

(13)

式中$ {\tilde D_{\rm{s}}}(\mathit{\boldsymbol{L}}, {p_{{\rm{L}}x}}, {T_{\rm{R}}}, {T_{\rm{I}}}) $为局部倾斜叠加的傅里叶逆变换。

Aki等^[24]给出的Hilbert变换公式如下

$ \begin{array}{*{20}{l}} {{f_{\rm{H}}}(t) = H[f(\omega )] = \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {f(\omega )} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}[ - {\rm{i}}\frac{\pi }{2}{\rm{sgn}}(\omega )]{\rm{exp}}( - {\rm{i}}\omega t){\rm{d}}\omega } \end{array} $

(14)

此时有

$ \begin{array}{*{20}{l}} {2\pi {A_{\rm{I}}}{{\tilde D}_{{\rm{S,H}}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},{T_{\rm{R}}},{T_{\rm{I}}})}\\ {\;\;\;\,\begin{array}{*{20}{l}} { = \int {\rm{i}} \omega {\rm{d}}\omega {D_{\rm{S}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},\omega ){\rm{exp}}( - \omega {T_{\rm{I}}}) \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{exp}}( - {\rm{i}}\omega {T_{\rm{R}}})( - {\rm{i}}{A_{\rm{I}}})} \end{array}} \end{array} $

(15)

式中$ {\tilde D_{{\rm{S, H}}}}(\mathit{\boldsymbol{L}}, {p_{{\rm{L}}x}}, {T_{\rm{R}}}, {T_{\rm{I}}}) $为局部加窗倾斜叠加的Hilbert变换。

将式(13)、式(15)代入式(10)，可以得到最终的成像公式

$ \begin{array}{l} I(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}}) = \frac{{\Delta L{\omega _{\rm{r}}}}}{{2\sqrt 2 {\pi ^{3/2}}{\omega _0}}}\sum\limits_L {\int {\frac{{{\rm{d}}{p_{{\rm{s}}x}}}}{{{p_{{\rm{s}}z}}}}} } \frac{{{\rm{d}}{p_{{\rm{L}}x}}}}{{{p_{{\rm{L}}z}}}}\frac{{{\rm{cos}}{\theta _{\rm{s}}}}}{{{V_{\rm{s}}}}}\frac{{{\rm{cos}}{\theta _{\rm{L}}}}}{{{V_{\rm{L}}}}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [{A_{\rm{R}}}{{\tilde D}_{\rm{s}}}(\mathit{\boldsymbol{L}},{p_{{\rm{L}}x}},{T_{\rm{R}}},{T_{\rm{I}}}) + {A_{\rm{I}}}{{\tilde D}_{{\rm{S,H}}}}(L,{p_{{\rm{L}}x}},{T_{\rm{R}}},{T_{\rm{I}}})] \end{array} $

(16)

从式(16)可知，时间域高斯束偏移通过改变积分顺序和引入相应变换，将频率域计算转换到时间域进行，降低了积分维度，提升了计算效率。

1.2 基于相速度的各向异性射线追踪算法

高斯束方法可分为高斯束的求解和成像结果的叠加两大步骤。高斯束的求解主要通过各向异性射线追踪来实现：运动学射线追踪通常用于求取中心射线的路径、旅行时；动力学射线追踪则用来计算振幅、相位等动力学相关信息。

Zhu^[9]给出了基于相速度的各向异性介质的运动学射线追踪方程

$ \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. $

(17)

式中：V为相速度；v为群速度；i为射线序号。

在TI介质中：$ V = \frac{1}{v}\frac{{\partial v}}{{\partial {p_i}}} + {v^2}{p_i} $。上述方程组可通过四阶Runge-Kutta进行求解。与基于弹性参数的各向异性运动学射线追踪方程组^[25]相比，式(17)不仅具有较高的计算效率，而且避免了在每一步射线追踪时特征值的计算。

在各向异性介质中，射线中心坐标系不再正交。Hanya^[26]通过引入一个沿射线的权重消除了非正交性带来的误差，并给出基于弹性参数的各向异性动力学射线追踪方程组。然而，其求解较困难，即计算效率相对低下。Zhu等^[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. $

(18)

式中：P_j和Q_j为各向异性动力学射线追踪参数；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}{\rm{ln}}v}}{{\partial {y_j}\partial {y_k}}}}\\ {{C_{jk}} = \frac{{{\partial ^2}{\rm{ln}}v}}{{\partial {y_j}\partial {q_k}}}}\\ {{D_{jk}} = \frac{{\partial {V_k}}}{{\partial {q_j}}}} \end{array}} \right. $

(19)

式中：y_j和y_k是射线中心坐标系(y_j, y_k, τ)中的坐标；$ {q_j} = \frac{{\partial \tau }}{{\partial {y_j}}} $。

上述基于相速度的动力学方程组，只需简单计算群速度和相速度的偏导数，易于实现，且还消除了求解地下介质弱各向异性参数时的不确定性，具有很强适用性。然而，其系数依然较复杂。因此，本文在传播方向上考虑各向异性的影响，沿相速度方向求高斯束振幅，相关系数简化为

$ \left\{ {\begin{array}{*{20}{l}} {{A_{jk}} = \frac{1}{v}\frac{{{\partial ^2}v}}{{\partial {y_j}\partial {y_k}}}}\\ {{B_{jk}} = 0}\\ {{C_{jk}} = 0}\\ {{D_{jk}} = \frac{{\partial {V_k}}}{{\partial {q_j}}}} \end{array}} \right. $

(20)

2 模型试算 2.1 Hess模型

首先通过Hess模型测试本文方法在VTI介质中的适用性。图 1展示Hess模型的速度场和各向异性参数场(ε、δ)，可见模型中包含高速岩丘、断层、尖灭等复杂构造。模型纵向和横向采样点分别为1500和3617，网格间距均为20ft。地震记录共720炮，每炮656道接收，记录长度为7.992s，间隔为6ms。

图 2是Hess模型试算结果。可见因忽略了各向异性因素的影响，在各向同性时间域高斯束偏移结果(图 2a)中存在大量的成像噪声，反射波也难以准确归位，同相轴聚焦性也不佳，总体成像质量较差。而通过基于相速度的VTI频率域高斯束偏移(图 2b)、基于弹性参数的VTI时间域高斯束偏移(图 2c)和本文方法(图 2d)所得结果中反射界面能正确归位，绕射波能较好收敛，断层、尖灭等复杂构造也得到精细成像，整体剖面特征显得更清晰，信噪比更高。

为便于对比，将图 2中红框所示区域放大显示。从各向同性时间域高斯束偏移局部放大结果(图 3a)可见，剖面中绕射波不能很好收敛，同相轴的连续性和聚焦性较差，红色箭头所示的各向异性层也未能准确成像；而在基于相速度的VTI频率域高斯束偏移(图 3b)、基于弹性参数的VTI时间域高斯束偏移(图 3c)和本文方法(图 3d)结果中，反射波能正确归位，高速岩丘、断层等复杂构造得到了清晰刻画，箭头所指的各向异性层位也得到准确成像，同相轴连续性和聚焦性较好，成像质量更高。表 1是四种偏移方法的计算效率对比。从单炮计算时间对比可知，各向异性类算法因要进行各向异性参数的计算，故计算时间均高于各向同性算法。但相比于其他各向异性算法，本文方法在确保成像精度的前提下，在计算效率方面具有明显优势。

2.2 TTI逆冲断层

选取TTI逆冲断层对本文方法进行验证。模型速度场和各向异性参数场如图 4所示，可见模型由两个水平反射层和一个逆冲断层组成，纵、横向网格点分别为300和1201，网格间距均为10m。采用TTI有限差分正演算法合成地震记录，地震子波采用主频为25Hz雷克子波。正演地震记录共有300炮，炮间隔为40m；每炮1201道接收，道间隔为10m；时间采样点为4001，间隔1ms。

首先采用各向同性时间域高斯束进行偏移(图 5a)，因未考虑各向异性参数的影响，导致剖面中存在大量成像噪声和绕射波，底部水平反射界面有轻微上翘，逆冲断层两侧的同相轴较难聚焦。审视基于相速度的VTI时间域高斯束偏移结果(图 5b)，虽然可见绕射波有一定的收敛，逆冲断层两侧的同相轴也更聚焦和连续，但由于忽略了倾斜角的影响，导致反射界面未能准确成像，剖面中还存在一些成像噪声，底部水平同相轴的连续性较差。而在基于相速度的TTI频率域(图 5c)和时间域(5d)高斯束偏移结果中，绕射波得到充分收敛，反射界面都能准确归位，逆冲断层得到精细刻画，底部水平同相轴的连续性和聚焦性更好，信噪比更高，成像剖面整体更清晰。

3 结束语

本文将基于相速度的各向异性射线追踪算法应用于时间域高斯束偏移，并在传播方向上考虑各向异性参数影响，优化了动力学射线追踪相关系数，提出一种更高效的各向异性声波时间域高斯束偏移方法。与传统基于弹性参数各向异性射线追踪算法相比，基于相速度的射线追踪算法不仅能避免每步射线追踪时特征值的求解问题，而且还消除了求取地下介质弱各向异性参数时的不确定性，具有更高计算效率和更强适用性。模型试算结果表明：相比于其他各向异性算法，在确保成像精度前提下，本文方法在计算效率方面具明显优势。后续拟将该方法推广到三维，以进一步拓展其适用性。