地球物理学报  2020, Vol. 63 Issue (2): 652-665   PDF    
二维弹性多波时空域高斯束偏移方法
胡自多1, 吕庆达2,3, 韩令贺1, 刘威1, 黄建平2, 杨继东2, 李振春2     
1. 中国石油勘探开发研究院西北分院, 兰州 730020;
2. 中国石油大学(华东)地球物理系, 青岛 266555;
3. 中国石油化工股份有限公司石油物探技术研究院, 南京 211003
摘要:随着我国勘探开发难度逐步增大,勘探目标开始向裂缝油气藏、岩性油气藏等复杂探区转移,研究高精度、适应性强的多波多分量深度偏移算法在后续的地震解释、属性分析及储层预测中具有重要意义.针对多波多分量地震数据,本文提出了一种二维弹性波时空域高斯束偏移方法.时空域高斯束沿中心射线传播时能够面向成像目标描述局部波场,且对振幅和频率可调制的Gabor基函数有天然的适应性,因而将基于Gabor分解的子波重构方法应用于震源波场构建,从而得到任意点源函数产生的时空域高斯束波场.该方法由于直接在时间域进行计算,可以避开频率域中出现的假频和边缘截断效应等问题.基于各向同性弹性波动方程的Kirchhoff-Helmholtz积分解,利用矢量时空域高斯束传播算子构建格林函数和格林位移张量,并结合上行射线追踪策略,实现了检波点波场的反向延拓.针对矢量波成像问题,本文借鉴弹性波逆时偏移方法从矢量延拓波场中分离出纯纵波分量和纯横波分量,进而采用修改后的内积成像条件产生具有明确物理意义的PP、PS成像结果,避免了转换波成像的极性反转问题.最后利用简单两层模型和不含盐体构造的部分Sigsbee2a模型的成像结果,并将其与应用近似纵横波成像条件、标量和矢量势成像条件的偏移剖面进行对比,验证了本文方法的正确性和有效性.
关键词: 弹性时空域高斯束      多分量偏移      格林函数      上行射线追踪      内积成像条件     
Elastic space-time Gaussian beam method for seismic depth imaging
HU ZiDuo1, Lü QingDa2,3, HAN LingHe1, LIU Wei1, HUANG JianPing2, YANG JiDong2, LI ZhenChun2     
1. Research Institute of Petroleum Exploration and Development-Northwest PetroChina, Lanzhou 730020, China;
2. Department of Geophysics China University of Petroleum(East China), Qingdao 266555, China;
3. Sinopec Geophysical Research Institute, Nanjing 211103, China
Abstract: As the difficulty of exploration and development in China gradually increases the exploration target begins to transfer to complex exploration areas such as fractured reservoirs and lithologic reservoirs. And it is of vital significance for subsequent seismic interpretation, seismic attribute analysis, and reservoir prediction to develop a multi-wave, multi-component depth migration algorithm which is highly accurate and strongly robust. In order to process common-shot multi-wave and multi-component seismic data, we develop a 2D elastic imaging method using vector space-time Gaussian beam propagation operators. The space-time Gaussian beam can describe the local wave-field of target-oriented areas over the whole central ray path and naturally adapt to the Gabor basis function which is modulated in terms of amplitude and frequency. Therefore, the source wavelet refactored method based on Gabor decomposition is applied to forward wavefield construction and then we can obtain a space-time Gaussian beam wave-field generated by an arbitrary point source function. This method can avoid problems such as aliasing and edge truncation effects occurring in the frequency domain due to its direct calculation in the time domain. Adhering to the framework of the acoustic space-time Gaussian beam method, we perform the up-going ray tracing from subsurface imaging points to the receiver surface for the construction of the backward wavefield. The P-and S-wave-fields are simultaneously extrapolated based on the Kirchhoff-Helmholtz integral solution of the isotropic elastodynamic equation. The Green's functions and Green's displacement tensors are represented with vector space-time Gaussian beam summation. We separate the extrapolated wavefields into compressional and shear modes using the Helmholtz decomposition, and then implement a modified dot-product imaging. This approach permits to produce clear PP images and avoid polarity reversal problems for PS images. Numerical experiments on a simple two-layer model and a partial Sigsbee2a model without salt structure have verified the accuracy and validity of the proposed method.
Keywords: Elastic space-time Gaussian beam    Multi-component migration    Green function    Up-going ray tracing    Dot-product imaging condition    
0 引言

随着地震勘探技术的进步与发展,多波多分量地震数据处理技术逐渐成为了勘探地球物理界研究的热点和难点.与单分量地震数据相比,多波多分量数据含有更为丰富的波场信息,更能代表地震波实际传播的物理规律,因而相比声波偏移,弹性波偏移能够形成多波成像数据,提供更为精确的地下构造信息.目前,弹性波偏移方法分为两类:标量延拓法和矢量延拓法.常规标量法是基于标量波动方程,先将多分量地震波场分离为纵波和横波(Wapenaar et al., 1990; Sun, 1999),后使用现有声波偏移算法成像(Wapenaar and Haimé, 1990; Sun and McMechan, 2001; Han et al., 2013).该方法的最大优势是无需推导弹性介质的矢量波场延拓算子和相应成像条件,但其存在以下两种弊端:(1)该方法难以将相互耦合的不同波型能量彻底分离,尤其是各向异性介质,残余的非本型波能量将在成像结果中产生大量的串扰噪声,导致解释困难(Fomel and Backus, 2003; Nickel and Sonneland, 2004);(2)该方法错误地将垂直分量和水平分量近似为纵波和横波(Yan and Sava, 2008),忽略了多波多分量数据的矢量特性,不能充分利用弹性地震波场的有效信息.

直接使用多分量数据作为边界条件求解弹性波动方程实现弹性波成像的矢量延拓法是解决上述问题的有效途径.Pao和Varatharajulu(1976)首先提出了非均匀介质的弹性Kirchhoff-Helmholtz积分,Kuo和Dai(1984)Dai和Kuo(1986)将其应用于矢量波场延拓,提出各项同性弹性介质的多分量Kirchhoff偏移方法,为射线类矢量波成像奠定了理论基础.Sena和Toksöz(1993)将该方法拓展至各向异性介质中.基于Claerbout提出的沉降观测(survey-sinking)成像条件和黏弹性Kirchhoff积分,Hokstad(2000)实现了黏弹介质多分量Kirchhoff偏移方法.这些方法均通过射线追踪获取本型波和转换波的走时表,以此沿所选波型的走时轨迹叠加矢量数据计算成像振幅.毋庸置疑,射线类方法在研究地震波的正向和反向传播的问题上是一种非常有效且实际应用较为广泛的方法.但对于复杂地质构造,射线类方法受限于零阶射线近似理论,存在焦散区、阴影区和临界区振幅奇异值(Červený and Pšenčík, 1983a)以及多次波至(Gray et al., 2001)等问题.而基于双程波动方程的弹性波逆时偏移采用全波算子,不受成像倾角、波传播路径问题的限制,最初是由Sun和McMechan(1986)提出,随后Chang和McMechan(1987, 1994)利用有限差分算法在深度域实现弹性波逆时偏移,并将其推广到三维空间,为波动类矢量波成像奠定了理论和算法基础.但成像结果受纵横波串扰噪声的影响,且不具有明确的物理意义,无法进行后续解释.针对此问题,常用的一种方法是首先使用耦合的弹性波动方程延拓矢量波场,然后利用散度和旋度算子将其分离为纯P、S波,再应用成像条件得到具有明确物理意义的成像结果(Yan and Sava, 2008; Du et al., 2012, 2014; Duan and Sava, 2015).弹性波逆时偏移虽然具有极高的成像精度,但内存需求量和计算量庞大,且成像结果伴有大量的低频噪声干扰.

基于旁轴射线理论的高斯束是波动方程的高频渐进解(Červený et al., 1982b; Červený, 2001; Popov, 1982, 2002),其利用复值振幅、走时实现动力学射线追踪和波场计算,从而解决了经典射线理论焦散区振幅奇异、阴影区无波场能量等问题(Červený and Pšenčík, 1983ab1984; Červený, 1983).此外,不同出射角度的高斯束独立向下传播实现波场延拓,可获取复杂地质构造的多路径、多走时信息,从而对高陡盐丘构造和复杂断层带准确成像.Hill(1990, 2001)最先提出高斯束偏移的基本框架;此后,高斯束偏移被拓展到共炮域(Nowack et al., 2005; Gray, 2005)、各向异性介质(Alkhalifah, 1995; Zhu et al., 2007; Protasov, 2015; 韩建光等, 2017; 吕庆达等, 2018)、起伏地表(Gray, 2005; 岳玉波等, 2010, 2012; Yang et al., 2014; 黄建平等, 2016)、数据驱动域(Hu and Stoffa, 2009; Yang and Zhu, 2018)以及真振幅偏移(Gray and Bleistein, 2009; Protasov and Tcheverda, 2011).新型地震波束偏移方法也在不断发展并成熟,比如聚焦束偏移(Nowack, 2008, 2009; Yang et al., 2015a)、菲涅尔束偏移(Huang et al., 2015; 杨继东等, 2015)、快速束偏移(Gao et al., 2006; 王旭谦等, 2017)以及复值束偏移(Zhu, 2013)等.时空域高斯束偏移是利用高斯束方法直接在时间空间域描述波束传播,从而实现地震波场的正反向外推,典型方法包括高斯波包偏移、高斯束逆时偏移以及Delta波包偏移.基于Gabor框架函数对观测波场的稀疏分解,高斯波包偏移取得较大进展(Klimeš, 1989; Žáček, 2005; Geng et al., 2014; 李辉等, 2012, 2014).Popov等(2010)从观测波场逆时外推的角度考虑,提出了Kirchhoff积分偏移框架下的高斯束逆时偏移,但是计算量极大;针对该问题,Yang等(2015b)通过引入Gabor变换对震源波场的稀疏分解,减少正向延拓波场计算量,提出了一种高效的时空域高斯束偏移.石秀林等(2016)利用高斯波束的傅立叶逆变换-delta波包的叠加表达时间域格林函数,进而实现了基于delta波包叠加的高精度深度偏移算法;Yang等(2018a)在此基础上将delta波包应用于线性正演,发展了基于高斯束叠加的时间域最小二乘偏移方法,有效提高了成像分辨率.

以上这些针对高斯束偏移方法的研究均用来处理纵波数据,无法直接实现多波多分量数据成像.Huang等(2013)利用矢量高斯束构建震源波场和检波点波场,但是其直接将重构波场的垂直分量和水平分量近似为P波和S波,这在理论上说不通,而且会掺杂串扰噪声.基于弹性波逆时偏移中Helmholtz分解思想,在应用成像条件时利用散度和旋度算子将重构波场分离为纯P波和纯S波分量(Yang et al., 2015c; 栗学磊和毛伟健, 2016),有效避免了反射能量的遗漏,但需通过判断入射波入射角的正、负以引入符号函数对PS转换波成像结果做极性校正.

本文在前人的研究基础上,针对多波多分量地震数据,提出一种二维弹性波时空域高斯束偏移方法.根据波动方程的线性叠加性,任何点源函数均可被分解成一系列Gabor基子波的加权叠加形式,而时空域高斯束适用于可调制的Gabor函数,因而将不同时间延迟的P波型时空域高斯束按权系数线性叠加,从而得到任意点源函数产生的震源波场.基于各向同性弹性波动方程的Kirchhoff-Helmholtz积分解(Pao and Varatharajulu, 1976),利用矢量时空域高斯束传播算子构建格林函数和格林位移张量,并结合上行射线追踪策略,实现检波点波场的反向延拓.针对矢量波成像问题,本文借鉴弹性波逆时偏移方法从矢量延拓波场中分离出纯纵波分量和纯横波分量,进而采用修改后的内积成像条件产生了准确的偏移剖面,且避免了转换波成像的极性反转问题.最后利用简单两层模型和不含盐体构造的部分Sigsbee2a模型的偏移成像结果,证明了本文所提出方法的准确性和优越性.

1 理论方法

弹性波时空域高斯束偏移包括两个步骤:(1)构建震源波场和检波点波场;(2)应用成像条件提取成像值.

1.1 矢量时空域高斯束

利用时空射线理论可直接寻找广义点源弹性动力学方程的渐近解(Červený et al., 1982a).常规射线理论均是基于弹性动力学方程的最为简单的表达形式,描述的是高频信号沿着射线路径的传播问题.但当其接触到焦散区(近似信号的Hilbert变换形式)或传播到界面(近似信号的线性叠加或Hilbert变换形式),信号形式通常会发生改变.因而,当计算广义点源产生的非稳波场时,需要考虑时空射线理论.Katchalov和Popov(1988)推导出沿中心时空射线传播的矢量时空域高斯束的表达式为

(1)

式中,UP(s, n, t-τ)和U S(s, n, t-τ)分别代表纯纵波震源和纯横波震源激发产生的位移矢量场;(s, n)为射线中心坐标;τ是中心时空射线走时,通过动力学射线追踪获取;t代表时间;A(t)和θ(t)均为充分光滑的函数(A(t)在有限时间区间内不为零),且有θ′=dθ/dt≠0, θ′(t)≤const.<0,使得当ξ→∞,瞬时频率ω=-ξθ′(t)满足高频假设;θ0为严格满足中心射线程函方程、输运方程的θ值,θ′0=dθ0/dtα(s)、β(s)和ρ(s)分别为P波速度、S波速度以及密度值;P(s)和Q(s)是动力学射线追踪(Popov, 2002; Červený, 2001)的复值基本解;ev(s)=(dr(s)/ds)为相应极化矢量,表示为

(2)

图 1所示,P波主分量为t,次分量为n;S波主分量为n,次分量为- t (Červený and Pšenčík, 1983b).假设点源函数f(t)的振幅和频率为可调制函数,考虑Gabor基子波,公式为

(3)

图 1 主分量与次分量示意图 (a) P波;(b) S波. Fig. 1 Principal and additional components (a) P-wave; (b) S-wave.

将式(2)、(3)两式代入到(1)式,可得与其相匹配的时空域高斯束为

(4)

根据Katchalov和Popov(1990),利用矢量时空域高斯束,可将由(s0, 0)(s0代表初始弧长,通常设为0)处全方位不同波震源引起的(s, n)处的位移矢量场通过一系列以(s0, 0)为初始点,具有不同出射角且对(s, n)有贡献的时空域高斯束的叠加来表示,即:

(5)

其中,ψP(φ)和ψS(φ)分别为P型和S型时空域高斯束初始振幅,通过对比均匀介质中位移矢量的解析解与式(5)的高频渐进解求得:

(6)

式中,i是虚数单位;ε是复值波束参数(Červený et al., 1982b; Červený and Pšenčík, 1983b; Müller, 1984),有:

(7)

向量[q1(s), p1(s)]和[q2(s), p2(s)]分别为动力学射线追踪的初始平面波解和初始点源解.

1.2 重构子波构建震源波场

假设在二维各向同性完全弹性介质中,震源为全方位的P波线源,则此时介质中仅存在P波和转换S波.根据李辉等(2012),任何点源函数均可被分解成一系列Gabor基子波gi(t-τi)的线性叠加形式,即:

(8)

式中,Re为取实部符号;τk为第k个Gabor基子波函数的时间延迟;ak为复值权系数,通过求解如下复值线性方程组得到:

(9)

其中,(gi, gj)和[gi, f(t)]表示函数内积.

图 2为原始雷克子波与重构雷克子波,主频为30 Hz.由图可见,只要选取合适的Gabor基函数和子波参数,就可通过式(9)线性叠加来精确重构任意点源函数f(t).将不同时间延迟的P波型时空域高斯束按权系数线性叠加,可得到任意点源函数产生的时空域高斯束波场.进而结合式(4)、(5)、(6)三式,推导可得到任意点源函数在炮点x s处产生的正向波场为

(10)

图 2 30 Hz雷克子波分解示意图 (a)原始子波和重构子波;(b)分解的Gabor基函数. Fig. 2 Sketch of Gabor decomposition for Ricker wavelet (a) The original Riker wavelet; (b) Reconstructed wavelet, lower is the Gabor base function.

与弹性波频率域高斯束相比,该方法由于直接在时间域进行计算,可以避开频率域中出现的假频和边缘截断效应等问题.此外,借助于初始子波的傅里叶变换形式,也可构造正向延拓波场(Popov et al., 2010; Huang et al., 2017).

1.3 弹性波波场反向延拓

本小节基于上行射线追踪策略(图 3)推导弹性波时空域高斯束反向波场延拓算子.弹性波Kirchhoff-Helmholtz积分方程相当于弹性动力学方程的积分表示形式(Pao and Varatharajulu, 1976; Aki and Richards, 2002).假设S为自由地表,在自由应力边界条件下,忽略体力项,t时刻检波点xr处的多分量地震波场u (x r, t0)反向延拓的在成像点x处的弹性位移矢量场可以表示为

(11)

图 3 弹性波时空域高斯束偏移实现示意图 Fig. 3 Scheme of elastic space-time Gaussian migration

这里符号T表示转置;m是边界S的外法向向量;Σ (x, xr)为格林应力张量,在各向同性介质中可表示为

(12)

其中,lmnk为笛卡尔坐标分量;δ是Kronecker Delta函数;λμ是弹性拉梅系数;G (x, t; xr)表示格林位移张量,具有高频近似形式,即:

(13)

式中,gP(x, t; xr, t0)和g S(x, t; x r, t0)分别是式(5)中g P(s, n, t)和g S(s, n, t)由射线中心坐标转换到笛卡尔坐标得到,xrx分别对应于(s0, 0)和(s, n),即:

(14)

κP(xr)和κS(xr)表示检波点xr处的应力方向,即:

(15)

其中φ是式(5)中的出射角.

式(12)中Gm, n, l表示∂Gxm, xn/∂xl,是格林位移张量的空间偏导数,在高频渐进条件下可近似为其主阶项(Popov et al., 2010),即:

(16)

式中,pl是慢度向量的xl分量.考虑m=(0, 0, -1),即水平记录地表,将式(12)—(15)代入式(11)中,推导可得到反向延拓的矢量波场为

(17)

其中:

(18)

Φ为二阶2×2系数矩阵,其分量展写为

(19)

式中,γ=β(xr)/α(xr)是检波点xr处的纵横波速度比;pxPpzP分别为P波高斯束在中心射线的水平、垂直慢度;pxSpzS分别为S波高斯束在中心射线的水平、垂直慢度.

1.4 成像条件

弹性波场为矢量波场,常规的互相关成像条件不适用于矢量波成像,若直接将各分量成像,会因为多种混合的波型产生串扰假象,且难以给出明确的物理含义.Yan和Sava(2008)利用散度、旋度算子从矢量波场中分离出纯纵波分量和纯横波分量,并提出标量和矢量势成像条件,但是该方法存在PS转换波成像极性反转问题(图 4),因而常规弹性波高斯束偏移需要通过判断入射波入射角的正、负以引入符号函数对PS转换波成像结果做极性校正(Yang et al., 2015c; 栗学磊等, 2016; 李振春等, 2018).为解决矢量波成像问题,本文拟采用Zhu(2017)提出的内积成像条件,产生具有准确相位、振幅的PP、PS偏移剖面,公式为

(20)

图 4 反射波极性分析 (a) PP;(b) PS.表示偏振方向垂直向里,表示偏振方向垂直向外. Fig. 4 Polarity analysis of PP and PS reflection (a) PP; (b) PS.denotes inward polarity direction.denotes outward polarity direction.

其中“·”表示内积;usP(x, t)、urP(x, t)、urS(x, t)分别为

(21)

式中,uP(x, t)和uS(x, t)为纯纵波分量和纯横波分量;下标s和r分别表示炮点和检波点;uD(x, t; xs)和uB(x, t; xr)分别为利用(10)和(17)两式计算的成像点x处的正向、反向延拓波场;▽ (▽ ·)和▽ × ▽ ×是用于提取纯纵波、纯横波分量的算子;t∈[T1, T2]为相关时窗,该参数可以为特定的问题做具体选择,比如在面向目标偏移中其取值不应超过震源波场穿过目标区域的延续时间.

该成像条件虽然解决了波场分解的相位和振幅失真问题,避免了PS转换波成像的极性反转问题,但其涉及到一个余弦函数cosΔφ的乘积,导致反射界面的错误估测(Yang et al., 2018b).内积,从代数上来讲是两个向量对应分量乘积之和,从几何上来讲是两个向量的模和两者夹角余弦这三者的乘积.据此,式(20)可简写为

(22)

其中“||”表示矢量场的模;ΔφPP(x)是usP(x, t)与urP(x, t)的夹角,ΔφPS(x)是usP(x, t)与urP(x, t)的夹角,如图 3所示,取值为(Du et al., 2017):

(23)

式中,φsP(x)为下行P波时空域高斯束在成像点x处的旁轴射线传播角度;φrP(x)为上行P波时空域高斯束在成像点x处的旁轴射线传播角度;φrS(x)为上行S波时空域高斯束在成像点x处的旁轴射线传播角度.

因此,为消除cosΔφ对最终偏移剖面产生的振幅影响和常规弹性波高斯束PS转换波成像极性反转问题,并结合式(20)—(23),这里直接给出修改后的内积成像条件为

(24)

其中,ε2是为保证运算稳定性加到分母上的较小值.

1.5 算法流程

基于时空射线理论形成的时空域高斯束理论与常规的频率-空间域高斯束理论的差别在于引入时域渐进射线级数和定义时空射线.如图 5所示,时空射线在空间上的投影即常规射线(渐进射线或旁轴射线),其中实线表示时空射线,虚线表示常规射线.但在实际应用中,我们无需具体计算时空射线,其时间坐标可根据常规射线的走时来定义,所以只需通过运动学射线追踪计算常规的空间射线即可.因此,我们设计了如图 6所示的共炮域偏移算法流程,包括五重循环(炮点、成像点、炮点高斯束、成像点高斯束和检波点波场),成像过程清晰明了;炮点下行射线追踪和成像点上行射线追踪通过经典四阶Runge-Kutta数值算法求解运动学和动力学射线追踪方程组实现,从而获取下行高斯束和上行高斯束的走时、振幅以及动力学参量(式(7));当上行射线到达地表的终点不与检波点重合时,我们根据Hale(1992)中提及的最近点搜索算法判断距离检波点最近的射线时间步,从而构建弹性波时空域高斯束反向波场延拓算子.

图 5 时空射线与常规射线的映射关系(Červený et al., 1982a) 实线表示时空射线,虚线表示常规射线. Fig. 5 The projection relation between space-time rays and space rays (Červený et al., 1982a) Solid lines denote space-time rays and dotted lines denote general rays.
图 6 弹性波时空域高斯束偏移算法流程示意图 Fig. 6 Algorithm flow chart of elastic space-time Gaussian beam migration
2 数值算例

本文利用2D模型数据的偏移成像结果,来验证所提出方法的准确性.为简单起见,文中模型通过取纵横波速比为设定横波波速.模型记录合成采用弹性波交错网格有限差分方法计算.

2.1 两层模型

两层介质速度模型如图 7所示,大小为3 km×3 km,反射界面在深度z=1500 m.震源为20 Hz雷克子波,坐标(1500 m, 0);检波器放置于地表 0 m到3 km处,道间距10 m,接收记录长度2.0 s,时间采样间隔0.5 ms.多分量地震记录如图 8所示,(a)为z分量,(b)为x分量.x分量中,转换横波能量比反射纵波能量强,且由于纵波入射角正负符号的不同,出现同相轴极性反转现象,零度角入射时不产生转换横波(图 8a);z分量中,反射纵波能量比转换横波能量强,但不存在极性反转现象(图 8b).

图 7 两层介质速度模型 Fig. 7 Two-layer medium velocity model
图 8 两层模型多分量正演炮记录 (a) x分量;(b) z分量.黑色标识箭头标识反射PP波和转换PS波. Fig. 8 The common-shot multicomponent records for the two-layer model (a) x components; (b) z components.The black arrows indicate the reflected PP-waves and converted PS-waves.

为验证提出方法的有效性,本文利用Huang等(2013)中近似纵横波成像条件、栗学磊和毛伟健(2016)中标量和矢量势的成像条件及(24)式进行对比测试,试算结果如图 9所示.由图 9可以看出,Huang等(2013)使用反向延拓波场的zx分量来近似代替P、S波获得PP、PS偏移剖面(图 9ab),但这种近似不仅遗漏了P波x分量、S波z分量的有效成像能量,还具有纵横波未分离彻底产生的串扰噪声,如图 9a红色圆框、9b红色箭头所示.采用标量和矢量势成像条件和本文提出的式(24)均能够将PP、PS波完全分离,具有很好的串扰压制效果.此外,未做极性校正的PS成像,在纵波垂直入射点位置,附近能量消失,同相轴断开,垂直入射点两侧同相轴极性相反(见图 9bd黑色箭头).本文使用前述(24)式成像条件得到的成像剖面清晰,有效地避免了PS转换波成像的极性反转问题(图 9ef).

图 9 两层介质模型应用不同成像条件的PP、PS偏移结果对比 (a—b)近似纵横波成像条件;(c—d)标量和矢量势成像条件;(e—f)本文式(24)成像条件. Fig. 9 Comparison of PP and PS imaging results for two-layer model using three imaging conditions (a—b) Huang et al.(2013)′s method; (c—d) Scalar and vector potential imaging condition; (e—f) Equation (24) (this paper).
2.2 Sigsbee2a模型

为进一步验证本文方法的有效性,在此对不包含盐体构造的部分Sigsbee2a模型进行测试,其纵波速度如图 10a所示,横波速度通过建立,由于动力学射线追踪需要速度的二阶偏导,因而实现偏移时应用Seismic Unix软件中smooth2程序对速度场做了十点平滑(图 10b).该模型网格大小为801×601,纵横向间隔均为8 m,具有典型正断层和绕射点,具有一定的构造代表性.弹性波多分量数据集使用交错网格有限差分正演模拟获得,该数据集共9炮,炮间隔800 m,每炮801道,全排列接收,道间距8 m,接收记录长度8.0 s,时间采样间隔0.5 ms.

图 10 部分Sigsbee2a速度模型(P波) (a)真实速度;(b)偏移速度. Fig. 10 The partial Sigsbee2a model(P wave) (a) The true velocities; (b) The migration velocities.

图 11abcdef分别为利用Huang等(2013)中近似纵横波成像条件、栗学磊和毛伟健(2016)中标量和矢量势的成像条件及前述(24)式成像条件对该数据集进行处理得到的弹性波偏移结果,图 12为应用不同成像条件的PS偏移结果在两块区域范围内的局部放大对比.从成像剖面看出,三种方法的PP成像均能够对典型的正断层和绕射点清晰地成像,具有较高的分辨率.但使用近似纵横波成像条件偏移的剖面(图 11a)显得更为嘈杂,这是由于波场未分离引起的串扰造成的.而由于转换波成像极性反转问题的存在,如图 11b图 12ad图 11d图 12be所示,同相轴不连续,多炮成像叠加剖面严重受损,成像精度大幅度下降.应用本文提出的成像条件,PS转换波成像结果清晰,分辨率明显提高,同相轴连续无错断(图 11f图 12cf).

图 11 Sigsbee2a模型应用不同成像条件的PP、PS偏移结果对比 (a—b)近似纵横波成像条件;(c—d)标量和矢量势成像条件;(e—f)本文式(24)成像条件. Fig. 11 Comparison of PP and PS imaging results for Sigsbee2a model using the three imaging conditions (a—b) Huang et al.(2013)′s method; (c—d) scalar and vector potential imaging condition; (e—f) Equation (24) (this paper).
图 12 Sigsbee2a模型应用不同成像条件的PS偏移结果局部放大对比 (a, d)近似纵横波成像条件;(b, e)标量和矢量势成像条件;(c, f)本文式(24)成像条件. Fig. 12 Comparison of locally zoomed PS migration results for Sigsbee2a model using three imaging conditions (a, d) Huang et al.(2013)′s method; (b, e) Scalar and vector potential imaging condition; (c, f) Equation (24) (this paper).
3 结论与讨论

针对多波多分量地震数据,本文提出了一种二维弹性波时空域高斯束偏移方法.时空域高斯束对振幅和频率可调制的Gabor基函数有天然的适应性,而任何点源函数均可被分解为一系列Gabor基子波的加权叠加形式,因而将不同时间延迟的P波型时空域高斯束按权系数线性叠加,可得到任意点源函数产生的震源波场.基于弹性波Kirchhoff-Helmholtz积分方程,利用矢量时空域高斯束传播算子构建格林函数和格林位移张量,从而实现检波点波场的反向延拓.针对转换波成像极性反转问题,本文采用修改后的内积成像条件产生具有准确相位、振幅的PP、PS偏移剖面.模型测试结果表明:同Huang等(2013)中近似纵横波成像条件、栗学磊和毛伟健(2016)中标量和矢量势的成像条件相比,本文方法获得的PP、PS成像剖面信噪比得到极大提高,且具有准确的相位和振幅.由此可见,二维弹性多波时空域高斯束偏移方法是一种高精度、适应性强的有效地震叠前深度成像工具之一.

共炮点时空域弹性高斯束偏移算法包括炮点、成像点、炮点高斯束、成像点高斯束以及检波点波场等五重循环,如图 6所示,这与频率域高斯束偏移在粗网格上计算成像结果不同,时空域高斯束偏移在实现时要在细网格上计算单炮偏移孔径内所有的成像点.这意味着时空域高斯束偏移将比频率-空间域高斯束偏移需要更多的计算机资源.就本文来讲,考虑到时空域偏移不对输入数据进行局部倾斜叠加处理的特点,所以直接移植了经典Kirchhoff型偏移以及Popov在2010年提出的严格遵守Kirchhoff积分理论的高斯束叠加偏移的实现方式,即单道互相关叠加成像.虽然这样做的好处是只需要在内存中保留模型,但是需要面对时间域褶积以及大量使用I/O所带来的计算成本.对于小尺度问题,这些问题可能表现得并不明显,但是对于大尺度问题,这些问题可能对时空域偏移的实用化形成瓶颈.因此,需要从整体成像算法框架的角度对时空域偏移算法进行顶层设计.这将是作者下一步的工作方向.

致谢  感谢审稿专家对本文给出的指导和建议!
References
Aki K, Richards P G. 2002. Quantitative Seismology. 2nd ed. Herndon: University Science Books.
Alkhalifah T. 1995. Gaussian beam depth migration for anisotropic media. Geophysics, 60(5): 1474-1484. DOI:10.1190/1.1443881
Červený V. 1983. Synthetic body wave seismograms for laterally varying layered structures by the Gaussian beam method. Geophysical Journal International, 73(2): 389-426. DOI:10.1111/j.1365-246X.1983.tb03322.x
Červený V. 2001. Seismic Ray Theory. Cambridge: Cambridge University Press.
Červený V, Molotkov I A, Pšenčík L, et al. 1982a. Space-time ray method for seismic wave fields. Studia Geophysica et Geodaetica, 26(4): 342-351. DOI:10.1007/BF01639635
Červený V, Popov M M, Pšenčík I. 1982b. Computation of wave fields in inhomogeneous media-Gaussian beam approach. Geophysical Journal International, 70(1): 109-128. DOI:10.1111/j.1365-246X.1982.tb06394.x
Červený V, Pšenčík L. 1983a. Gaussian beams and paraxial ray approximation in three-dimensional elastic inhomogeneous media. Journal of Geophysics, 53(1): 1-15.
Červený V, Pšenčík L. 1983b. Gaussian beams in two-dimensional elastic inhomogeneous media. Geophysical Journal International, 72(2): 417-433. DOI:10.1111/j.1365-246X.1983.tb03793.x
Červený V, Pšenčík L. 1984. Gaussian beams in elastic 2-D laterally varying layered structures. Geophysical Journal International, 78(1): 65-91. DOI:10.1111/j.1365-246X.1984.tb06472.x
Chang W F, McMechan G A. 1987. Elastic reverse-time migration. Geophysics, 52(10): 1365-1375. DOI:10.1190/1.1442249
Chang W F, McMechan G A. 1994. 3-D elastic prestack reverse-time depth migration. Geophysics, 59(4): 597-609. DOI:10.1190/1.1443620
Dai T F, Kuo J T. 1986. Real data results of Kirchhoff elastic wave migration. Geophysics, 51(4): 1006-1011. DOI:10.1190/1.1442139
Du Q Z, Gong X F, Zhang M Q, et al. 2014. 3D PS-wave imaging with elastic reverse-time migration. Geophysics, 79(5): S173-S184. DOI:10.1190/geo2013-0253.1
Du Q Z, Guo C F, Zhao Q, et al. 2017. Vector-based elastic reverse time migration based on scalar imaging condition. Geophysics, 82(2): S111-S127. DOI:10.1190/geo2016-0146.1
Du Q Z, Zhu Y T, Ba J. 2012. Polarity reversal correction for elastic reverse time migration. Geophysics, 77(2): S31-S41. DOI:10.1190/geo2011-0348.1
Duan Y T, Sava P. 2015. Scalar imaging condition for elastic reverse time migration. Geophysics, 80(4): S127-S136. DOI:10.1190/geo2014-0453.1
Fomel S, Backus M M. 2003. Multicomponent seismic data registration by least squares. //73rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 781-784.
Gao F C, Zhang P, Wang B, et al. 2006. Fast beam migration-a step toward interactive imaging. //74th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2470-2471.
Geng Y, Wu R S, Gao J H. 2014. Gabor-frame-based Gaussian packet migration. Geophysical Prospecting, 62(6): 1432-1452. DOI:10.1111/1365-2478.12152
Gray S H. 2005. Gaussian beam migration of common-shot records. Geophysics, 70(4): S71-S77. DOI:10.1190/1.1988186
Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23. DOI:10.1190/1.3052116
Gray S H, Etgen J, Dellinger J, et al. 2001. Seismic migration problems and solutions. Geophysics, 66(5): 1622-1640. DOI:10.1190/1.1487107
Hale D. 1992. Computational Aspects of Gaussian Beam Migration. Washington, DC: Colorado School of Mines, Center for Wave Phenomena Report: 139.
Han J G, Wang Y, Chen C X, et al. 2017. Multiwave Gaussian beam prestack depth migration for 2D anisotropic media. Progress in Geophysics, 32(3): 1130-1139.
Han J G, Wang Y, Lu J. 2013. Multi-component Gaussian beam prestack depth migration. Journal of Geophysics and Engineering, 10(5): 055008. DOI:10.1088/1742-2132/10/5/055008
Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428. DOI:10.1190/1.1442788
Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics, 66(4): 1240-1250. DOI:10.1190/1.1487071
Hokstad K. 2000. Multicomponent Kirchhoff migration. Geophysics, 65(3): 861-873. DOI:10.1190/1.1444783
Hu C S, Stoffa P L. 2009. Slowness-driven Gaussian-beam prestack depth migration for low-fold seismic data. Geophysics, 74(6): WCA35-WCA45. DOI:10.1190/1.3250268
Huang J, Duan X, Yuan M, et al. 2013. Elastic Gaussian beam migration method for rugged topography. //75th EAGE Conference & Exhibition, Tu P06.
Huang J P, Yang J D, Li Z C, et al. 2016. An amplitude-preserved Gaussian beam migration based on wave field approximation in effective vicinity under irregular topographical conditions. Chinese Journal of Geophysics (in Chinese), 59(6): 2245-2256.
Huang J P, Yang J D, Liao W Y, et al. 2015. Common-shot Fresnel beam migration based on wave-field approximation in effective vicinity under complex topographic conditions. Geophysical Prospecting, 264(3): 554-570.
Huang J P, Yuan M L, Zhang Q, et al. 2017. Reverse time migration with elastodynamic Gaussian beams. Journal of Earth Science, 28(4): 695-702. DOI:10.1007/s12583-015-0609-9
Katchalov A P, Popov M M. 1988. Gaussian beam methods and theoretical seismograms. Geophysical Journal International, 93(3): 465-475. DOI:10.1111/j.1365-246X.1988.tb03874.x
Katchalov A P, Popov M M. 1990. Application of the method of summation of Gaussian beams to the calculation of theoretical seismograms. Journal of Soviet Mathematics, 50(4): 1727-1742. DOI:10.1007/BF01097103
Klimeš L. 1989. Gaussian packets in the computation of seismic wavefields. Geophysical Journal International, 99(2): 421-433. DOI:10.1111/j.1365-246X.1989.tb01699.x
Kuo J T, Dai Y. 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver. Geophysics, 49(8): 1223-1238. DOI:10.1190/1.1441751
Li H, Feng B, Wang H Z. 2012. A new wave field modeling method by using Gaussian Packets. Geophysical Prospecting for Petroleum (in Chinese), 51(4): 327-337.
Li H, Wang H Z, Feng B, et al. 2014. Efficient pre-stack depth migration method using characteristic Gaussian packets. Chinese Journal of Geophysics (in Chinese), 57(7): 2258-2268. DOI:10.6038/cjg20140720
Li X L, Mao W J. 2016. Multimode and multicomponent Gaussian beam prestack depth migration. Chinese Journal of Geophysics (in Chinese), 59(8): 2989-3005.
Li Z C, Liu Q, Han W G, et al. 2018. Angle domain converted wave Gaussian beam migration in VTI media. Chinese Journal of Geophysics (in Chinese), 61(4): 1471-1481.
Lü Q D, Huang J P, Li Z C, et al. 2018. Efficient Gaussian beam method in time domain for anisotropic media. Progress in Geophysics, 33(6): 2428-2434.
Müller G. 1984. Efficient calculation of Gaussian-beam seismograms for two-dimensional inhomogeneous media. Geophysical Journal International, 79(1): 153-166. DOI:10.1111/j.1365-246X.1984.tb02847.x
Nickel M, Sonneland L. 2004. Automated PS to PP event registration and estimation of a high-resolution VP-VS ratio volume. //74th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 869-872.
Nowack R L. 2008. Focused gaussian beams for seismic imaging. //78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2376-2380.
Nowack R L. 2011. Dynamically focused gaussian beams for seismic imaging. International Journal of Geophysics, 2011: 316581.
Nowack R L, Sen M K, Stoffa P L. 2005. Gaussian beam migration for sparse common-shot and common-receiver data. //75th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1114-1117.
Pao Y H, Varatharajulu V. 1976. Huygens' principle, radiation conditions, and integral formulas for the scattering of elastic wave. The Journal of the Acoustical Society of America, 59(6): 1361-1371. DOI:10.1121/1.381022
Popov M M. 1982. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4(1): 85-97. DOI:10.1016/0165-2125(82)90016-6
Popov M M. 2002. Ray Theory and Gaussian Beam Method for Geophysicists, Lecture Notes University of Bahia. Salvador, Brazil:Universidade Federal da Bahia, Salvador, Brazil.
Popov M M, Semtchenok N M, Popov P M, et al. 2010. Depth migration by the Gaussian beam summation method. Geophysics, 75(2): S81-S93. DOI:10.1190/1.3361651
Protasov M I. 2015. 2-D Gaussian beam imaging of multicomponent seismic data in anisotropic media. Geophysical Journal International, 203(3): 2021-2031. DOI:10.1093/gji/ggv408
Protasov M I, Tcheverda V A. 2011. True amplitude imaging by inverse generalized Radon transform based on Gaussian beam decomposition of the acoustic Green's function. Geophysical Prospecting, 59(2): 197-209. DOI:10.1111/j.1365-2478.2010.00920.x
Sena A G, Toksöz M N. 1993. Kirchhoff migration and velocity analysis for converted and nonconverted waves in anisotropic media. Geophysics, 58(2): 265-276. DOI:10.1190/1.1443411
Shi X L, Sun J G, Sun H, et al. 2016. The time-domain depth migration by the summation of delta packets. Chinese Journal of Geophysics (in Chinese), 59(7): 2641-2649.
Sun R. 1999. Separating P- and S-waves in a prestack 2-Dimensional elastic seismogram. //61st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 6023.
Sun R, McMechan G A. 1986. Pre-stack reverse-time migration for elastic waves with application to synthetic offset vertical seismic profiles. Proceedings of the IEEE, 74(3): 457-465. DOI:10.1109/PROC.1986.13486
Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics, 66(5): 1519-1527. DOI:10.1190/1.1487098
Wang X Q, Yang J D, Huang J P, et al. 2017. Fast Gaussian beam migration method based on matching-pursuit decomposition and events-picking in tau-p domain. Progress in Geophysics, 32(2): 745-752.
Wapenaar C P A, Haimé G C. 1990. Elastic extrapolation of primary seismic P-and S-waves. Geophysical Prospecting, 38(1): 23-60. DOI:10.1111/j.1365-2478.1990.tb01833.x
Wapenaar C P A, Herrmann P, Verschuur D J, et al. 1990. Decomposition of multicomponent seismic data into primary P-and S-wave responses. Geophysical Prospecting, 38(6): 633-661. DOI:10.1111/j.1365-2478.1990.tb01867.x
Yan J, Sava P. 2008. Isotropic angle-domain elastic reverse-time migration. Geophysics, 73(6): S229-S239. DOI:10.1190/1.2981241
Yang J D, Huang J P, Wang X, et al. 2014. Amplitude-preserved Gaussian beam migration based on wave field approximation in effective vicinity under rugged topography condition. //84th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3852-3856.
Yang J D, Huang J P, Wang X, et al. 2015. Prestack Fresnel beam migration method under complex topographic conditions. Chinese Journal of Geophysics (in Chinese), 58(10): 3758-3770.
Yang J D, Huang J P, Wang X, et al. 2015a. An amplitude-preserved adaptive focused beam seismic migration method. Petroleum Science, 12(3): 417-427. DOI:10.1007/s12182-015-0044-7
Yang J D, Huang J P, Wang X, et al. 2015b. Prestack depth migration method using the time-space Gaussian beam. //85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4303-4307.
Yang J D, Huang J P, Wang X, et al. 2015c. Common-shot elastic Gaussian beam depth migration. //85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2159-2164.
Yang J D, Zhu H J. 2018. A practical data-driven optimization strategy for Gaussian beam migration. Geophysics, 83(1): S81-S92. DOI:10.1190/geo2017-0314.1
Yang J D, Zhu H J, McMechan G, et al. 2018a. Time-domain least-squares migration using the Gaussian beam summation method. Geophysical Journal International, 214(1): 548-572. DOI:10.1093/gji/ggy142
Yang J D, Zhu H J, Huang J P, et al. 2018b. 2D isotropic elastic Gaussian-beam migration for common-shot multicomponent records. Geophysics, 83(2): S127-S140. DOI:10.1190/geo2017-0078.1
Yue Y B, Li Z C, Qian Z P, et al. 2012. Amplitude-preserved Gaussian beam migration under complex topographic conditions. Chinese Journal of Geophysics (in Chinese), 55(4): 1376-1383.
Yue Y B, Li Z C, Zhang P, et al. 2010. Prestack Gaussian beam depth migration under complex surface conditions. Applied Geophysics (in Chinese), 7(2): 143-148. DOI:10.1007/s11770-010-0238-0
Žáček K. 2005. Gaussian packet pre-stack depth migration of the Marmousi data set. //75th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1822-1825.
Zhu H J. 2017. Elastic wavefield separation based on the Helmholtz decomposition. Geophysics, 82(2): S173-S183. DOI:10.1190/geo2016-0419.1
Zhu T F. 2013. Complex-beam migration: formulation and comparisons. //83rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3687-3691.
Zhu T F, Gray S H, Wang D L. 2007. Prestack Gaussian-beam depth migration in anisotropic media. Geophysics, 72(3): S133-S138. DOI:10.1190/1.2711423
韩建光, 王赟, 陈传绪, 等. 2017. 二维各向异性多波高斯束叠前深度偏移. 地球物理学进展, 32(3): 1130-1139.
黄建平, 杨继东, 李振春, 等. 2016. 基于有效邻域波场近似的起伏地表保幅高斯束偏移. 地球物理学报, 59(6): 2245-2256.
李辉, 冯波, 王华忠. 2012. 波场模拟的高斯波包叠加方法. 石油物探, 51(4): 327-337. DOI:10.3969/j.issn.1000-1441.2012.04.002
李辉, 王华忠, 冯波, 等. 2014. 特征高斯波包叠前深度偏移方法. 地球物理学报, 57(7): 2258-2268.
栗学磊, 毛伟健. 2016. 多波多分量高斯束叠前深度偏移. 地球物理学报, 59(8): 2989-3005.
李振春, 刘强, 韩文功, 等. 2018. VTI介质角度域转换波高斯束偏移成像方法研究. 地球物理学报, 61(4): 1471-1481.
吕庆达, 黄建平, 李振春, 等. 2018. VTI介质时间域高斯束偏移. 地球物理学进展, 33(6): 2428-2434.
石秀林, 孙建国, 孙辉, 等. 2016. 基于delta波包叠加的时间域深度偏移. 地球物理学报, 59(7): 2641-2649.
王旭谦, 杨继东, 黄建平, 等. 2017. 基于匹配追踪分解和tau-p域拾取的快速高斯束偏移方法. 地球物理学进展, 32(2): 745-752.
杨继东, 黄建平, 王欣, 等. 2015. 复杂地表条件下叠前菲涅尔束偏移方法. 地球物理学报, 58(10): 3758-3770. DOI:10.6038/cjg20151026
岳玉波, 李振春, 钱忠平, 等. 2012. 复杂地表条件下保幅高斯束偏移. 地球物理学报, 55(4): 1376-1383. DOI:10.6038/j.issn.0001-5733.2012.04.033
岳玉波, 李振春, 张平, 等. 2010. 复杂地表条件下高斯波束叠前深度偏移. 应用地球物理, 7(2): 143-148.