基于Born散射理论的二维黏声介质高斯波束正演

孙辉, 岳玉波, 李猛

1 西南交通大学地球科学与环境工程学院， 成都 611756 2 西南交通大学高速铁路线路工程教育部重点实验室， 成都 610031 3 东方地球物理公司物探技术研究中心， 涿州 072751

0 引言

地震波场数值模拟是勘探地震学重要组成部分，可以揭示地震波的传播规律，有效指导地震资料的采集、处理和解释，同时还是多种地震反演方法的基础(Long et al., 2019; Wang and Zhan, 2020).传统地震正演方法通常忽略介质黏滞性对地震波传播的影响，但这一性质会显著影响地震波的动力学特征，包括能量衰减、频带变窄、相位畸变等(Walcott, 1970; White, 1975；Zong et al., 2020).因此研究快速、高精度黏声介质地震正演方法具有重要意义.

地震波场数值模拟方法主要分为两大类.一类是基于波动方程的数值解法，这类方法虽然计算效率不高但具有较高的计算精度，且可以计算全波场，因此众多学者针对这一类方法在黏声介质中的应用展开了研究(Arntsen et al., 1998; Operto et al., 2007; Abubakar and Habashy, 2013; 吴玉等，2017；Li et al., 2019; Wang et al., 2020).另一类方法是射线类地震波场数值模拟方法，这一类方法基于射线理论，具有较高的计算效率，但传统射线类正演方法存在计算阴影区，在焦散区无法获得准确的振幅并且难以处理多至走时问题(Jaramillo and Bleistein, 1999; Duquet et al., 2000; Fawcett, 2001; Chapman, 2004; Shekar and Tsvankin, 2014).

虽然射线Born正演相对于传统射线类地震正演方法的计算精度更高(Bleistein et al., 2001; 符立耘，2010)，但由于该方法采用射线类方法计算格林函数，仍然无法从根本上解决射线类正演方法所存在的问题.作为解决办法，可以使用基于高斯束叠加的格林函数代入Born正演公式中.由高斯束叠加获取的格林函数可以处理多路径及焦散问题，其有效性已在相关偏移成像方法中得到了验证(Hill, 1990, 2001; Gray and Bleistein, 2009; Liu and Palacharla, 2010; 曹文俊等；2013；杨继东等，2016；韩建光等，2017；杨晶和黄建平，2018；Sun et al.,2020).Shekar和Tsvankin(2014)将基于高斯束叠加的格林函数和Kirchhoff散射理论相结合提出了一种针对各向异性衰减介质的地震波正演方法.Huang等(2016)以及岳玉波等(2019)将由高斯束叠加获取的格林函数融入了Born散射理论中，提出了一种针对声波介质的地震波正演方法，并取得了良好的数值模拟结果.

本文提出了一种针对黏声介质的二维高斯波束Born正演方法，在Born正演理论的框架下采用由高斯束叠加获得的格林函数，在高斯束传播算子中考虑了介质黏滞性信息，并将其融入了wavelet-bank平面波合成方法中.在接下来的章节中，首先推导出基于高斯束的Born正演公式，接着介绍正演方法的实现策略，最后通过两个数值模拟对本文提出方法的计算精度和计算效率进行验证.

1 基本原理

1.1 2D黏声介质中基于高斯束的格林函数计算

(1)

其中G表示格林函数；x为散射点；x′为源点；p=(px,pz)为慢度，px和pz分别为慢度在水平方向和垂直方向上的分量；ω为角频率.高斯波束uGB(x,x′,p,ω)可以进一步表示为:

uGB(x,x′,p,ω)=

(2)

其中s、n分别为射线中心坐标系中的坐标；p(s)、q(s)为动力学射线参数，可以通过求解动力学射线方程组获得；τ(s)为沿着射线路径的走时.

令：

(3)

式(1)中格林函数可以表示为：

(4)

其中：

(5)

AGB(x,x′)和τGB(x,x′)分别为复振幅和复走时，下标R和I表示复数的实部和虚部.

黏声介质有多种模型，本文采用的是常Q模型，此时介质黏滞性强弱由与频率无关的品质因子Q决定.地震波在黏声介质中的传播可以看作为以一个复速度在声波介质中传播，当满足条件Q≫1,此时的复速度V(x)可以表示为(Aki and Richards，2002)：

(6)

其中ωr为参考频率.

(7)

其中：

(8)

2D黏声介质中基于高斯束的Born正演公式

Born散射理论可以通过省略高阶项实现针对一次散射波的正演，此时的Born正演公式为(Moser，2012)：

u(xr,xs,ω)=

(9)

其中xs和xr分别表示源点和接收道的位置；F(ω)为震源函数；x为散射点；D为散射点的集合；c0(x)为背景速度；c1(x)为扰动速度；G(x,xs,ω)和G(x,xr,ω)分别为源点和接收点的格林函数，它们的表达式为：

(10)

将式(10)代入式(9)可得黏声介质中基于高斯束的Born正演公式：

u(xr,xs,ω)=

(11)

-iωpLx(xr-L)]，

(12)

其中pLx和pLz为中心道初始慢度在水平和垂直方向上的参数.

结合公式(13)：

(13)

式(11)可以进一步转化为：

×U(L,xs,psx,pLx,ω)exp[-iωpLx(xr-L)

(14)

其中ΔL为窗间隔；ωr为参考频率；w0为初始波束宽度；U(L,xs,psx,pLx,ω)的表达式为：

U(L,xs,psx,pLx,ω)=

(15)

式(15)的含义是将反射率信息映射为局部平面波，接着通过式(14)将局部平面波转化为窗内接收道的地震记录.式(15)为整个计算的核心也决定了正演算法的计算效率，下一小节将着重讨论式(15)的计算方法.

1.2 基于wavelet-bank的时间域实现策略

直接按照式(15)在频率域中对U(L,xs,psx,pLx,ω)进行计算效率并不高，文中将给出高效的时间域计算方法.将式(5)和式(7)代入到式(15)中可得：

U(L,xs,psx,pLx,ω)=

(16)

对U(L,xs,psx,pLx,ω)进行傅里叶反变换可将其转换到时间域：

(17)

(18)

首先针对τI和τQ建立规则的序列：

(19)

(20)

(21)

随后通过双线性插值以及高斯束的振幅和走时将散射点的反射率信息映射到Rjk(t)和Ijk(t)中：

(22)

(23)

整个正演算法的实现流程如图1所示.

2 模型试算

在本节中，将使用一个水平层状黏声模型和一个复杂黏声模型对所提出新正演方法的计算能力进行验证.

2.1 水平层状黏声模型计算

图2为水平层状黏声模型的的速度场分布和反射率信息.模型横向和纵向的网格点数均为301；横向和纵向网格间距分别为20 m和10 m.从上到下各层的声波速度分别为2000 m·s-1、2500 m·s-1、3000 m·s-1；Q值分别为60、50以及40.炮点位于模型顶部中间位置，接收道以炮点为中心向两侧均匀分布201道，道间距为20 m，正演使用的子波为雷克子波，主频为20 Hz，时间采样间隔为4 ms，时间采样点数为750.

图3为水平层状黏声模型采用不同方法获得的正演结果，图3a为声波高斯束Born正演结果，图3b为黏声高斯束Born正演结果.图4进一步展示了两种方法在零偏移距处的波形对比.通过对比可以发现：没有考虑介质黏滞性的声波正演结果中与地层对应的同相轴能量更强，黏声介质高斯束Born正演方法可以很好的将介质黏滞性对地震波场的影响反映到正演结果中.在测试中，由于考虑介质黏滞性的相关计算更加耗时，黏声正演结果总共耗时2.2 s，比声波正演多消耗了约50%的时间.

图1 黏声介质高斯束Born正演实现流程Fig.1 Flow chart of the Born forward modeling for visco-acoustic media using Gaussian beam

2.2 各向异性复杂模型计算

图5为复杂黏声模型的的速度场分布和反射率分布.模型横向网格点数为250；纵向网格点数为750；横向和纵向网格间距分别为10 m和5 m.模型为常Q模型，Q的取值为50.炮点位于横向1000 m处，接收道以炮点为中心向两侧均匀分布50道，道间距为15 m，正演使用的子波为雷克子波，主频为15 Hz，时间采样间隔为4 ms，时间采样点数为750.

图6为复杂黏声模型采用不同方法获得的正演结果，图6a为黏声有限差分方法正演结果，图6b为黏声高斯束Born正演结果.通过对比可以发现：有限差分方法可以模拟直达波等全波场，黏声高斯束Born正演方法只针对一次散射波进行模拟.两种方法得到的一次散射波场结果几乎相同.就计算时间而言，黏声有限差分正演方法共耗时101.4 s，黏声高斯束Born正演方法共耗时1.3 s，可以看出本文提出的正演方法在计算效率上具有很大的优势.

3 结论

本文提出了一种针对二维黏声介质一次散射波场模拟的高斯束Born正演方法.该方法通过叠加不同初射方向的高斯束获取格林函数，克服了传统射线类方法计算格林函数的缺点，保证了格林函数以及波场模拟的精度.同时基于黏声高斯束所包含的信息，新的正演方法采用wavelet-bank方法快速有效的将反射率信息映射成为了窗内的局部平面波，保证了正演算法的计算效率.通过两个数值模型测试表明，本文提出的正演方法具有良好的计算精度，并且计算效率相较于黏声介质有限差分地震波正演要高很多.

图2 水平层状黏声模型(a) 声波速度分布； (b) 反射率分布.Fig.2 Layered visco-acoustic model(a) Acoustic velocity distribution； (b) Reflectivity distribution.

图3 水平层状黏声模型不同方法正演结果(a) 声波高斯束Born正演结果； (b) 黏声高斯束Born正演结果.Fig.3 Modeling results with different modeling methods for layered visco-acoustic model(a) Result of acoustic Gaussian beam Born modeling； (b) Result of visco-acoustic Gaussian beam Born modeling.

图4 零偏移距不同正演方法波形对比虚线为图3a中零偏移距波形，实线为图3b中零偏移距波形.Fig.4 Waveform comparison of zero-offset traces between different modeling methodsDashed line marks the waveform in Fig.3a, solid line is the waveform in Fig.3b.

图5 复杂黏声模型(a) 声波速度分布； (b) 反射率分布.Fig.5 Complex visco-acoustic model(a) Acoustic velocity distribution； (b) Reflectivity distribution.

图6 复杂黏声模型不同方法正演结果(a) 黏声有限差分正演结果； (b) 黏声高斯束Born正演结果.Fig.6 Forward modeling results with different methods for layered visco-acoustic model(a) Result of visco-acoustic finite-difference modeling；(b) Result of visco-acoustic Gaussian beam Born modeling.