石油地球物理勘探  2020, Vol. 55 Issue (4): 813-820  DOI: 10.13810/j.cnki.issn.1000-7210.2020.04.013
0
文章快速检索     高级检索

引用本文 

秦宁. 声波各向异性时间域高斯束叠前深度偏移. 石油地球物理勘探, 2020, 55(4): 813-820. DOI: 10.13810/j.cnki.issn.1000-7210.2020.04.013.
QIN Ning. Time-domain Gaussian beam prestack depth migration for acoustic anisotropic media. Oil Geophysical Prospecting, 2020, 55(4): 813-820. DOI: 10.13810/j.cnki.issn.1000-7210.2020.04.013.

本项研究受国家科技重大专项“渤海湾盆地精细勘探关键技术”(2016ZX05006)、中国石化优秀青年科技创新基金项目“潜山内幕散射波叠前深度偏移技术研究”(P18028-1)及胜利油田科技攻关项目“地震层间多次波预测及压制技术研究”(YKW1901)联合资助

作者简介

秦宁  高级工程师, 博士, 1985年生; 2008年获中国石油大学(华东)勘查技术与工程专业学士学位, 2013年获该校地质资源与地质工程专业博士学位, 2015年从胜利石油管理局博士后科研工作站出站; 现就职于中国石化胜利油田物探研究院, 主要从事地震数据处理、速度建模及深度偏移等方面研究

秦宁, 山东省东营市东营区北一路210号胜利物探研究院处理二室, 257022。Email:geoqin@163.com

文章历史

本文于2019年12月19日收到,最终修改稿于2020年5月24日收到
声波各向异性时间域高斯束叠前深度偏移
秦宁     
中国石化胜利油田物探研究院, 山东东营 257022
摘要:作为一种改进的射线类偏移方法,高斯束偏移不仅具备射线类偏移灵活、高效的特性,而且还拥有接近波动类偏移的成像精度。与传统频率域高斯束偏移相比,时间域高斯束偏移通过调整成像公式,把频率域积分转换为成像时刻的函数,提升了计算效率。随着地震采集方位角的变宽和观测排列的加长,各向异性对地震偏移成像的影响变得不容忽略。文中将基于相速度的各向异性射线追踪算法应用于时间域高斯束偏移,并对动力学射线追踪相关系数进行优化,形成一种更高效的各向异性声波时间域高斯束偏移方法。模型试算表明,相比于其他各向异性算法,在保证成像精度前提下,所提方法具有更高计算效率和更强适用性,尤其适用于各向异性复杂构造成像。
关键词时间域    高斯束    深度偏移    TI介质    射线追踪    
Time-domain Gaussian beam prestack depth migration for acoustic anisotropic media
QIN Ning     
Geophysical Research Institute, Sinopec Shengli Oilfield Company, Dongying, Shandong 257022, China
Abstract: As an improved ray migration method, Gaussian beam migration not only inherits the efficiency and flexibility of ray methods, but also has the imaging accuracy close to wave equation migration methods. Compared with traditional Gaussian beam migration in frequency domain, the method in time domain greatly improves computational efficiency by adjusting the imaging formula and converting frequency integral into the function of imaging time. As seismic azimuth becomes larger and recording spread becomes longer, the influence of anisotropy on seismic migration can't be ignored any longer. This paper introduces the anisotropic ray tracing algorithm based on phase velocity to time-domain Gaussian beam migration and optimizes the dynamic ray tracing correlation coefficient, and finally proposes a more efficient method for acoustic anisotropic media. Model tests have demonstrated that the method proposed has higher computational efficiency and better applicability while maintaining higher imaging accuracy than similar anisotropic algorithms, especially it can be used for imaging complex anisotropic structures.
Keywords: time domain    Gaussian beam    depth migration    TI media    ray tracing    
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(·)为符号函数;Ur(x, xs; ω)和Us*(x, xs; ω)分别为检波点和震源点波场,基于格林函数的渐近式如下[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)

式中:VsVr分别为震源点、检波点处的速度;θsθr分别为震源点、检波点处的出射角;AsArτsτr分别为通过渐近射线理论求得的震源点、检波点处的实值振幅和旅行时;G*(x, xs; ω)和G*(x, xr; ω)分别为震源点和检波点在成像点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)

其中TsTr分别为震源点和检波点到成像点x处的旅行时;pxpz为成像点处的水平、垂直慢度。

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

高斯窗函数具有如下性质[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为参考频率;w0为初始宽度。

局部倾斜叠加公式为[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)

式中:psxpsz分别为震源点处的水平慢度和垂直慢度;pLxpLz分别为束中心处水平慢度和垂直慢度;ALVLTL分别为束中心处振幅、速度和旅行时。

通过改变积分顺序,并令:A=AsALT=Ts+TLAR=Re(A),AI=Im(A);TR=Re(T),TI=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)

式中:PjQj为各向异性动力学射线追踪参数;AjkBjkCjkDjk为互相关系数,分别表示为

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

式中:yjyk是射线中心坐标系(yj, yk, τ)中的坐标;$ {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。

图 1 Hess模型 (a)速度场;(b)各向异性参数ε;(c)各向异性参数δ

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

图 2 Hess模型试算结果 (a)各向同性时间域高斯束偏移;(b)VTI频率域高斯束偏移;(c)VTI时间域高斯束偏移(基于弹性参数);(d)本文方法

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

图 3 局部结果对比 (a)各向同性高斯束偏移;(b)VTI频率域高斯束偏移;(c)VTI时间域高斯束偏移(基于弹性参数);(d)本文方法

表 1 计算效率对比
2.2 TTI逆冲断层

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

图 4 逆冲断层模型的速度场(a)和各向异性参数ε(b)、δ (c)及与对称轴倾角θ(d)

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

图 5 逆冲断层模型试算结果 (a)各向同性时间域高斯束偏移;(b)VTI时间域高斯束偏移;(c)TTI频率域高斯束偏移;(d)TTI时间域高斯束偏移
3 结束语

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

参考文献
[1]
Popov M M. A new method of computation of wave fields using Gaussian beams[J]. Wave Motion, 1982, 4(1): 85-97. DOI:10.1016/0165-2125(82)90016-6
[2]
Červený V, Vaněk R J. Expansion of a plane wave into Gaussian beams[J]. Studia Geophysica Et Geodaetica, 1982, 26(2): 120-131. DOI:10.1007/BF01582305
[3]
Hill N R. Gaussian beam migration[J]. Geophysics, 1990, 55(11): 1416-1428. DOI:10.1190/1.1442788
[4]
Hale D.Migration by the Kirchhoff, slant stack, and Gaussian beam methods[R]. CWP Annual Project Review Meeting, Colorado School of Mines, 1992, 22-34.
[5]
Alkhalifah T. Gaussian beam depth migration for ani-sotropic media[J]. Geophysics, 1995, 60(5): 1474-1484.
[6]
Hill N R. Prestack Gaussian-beam depth migration[J]. Geophysics, 2001, 66(4): 1240-1250.
[7]
Nowack R L, Sen M K, Stoffa P L, et al.Gaussian beam migration for sparse common-shot and common-receiver data[C]. SEG Technical Program Expanded Abstracts, 2003, 22: 1114-1117.
[8]
Gray S H. Gaussian beam migration of common shot records[J]. Geophysics, 2005, 70(4): S71-S77.
[9]
Zhu T.Kinematic and dynamic raytracing in aniso-tropic media-a phase-velocity formulation[C]. SEG Technical Program Expanded Abstracts, 2005, 24: 96-99.
[10]
Zhu T, Gray S H, Wang D. Prestack Gaussian-beam depth migration in anisotropic media[J]. Geophysics, 2007, 72(3): S133-S138.
[11]
Gray S H, Bleistein N. True-amplitude Gaussian-beam migration[J]. Geophysics, 2009, 74(2): S11-S23.
[12]
Bleistein N, Zhang Y, Xu S, et al. Migration/inver-sion:think image point coordinates, process in acquisition surface coordinates[J]. Inverse Problems, 2005, 21(5): 1715-1744. DOI:10.1088/0266-5611/21/5/013
[13]
李振春, 岳玉波, 郭朝斌, 等. 高斯波束共角度保幅深度偏移[J]. 石油地球物理勘探, 2010, 45(3): 360-365.
LI Zhenchun, YUE Yubo, GUO Chaobin, et al. Gau-ssian beam common angle preserved-amplitude migration[J]. Oil Geophysical Prospecting, 2010, 45(3): 360-365.
[14]
岳玉波.复杂介质高斯束偏移成像方法研究[D].中国石油大学(华东), 山东青岛, 2011.
YUE Yubo.Study on Gaussian Beam Migration Methods in Complex Medium[D]. China University of Petroleum (East China), Qingdao, Shandong, 2011.
[15]
段鹏飞, 程玖兵, 陈爱萍, 等. TI介质局部角度域高斯束叠前深度偏移成像[J]. 地球物理学报, 2013, 56(12): 4206-4214.
DUAN Pengfei, CHENG Jiubing, CHEN Aiping, et al. Local angle-domain Gaussian beam prestack depth migration in a TI medium[J]. Chinese Journal of Geophysics, 2013, 56(12): 4206-4214. DOI:10.6038/cjg20131223
[16]
高成, 孙建国, 齐鹏, 等. 2D共炮时间域高斯波束偏移[J]. 地球物理学报, 2015, 58(4): 1333-1340.
GAO Cheng, SUN Jianguo, QI Peng, et al. 2-D Gau-ssian-beam migration of common-shot records in time domain[J]. Chinese Journal of Geophysics, 2015, 58(4): 1333-1340.
[17]
杨晶, 黄建平. 基于汉宁窗函数滤波时间域高斯束成像方法[J]. 地球物理学进展, 2018, 33(3): 1161-1166.
YANG Jing, HUANG Jianping. Time domain Gaussian beam migration method based Hanning window filtering[J]. Progress in Geophysics, 2018, 33(3): 1161-1166.
[18]
吴娟, 白敏, 张华. 基于高斯束偏移的鬼波压制[J]. 石油地球物理勘探, 2018, 53(3): 443-448.
WU Juan, BAI Min, ZHANG Hua. Deghosting with Gaussian beam migration[J]. Oil Geophysical Prospecting, 2018, 53(3): 443-448.
[19]
石星辰, 毛伟建, 栗学磊. 矢量黏弹性衰减补偿高斯束偏移[J]. 地球物理学报, 2019, 62(4): 1480-1491.
SHI Xingchen, MAO Weijian, LI Xuelei. Viscoelastic Q-compensated Gaussian beam migration based on vector-wave imaging[J]. Chinese Journal of Geophy-sics, 2019, 62(4): 1480-1491.
[20]
吕考考, 徐基祥, 张才, 等. 数据驱动的控制束偏移方法[J]. 石油地球物理勘探, 2019, 54(4): 805-813.
LYU Kaokao, XU Jixiang, ZHANG Cai, et al. A data-driven controlled-beam migration method[J]. Oil Geophysical Prospecting, 2019, 54(4): 805-813.
[21]
张瑞, 黄建平, 李振春, 等. 相似系数阈值滤波的数据驱动控制束偏移[J]. 石油地球物理勘探, 2019, 54(6): 1267-1279.
ZHANG Rui, HUANG Jianping, LI Zhenchun, et al. A data-driven controlled beam migration based on the semblance threshold filtering[J]. Oil Geophysical Prospecting, 2019, 54(6): 1267-1279.
[22]
Yue Y, Sava P, Qian Z, et al. Least-squares Gaussian beam migration in elastic media[J]. Geophysics, 2019, 84(4): S329-S340.
[23]
Zhang Y, Xu S, Bleistein N, et al. True-amplitude, angle-domain, common-image gathers from one-way wave-equation migrations[J]. Geophysics, 2007, 72(1): S49-S58.
[24]
Aki K, Richards P G.Quantitative Seismology, Theory and Methods[M]. W H Freeman, 1980.
[25]
Červený V. Seismic rays and ray intensities in inhomogeneous anisotropic media[J]. Geophysical Journal International, 1972, 29(1): 1-13. DOI:10.1111/j.1365-246X.1972.tb06147.x
[26]
Hanyga A. Gaussian beams in anisotropic media[J]. Geophysical Journal International, 2007, 85(3): 473-504.