地球物理学报  2015, Vol. 58 Issue (4): 1317-1332   PDF    
基于高斯束与高斯波包的Gabor框架散射波模拟方法
李辉, 王华忠    
同济大学海洋与地球科学学院 波现象与反演成像研究组, 上海 200092
摘要:在给出真实模型和相应光滑背景模型的情况下,如何计算扰动模型(散射体)产生的散射波场是一个有实际意义的正演问题.在Gabor变换域描述散射体,且入射波场为短时宽带信号时,散射波场可以在频率域用高斯束或时间域用高斯波包描述.相对于波动方程方法,高斯束和高斯波包的计算效率更高;背景模型光滑时,高斯束和高斯波包方法的精度也接近波动方程方法.文中导出了声波假设下应用高斯束和高斯波包计算散射波的方法.测试分析了高斯波包的计算精度.给出了一般散射体的散射波模拟策略.同时针对一个理论模型完成了本文方法计算散射波的实验,实验结果表明高斯波包散射波计算方法是有效可行的.
关键词散射体     散射波场     Gabor函数     高斯束     高斯波包    
A scattering wave modeling method using Gaussian beam and Gaussian packet in the Gabor domain
LI Hui, WANG Hua-Zhong    
Wave Phenomena and Inversion Imaging Group (WPI), Tongji University, Shanghai 200092, China
Abstract: How to calculate the wavefield that scattered by a perturbed model (scatterers) in a smooth background model is a practical problem. If scatterers are described in the Gabor domain, and the incident waves are short-duration and broad-band, then the scatter wavefield can be calculated using Gaussian beam in the frequency domain or using Gaussian packet in the time domain. Compared with wave-equation, Gaussian beam/packet shares higher efficiency. Moreover, the accuracies of Gaussian beam/packet are comparable with wave-equation in the case of smooth media. The method that calculates scatter waves using Gaussian beam and Gaussian packet in an acoustic medium is proposed. The accuracy of the modeling method is also analyzed.The incident waves and Green's function can be both written as the expressions with amplitude and travel-time under WKBJ approximation. Substituting these expressions into the formula of Born approximation of scatter waves in the acoustic medium, then the analytic expression of scatter waves that are scattered by Gabor-scatterer is deduced.The initial Hessian of scatter Gaussian beam's travel-time with respect to space coordinates is the same as that in Klimeš's work (2012). We present the initial values of Gaussian beam's dominate frequency, ray parameters, and amplitude in 3-dimensions, respectively. The initial amplitude of scatter Gaussian beam in 2-dimensions is also presented.Similar to scatter Gaussian beam, the initial Hessian of scatter Gaussian packet's travel-time is the same as that of Klimeš's work (2012), too. Scatter Gaussian packet's dominant frequency and initial ray parameters can be calculated with the same formulas as the scatter Gaussian beam. Both formulas for calculating initial amplitudes of scatter Gaussian packets in 3-dimensions and 2-dimensions are deduced.Finally, the strategy of simulating scatter wave with Gaussian beam and Gaussian packet can be concluded as a flow chart in the main body of the paper. The wavefield scattered by a Gabor-scatterer is simulated with the finite-difference (FD for short) method and Gaussian packet. The two are almost the same in the time domain and very similar with each other in the frequency domain. The accuracy of scatter Gaussian packet is higher as the scatter angle becomes smaller, and vice versa. The scattered waves of a synthetic model are examined, and the common shot gathers calculated using FD method and Gaussian packet are comparable, except boundaries of shot gather. For comparing details of wavefields simulated using the above two methods, two groups of traces are extracted from common-shot-gathers. Comparison of waveforms show that amplitudes of different methods are almost the same while phases of wavelets are distinguishable. Accuracies of waveforms calculated with Gaussian packet are lower when offset becomes bigger.Calculation of the wavefield scattered by Gabor-scatterer using Gaussian beam and Gaussian packet are deduced for constant-density acoustic medium. Compared with scatter wave modeling in an anisotropic elastic medium, the number and initial amplitudes of scatter Gaussian beams/packets are both distinguishable here. Scatter Gaussian beam/packet's initial amplitudes are disparate in 2- and 3-dimensional media. With the proposed scatter theory, wavefields scattered by Gabor-scatterer calculated using the Gaussian packet and FD method are compared. They are comparable with each other in the time domain while there is a little difference in the frequency domain, and the reason has been analyzed in this article. Accuracy of simulation using Gaussian packet is influenced by scatter angles. The accuracy is higher when scatter angle becomes smaller and vice versa. Common-shot-gathers observed on the earth surface of a general model are simulated with the proposed Gaussian packet approach. Comparison of the presented scheme and FD method shows that our method is effective, although waveforms of common-shot-gathers simulated using Gaussian packet are slightly different from that simulated using the FD method.
Key words: Scatterer     Scattering wave     Gabor function     Gaussian beam     Gaussian packet    
1 引言

已知背景模型与其中传播的背景波场,同时给出散射体(模型扰动部分),如何计算散射体产生的散射波场是一个有应用价值的工作.散射波的正、反演与偏移成像有密切的联系.Cheng和Coen(1984)详细推导并分析了常速介质中散射波反演与F-K偏移(Stolt,1978)的关系,Bleistein等(1984)指出WKBJ近似下散射波反演的一次迭代与Kirchhoff偏移(Schneider,1978)是等价的.最小二乘偏移即从散射波反演中发展而来,其中的正演过程即为散射波模拟(Beydoun and Mendes, 1989任浩然,2011).反射全波形反演(RFWI)的研究方兴未艾(Xu et al., 2012Ma and Hale, 2013),其中关键步骤是利用反偏移计算目标函数的梯度,反偏移在背景介质中利用模型高波数成分(即散射体)模拟散射波.事实上,Mora(1989)在1980年代末就已经分析了反演可分为“偏移”和“层析”两部分,其中层析部分的模拟过程为散射波模拟.从上述分析可知,在勘探地震领域,数据处理的核心部分——偏移与反演均与散射波的正演有着密不可分的关系.所以散射波正演理论研究对勘探地震的发展从基础上存在显著的推动作用.鉴于散射波模拟有重要的理论和应用价值,本文致力于在Gabor域描述散射体时散射波模拟方法的研究.

在网格中描述散射体的空间分布,每一个散射体单元是一个网格点,在此基础上,传统散射波应用(最小二乘偏移等)中需要每个网格点的散射波,实际计算时应用格林函数计算每个检波点观测到的散射波(任浩然,2011).计算格林函数的方法有射线法(Bleistein and Cohen, 1979Bleistein et al.,1984Červený,1992)、波动方程法(任浩然,2011Ren et al., 2013)和高斯束/高斯波包法(Červený and Psencik,1983李辉等,2012)等.射线方法效率较高,但受到射线理论的约束,无法处理正演时的焦散及阴影区问题(;Červený,1992),且难以计算波场中的子波相位信息;波动方程的计算效率特别低;高斯束/高斯波包计算格林函数时需要对不同初始角度的高斯束/高斯波包进行积分,体现不出高斯束/高斯波包的效率优势.总之,在传统网格描述散射体的基础上,通过格林函数计算散射波的方法存在效率低等诸多问题,同时由于网格描述散射体时冗余性较高,进一步降低了散射波模拟的效率.

在Gabor域描述散射体时,每一个散射体单元变成一个Gabor函数,一个散射体单元散射出的散射波不再是格林函数.在一阶Born近似和零阶WKBJ近似假设下,Klimeš(2012)推导出Gabor函数散射体产生的散射波场可以用高斯束(频率域)或高斯波包(时间域)计算,在各向异性弹性介质中(模型参数为刚度系数和密度),一个Gabor函数表达的散射体单元所产生的散射波可利用最多5个高斯波包(时间域)或高斯束(频率域)描述.高斯束( ;Červený,et al, 1982Popov,1982)和高斯波包(Ralston,1983Klimeš 2004)是介于射线理论和波动理论之间的一类方法,效率显著高于波动方程,同时可避免射线法的焦散和阴影区等问题.

图 1示意了用网格描述散射体,背景波场入射(图 1b中的红色弧线为波前面)时一个网格表达的散射体单元(图 1中的蓝色点)所产生的散射波如图 1c中红色圆圈所示,可以用格林函数计算.类似地,散射体用Gabor函数作为框架描述时,一个Gabor框架函数表达的散射体单元(如图 2a所示)所产生的散射波如图 2c所示,是一个高斯波包.

图 1 网格描述的点散射体(a),入射背景波到达散射体(b),背景波场经过散射体之后产生散射波(c)
图中蓝色点是散射体,红色圆弧是背景波场的波前面,红色圆圈是散射波场.
Fig. 1 The one-grid-perturbation(a),the background wave-field(b), and the scattered Green′s function(c)
The blue dot,the red arc and the red circle are the perturbed model,the wave-front of the background wave, and the scattered Green′s function,respectively.

图 2 Gabor函数描述的散射体(a),入射背景波到达散射体(b),背景波场经过散射体之后产生散射波(c)
图中蓝色震荡图案表示散射体,红色圆弧是背景波场的波前面,红色震荡图案是散射波场(引自Klimeš, 2012).
Fig. 2 The Gabor-perturbation(a),the background wave-field(b),the scattered Gaussian packet(c)
The blue pattern,the red arc, and the red packet are the perturbed model,the wave-front of the background wave, and the scattered wave packet,respectively(Klimeš, 2012).

Klimeš(2012)工作中的模型参数为刚度系数与密度,相应的散射波包括纵波和两类横波等多种波现象,虽然该理论具有普适性,但也限制了此方法的实际应用,不利于偏移和反演,所以其文中仅做了理论分析工作.本文在常密度声波假设下,以声波速度为模型参数,参照Klimeš的思路推导了声波假设下的散射波正演理论,并推导了在频率域利用高斯束与在时间域利用高斯波包计算散射波场的具体方法.与Klimeš的工作不同,声波近似下一个Gabor函数表达的散射体最多产生一个散射波包.

既然单个Gabor函数描述的散射体所产生的散射波可用高斯束或高斯波包计算,那么在一阶Born近似下,任意空间分布的散射体所产生的散射波均可用高斯束或高斯波包计算.

本文第2节介绍了零阶WKBJ近似下Gabor散射体产生的散射波利用振幅和走时表达的形式,并在常密度声波介质中分析了此表达形式下散射波的特点.第3节和第4节在三维声波介质中分别推导了零阶WKBJ近似下散射波利用高斯束和高斯波包的具体计算公式,尤其是高斯束和高斯波包的初始振幅.第5节给出了二维声波介质中高斯束和高斯波包初始振幅的计算方式.数值实验部分首先利用本文方法和有限差分计算了Gabor函数表达的散射体所产生的散射波场,并对比二者分析了高斯波包计算的散射波精度.多Gabor函数散射波模拟的实验给出了地表数据的模拟策略.同时针对一个理论模型,应用高斯波包计算散射波的方法模拟出了地表观测波场,该计算结果与有限差分法模拟的结果对比,可发现两者几乎完全一致,只是观测边界附近有差异,说明了本文散射波模拟方法是有效可行的.

本文不复述高斯束和高斯波包的理论,对文中理论部分的理解需要一定程度关于上述两者的理论基础.

2 零阶WKBJ近似下Gabor散射体产生的散射波场

高斯束和高斯波包在描述波场传播时引入了零阶WKBJ近似.本节首先介绍零阶WKBJ近似下Gabor散射体产生的散射波场的具体表达形式,然后分析散射波场在常密度声波介质中的特点.为高斯束和高斯波包计算散射波的方法提供理论基础.

常密度声波方程中的一阶Born近似下散射波us的计算写成如下形式:

其中u0是在背景速度中传播的背景波场,G是背景速度中的格林函数,σ是模型扰动,或称散射体.散射体和背景速度c0及扰动后速度c的关系为

令散射体σ的空间分布为Gabor函数,满足如下公式:

其中σ0αα及对称矩阵α分别是Gabor函数的振幅、中心点、波数和高斯衰减系数.图 3示意了二维空间中一个Gabor函数的形态.从(3)式及图 3可知Gabor函数是加高斯衰减窗的简谐震荡函数.零阶WKBJ近似下,入射波场写成如下用振幅和走时函数近似的表达形式:

其中A和T是背景波场传播到′点的振幅和走时.三维介质中格林函数的零阶WKBJ形式为(Bleistein,2008)

其中L表示几何扩散,τ(′)是从点′到的走时.
图 3 二维空间中的Gabor函数示意图Fig. 3 A Gabor function in 2-dimensions

将散射体的Gabor函数表达式(3),背景波场与格林函数的零阶WKBJ近似(4)和(5)代入一阶Born近似下的散射波场积分公式(1)中,得到Gabor散射体的散射波积分形式(Klimeš, 2012)

其中

式(8)的表述中使用了爱因斯坦求和准则.Klimeš(2012)指出上述散射波的能量集中在主散射射线路径附近,主散射射线上的点为0Δ0到主散射射线附近点的矢量.公式(8)中变量在Klimeš(2012)的文章中有详细的解释,为了保持本文理论的完整性,我们把上述变量的物理意义表述在附录A中.

0位于主散射射线上,当Δ=时,散射波us() 的振幅最大,即式(8)中振幅衰减最小,此时有

从附录A中T,iα和τ,i0′的定义可知,T,iα是背景波场入射到Gabor函数中心点时射线参数矢量的第i个分量,记作Pi;(-τ,i0′)是主散射射线在Gabor函数中心点处射线参数矢量的第i个分量,记作pi0.从而(9)式也可写成

(10)式给出了散射波场主散射射线的初始慢度向量的计算公式.(10)式中频率ω未知,接下来推导其计算方式.由常密度声波介质中慢度向量的定义可知,慢度向量的模等于速度的倒数,Pi和pi0满足如下关系式:

cα是Gabor函数中心点α处的背景介质速度.结合(10)式和(11)式,得到Gabor函数中参数α与背景波场慢度向量 P、散射波主频ω之间的关系式

当入射波场方向与Gabor函数的矢量参数α垂直时,(12)式不存在关于1/ω的非平凡解,此时不产生散射波,否则(12)式存在唯一关于1/ω的非平凡解,此时散射波的主频为

公式(13)确定的频率成分对应的散射波场能量最大,这里称此频率为散射波的主频.唯一的散射波主频代入主散射射线初始射线参数的计算公式(10)中可得唯一的主散射射线初始射线参数,即主散射射线唯一,这点与各向异性弹性介质中不同(Klimeš, 2012).

至此推导出在零阶WKBJ近似下通过射线理论计算Gabor散射体产生散射波场的具体方式:首先利用公式(13)计算出散射波主频;再通过(10)式计算得到主散射射线的初始射线参数;射线追踪得到主散射射线路径;最后利用追踪出的射线路径,结合(6)、(7)和(8)计算出射线附近任意点任意频率的散射波.下文中我们将在此基础上指出上述散射波场的计算可利用高斯束和高斯波包实现.事实上,(13)式确定的主频是Gabor散射体产生的散射高斯波包的主频,(10)式确定的射线参数是高斯束和高斯波包计算散射波时中心射线的初始射线参数.

3 频率域的散射波计算方法——高斯束

本节讨论利用高斯束在频率域计算Gabor散射体所产生散射波的具体方法及实现策略.高斯束计算散射波场时可写成如下形式:

其中(s,q)是射线中心坐标系中的坐标,a(s)和T(s,q)是振幅项和复数相位项,其中复数相位项包含振幅在q方向的高斯衰减.

在各向异性弹性介质中(模型参数为刚度系数和密度),Klimeš(2012)提出了Gabor散射体产生的散射波场在频率域利用高斯束计算的方法.声波近似下高斯束计算散射波复数相位部分的公式与Klimeš的工作相同,为了保持本文理论过程的完整性,附录B中在第2节理论的基础上直接给出高斯束计算复数相位部分的具体公式,详细的推导过程参考Klimeš的工作.利用高斯束计算散射波的振幅部分,Klimeš(2012)的工作和本文方法不同,本节将对振幅部分的高斯束计算方法做详细推导.

由(6)、(7)式和(8)式可知,散射波场的振幅项可以写成

格林函数的几何扩散项L为(Bleistein,2008)

其中Q2是附录B中(B4)式所示动力学射线追踪系统的传播矩阵元素.公式(16)表达的几何扩散代入振幅项(15),经过简单推导,得到散射波振幅如下表达式:

其中Q是动力学射线参数,计算方式参考附录B中(B3)式,m的物理意义及计算公式参考附录B中(B12)式.另外,高斯束理论中的振幅计算公式有如下形式(Popov,1982):

对比Gabor散射体产生的散射波振幅计算公式(17)与高斯束的振幅计算公式(18),可得利用高斯 束计算Gabor散射体产生的散射波时的初始振幅为

至此,可总结利用高斯束计算Gabor散射体所产生的散射波的计算过程:

I:利用公式(13)计算散射波主频;

II:利用公式(10)计算散射波的主散射射线的初始射线参数;

III:参考附录B中公式(B5),计算动力学射线参数的初始值;

IV:参考附录B中公式(B15)和(B16),计算复数相位函数关于频率的二阶导数的初始值;

V:参考公式(19)计算高斯束初始振幅;

VI:运动学射线追踪得到中心射线路径;

VII:动力学射线追踪得到动力学射线参数;

VIII:利用公式(B13)和(B14)计算复数相位函数关于频率的二阶导数;

IX:利用公式(B1)和(B2)计算高斯束的复数相位部分;

X:利用公式(18)计算高斯束的振幅部分;

XI:利用高斯束振幅和复数相位部分,通过(14)式计算散射波场.

4 时间域的散射波计算方法——高斯波包

第3节中讨论了Gabor散射体产生的散射波场在频率域用高斯束计算的过程.对所有频率的高斯束按频率积分即可得到Gabor散射体产生的散射波场在时间域的表达形式,得到时间域计算散射波场的方式——高斯波包.高斯波包计算散射波场写成振幅和复数相位的表达形式:

其中aGP(,t)和T(,t)是振幅项和复数相位项,复数相位项包含振幅在空间中的高斯衰减.各向异性弹性介质(模型参数为刚度系数和密度)中,散射高斯波包复数相位部分的公式与常密度声波介质中相同,振幅部分不同.同样为了本文理论过程的完整性,附录C中直接给出Klimeš(2012)用高斯波包计算复数相位的公式.本节中我们将给出振幅部分的计算方法及相应的推导过程.已知初始振幅,高斯波包振幅在时空域的计算公式为

其中aGP(α,t)是高斯波包初始振幅,Q1和Q 2是公式(B4)所示动力学射线追踪传播矩阵的元素,MGPα是高斯波包相位函数在射线中心坐标系中q方向的二阶导数初值,详见(C12)式.

高斯波包(20)式中振幅部分(21)式的计算可通过高斯束的频率积分得到,积分结果的振幅在散射点的表达式即为高斯波包的初始振幅,具体形式为

其中aGB(α,ω0)是频率为ω0的高斯束初始振幅,M44α是高斯束相位函数在频率方向的二阶导数,参见附录B中的(B16)式.

结合高斯束初始振幅(19)、(B16)和(22)式,高斯波包的初始振幅可表达如下:

(23)式中k3(q)是Gabor函数的参数波数向量α在射线中心坐标系中沿射线方向的分量.用汉密尔顿常量H(Červený,2001)表示k3(q),可把(23)式所示的高斯波包初始振幅转换成如下形式:

本文中与(24)式对应的汉密尔顿常量为

把(25)式代入(24)式得到散射高斯波包振幅初值的最终表达式

至此Gabor散射体的散射波利用高斯波包的计算可由(20)式实现,高斯波包的初值公式通过(C15)—(C17)及(26)获取.具体计算过程可总结为如下若干步骤:

I:利用公式(13)计算散射波主频;

II:利用公式(10)计算散射波的主散射射线的初始射线参数;

III:参考附录B中公式(B5),计算动力学射线参数的初始值;

IV:参考附录B中公式(B15)和(B16),计算高斯束复数相位函数关于频率的二阶导数的初始值;

V:参考附录C中公式(C15)、(C16)和(C17),计算高斯波包复数相位函数在时空域二阶导数的初始值;

VI:参考公式(26)计算高斯波包初始振幅;

VII:运动学射线追踪得到中心射线路径;

VIII:动力学射线追踪得到动力学射线参数;

IX:参考附录C中公式(C5)—(C9)计算复数相位函数在时空域的二阶导数;

X:利用公式(C1)计算高斯波包的复数相位部分;

XI:利用公式(21)式计算高斯波包的振幅部分;

XII:利用高斯波包振幅和复数相位部分,通过(20)式计算散射波场.

5 二维介质中的散射波振幅

公式(19)和(26)分别给出了三维介质中高斯束和高斯波包的初始振幅,本节中将指出二维介质中振幅计算与三维的不同之处,并给出具体的初始振幅计算公式.

二维介质中格林函数的零阶WKBJ近似形式为(Bleistein,2008)

其中Q2是动力学射线追踪传播矩阵中的元素.通过与第3节中类似的推导过程,可得散射高斯束的初始振幅

同样参考第4节中的理论推导过程,得到散射高斯波包的初始振幅为

6 数值实验

本节在二维介质中对上述理论进行数值实验.有限差分法(FD)求解声波方程是一种比较成熟且应用广泛的数值模拟方法,所以我们以FD模拟的散射波为标准分析本文方法的有效性和模拟结果的精度.实验分成三部分.首先针对一个Gabor函数表达的散射体,在时间域利用高斯波包进行散射波正演,并把正演结果与FD模拟的散射波对比,以测试高斯波包计算散射波的有效性与精度.之后针对多Gabor函数表达的散射体,给出高斯波包计算散射波的具体策略,并将此策略下模拟的散射波与FD的计算结果进行对比,以证明高斯波包叠加计算散射波的有效性,同时指出一个完整的反射层对观测数据有贡献的散射体部分.最后针对一个理论模型通过高斯波包模拟出地表观测数据,并将模拟结果与FD对比.

6.1 散射波正演方法精度测试

本小节中对高斯波包正演过程的精度做测试,对比高斯波包计算的散射波和FD的模拟结果,并分析二者出现差异的原因.影响散射波模拟的因素较复杂,包括Gabor散射体参数、背景波场、观测方式、正演算子等.为测试本文方法计算散射波包振幅的精度,我们对相同的Gabor散射体、背景波场和观测方式,分别利用高斯波包和FD计算观测点接收到的散射波,在时间域和频率域对比观测到的波形以及散射波的振幅谱,以FD计算的结果为标准,分析高斯波包的计算精度.

我们设计了一个常速背景模型,模型参数如表 1所示.扰动模型在Gabor变换域描述,这里针对一个Gabor框架散射出的波包进行测试,Gabor函数如公式(3)所示,用于测试的Gabor框架参数见表 2.扰动后的模型,即扰动模型对应的速度与背景模型速度场之和如图 4所示,速度扰动以xα点为中心,震荡方向与z轴平行,速度最大扰动量约216 m·s-1. 背景波场是在背景介质中传播的点源波场,震源位置为(2240 m,24 m),高斯波包方法的背景波场利用(27)式所示零阶WKBJ近似下的格林函数计算,检 波点分布于深度为24 m的水平面上.FD和高斯波包计算的散射波如图 5所示,图 5a中FD计算的散射波与图 5b中高斯波包计算的散射波非常相似.为了进一步对比两种方法计算的散射波场,从接收到的散射波记录中抽出一道,对二者进行波形对比,该道的检波点位于x=2760 m处,是散射波包中心点的x坐标.图 6a在时间域对比了两种方法计算的散射波单道信号,从中可看出高斯波包和FD计算的散射波没有明显差异.图 6b是两种方法计算的散射波单道信号振幅谱对比,高斯波包计算的散射波是纯粹的窄带信号;而FD的散射波除相近频率范围的窄带信号之外,还有两个相对高频的信号成分,FD散射波信号包含多个窄带信号.此外FD散射信号的第一个窄带频谱振幅也小于高斯波包散射信号.这些差异将在下文中继续讨论.

表 1 背景模型参数Table 1 The parameters of the background model

表 2 Gabor扰动模型参数Table 2 The parameters of the perturbed Gabor function

图 4 扰动后的速度模型Fig. 4 The perturbed velocity model

图 5 检波点观测到的散射波.(a)FD计算的散射波;(b)高斯波包计算的散射波Fig. 5 The scattering waves simulated with FD method(a) and Gaussian packet(b)

图 6 不同方法计算散射波的单道信号对比.(a)时间域对比;(b)振幅谱对比
红色曲线是FD计算结果,蓝色曲线是高斯波包计算结果.
Fig. 6 The single trace signal of the scattering wave in time domain(a) and frequency domain(b)
The red curve is calculated with FD method and the blue curve is calculated with Gaussian packet.

接下来对比背景波场入射方向不同时两种方法计算的振幅,以考察背景波场对振幅计算精度的影响.从表 2图 4中可知,Gabor扰动模型的震荡方向与z轴平行.类似于反射波入射角定义方式,这里定义散射的入射角度——在Gabor中心点处入射的背景波场方向与Gabor函数震荡方向的夹角.入射角不同时,散射波包的振幅、频率均不同,此时高斯波包计算散射波场的精度也不同.令一组震源位于以(xα,zα)为圆心、半径为2490 m的圆上,如图 7所示,红色五角星标识震源的分布,震源和圆心连线与z轴的夹角范围为0°~60°,该范围内等间隔地分布11个震源,即相邻震源间的弧度为6×2π.以这11个震源在背景模型中产生的波场为背景波场,分别计算图 4中扰动模型的散射场,每个震源对应的检波点分布于与震源深度相同的水平面上.由于背景介质是常速模型,通过高斯波包计算散射波的理论可知观测面上散射波包中心点位于与震源轴对称的位置处,如图 7中蓝色三角形所示,对称轴是x=2500 m处的竖直直线.对散射波包中心点处接收到的散射波信号求包络,包络的极值被看作波包的振幅,11次实验测量出的散射波振幅分布及两种算法的相对误差展示在图 8中.无论FD还是高斯波包计算出的散射波振幅均出现类似反射波的AVA规律,随着入射角的增大,散射波的振幅也增大.入射角较小时,散射波振幅较小,高斯波包计算的散射波振幅误差也比较小;入射角增大,散射振幅增大,高斯波包计算的散射波振幅误差也随之增大.总体来说,除个别角度,高斯波包计算的散射波振幅相对误差小于9%,在可接受的范围之内.

图 7 振幅分析的炮检点分布示意图Fig. 7 The sketch map of the locations of sources and receivers

图 8 不同散射角的散射波振幅对比
(a)散射波绝对振幅对比,FD=有限差分,GP=高斯波包;(b)高斯波包散射波相对FD散射波振幅的相对误差.
Fig. 8 The scattering wave′s amplitude versus scattering angle
(a)The comparing of the absolute amplitudes,FD=Finite difference,GP=Gaussian packet;(b)The relative amplitude errors calculated with Gaussian packet.

上述高斯波包计算散射波的正演实验中,无论从时间域波形、频谱还是振幅的角度看,该散射波正演方法或多或少均存在一定误差.这是因为高斯波包正演理论中的几个假设:零阶WKBJ近似、一阶Born近似和入射波无限宽频.尤其是对入射波场的无限宽频假设,由于FD方法的入射波场不满足该假设,入射波形是有时间延续长度的子波,而非脉冲函数,致使散射波在频率域出现多个窄带信号(图 6b),同时也影响到FD计算的散射波包振幅与高斯波包振幅的差异.

6.2 多Gabor散射体模拟策略

这一部分针对多Gabor函数描述的散射体,利用高斯波包计算其散射出的散射波.既然单个Gabor函数描述的散射体所产生的散射波可用高斯束或高斯波包计算,那么在一阶Born近似下,任意空间分布的散射体所产生的散射波均可用高斯束或高斯波包计算.首先把散射体做Gabor分解,得到Gabor框架的线性组合,对每一个Gabor框架函数均利用高斯束或高斯波包计算出相应的散射波场,最终散射模型的散射波场即为Gabor框架函数产生散射高斯束或高斯波包的线性组合,其系数和Gabor框架的系数相同.具体策略是首先把散射体分解成Gabor函数线性叠加的形式

其中gi()和ci分别是第i个Gabor函数及其对应的分解系数.然后利用第4节中给出的算法计算每一个Gabor函数gi()产生的散射波,根据一阶Born近似把所有Gabor函数产生的散射波线性叠加即可得到整个散射体的散射波

其中uig()是gi()产生的散射波,可通过高斯波包方法计算.(31)式中散射波线性叠加的系数和(30)式中散射体线性叠加的系数相同.

按照上述散射波计算策略,这里测试了三个散射体,分别把本文方法与FD的模拟结果对比.背景介质模型与6.1节实验中相同,参数如表 1所示.震源和观测点均位于地表,震源的x坐标为1000 m,观测点位于整个地表 501个CDP点.

第一个散射体由两个Gabor函数组成,两个Gabor函数的空间中心点不同,分别位于(2500 m,2000 m)和(2500 m,2500m),其他参数均相同,相应的扰动后速度场如图 9a所示.从速度模型中看出,两个散射体的垂向位置不同,且在空间上完全分开.高斯波包和FD计算的散射波展示在图 9b和9c中,二者几乎完全一致.

图 9 两个Gabor函数表达的散射体对应的扰动后速度模型(a),高斯波包(b)和有限差分(c)模拟的散射波场Fig. 9 The perturbed velocity expressed with two Gabor functions(a),the scattering wave simulated with(b)Gaussian packet and (c)FD method

第二个散射体也由两个Gabor函数组成,中心点空间位置分别为(2000 m,2500 m)和(2500 m,2500 m),其他参数两者均相同.图 10a是扰动后的 速度场,从中得知两个Gabor函数相互重叠,空间上没有分开.高斯波包和FD模拟的散射波分别如图 10b图 10c所示,可见空间上没有分开的两个Gabor函数描述的散射体产生的 散射波利用高斯波包计算的结果与FD差别也很小.

图 10 两个Gabor函数表达的散射体对应的扰动后速度模型(a),高斯波包(b)和有限差分(c)模拟的散射波场Fig. 10 The perturbed velocity expressed with two Gabor functions(a),the scattering wave simulated with(b)Gaussian packet and (c)FD method

第三个散射体是深度为2500 m的水平反射层,扰动后速度模型如图 11a所示.利用特征Gabor分解方法(李辉等,2014)把此散射体分解成21个Gabor函数,通过高斯波包法计算出每个Gabor函数的散射波,按公式(31)叠加后的整个散射波场如图 11b所示,图 11c是利用FD计算的散射波场,对比两者可知在单层模型中本文方法仍然有效.用Gabor函数描述散射体,通过高斯波包正演算子计算散射波的一个特点是可以考察散射体哪些部分对观测到的散射波有贡献.因为整体散射波的计算由所有Gabor散射体的高斯波包叠加得到,所以对观测波场有贡献的高斯波包对应的Gabor散射体被认为是有效的.对图 11b中整体散射波有贡献的Gabor函数线性叠加结果如图 11d所示,由此可知,利用图 11b所示的单炮数据做偏移或反演时只能得到图 11d中的一部分结果(剖面、成像道集、扰动模型等),因此本文模拟方法在实际应用时也不需要对所有Gabor函数表达的模型扰动进行偏移或反演.

图 11 水平散射体扰动后的速度模型(a),高斯波包(b)和有限差分(c)模拟的散射波场,对观测波场有贡献的散射体成分(d)Fig. 11 The perturbed velocity with a horizontal reflectivity(a),the scattering wave simulated with(b)Gaussian packet and (c)FD method, and the segmental reflectivity that has contribution to received scattering wave
6.3 理论模型实验

针对一个理论模型,利用高斯波包模拟散射波的理论,结合公式(30)和(31)的散射体分解及散射波合成策略,模拟地表观测的单炮数据,并把模拟结果与FD对比.理论模型的速度参数如图 12a所示,模型参数和观测参数列于表 3中.背景模型在理论模型的基础上通过平滑得到,背景模型速度场如图 12b所示,标准模型与背景模型作差可得到速度差,利用公式(2)计算得到散射体的空间分布,准确速度差与散射体分别展示于图 12c图 12d中,散射体分布在反射界面附近.利用第2节和第4节的公式及策略,计算图 12d表示的散射体产生的散射波.正演实验中分别利用高斯波包和FD模拟了表 3所示观测系统中201炮数据的散射波场,图 13a、13b、13c是利用FD模拟的三个共炮点道集,炮点分别位于(1500 m,10 m)、(3000 m,10 m)、(4500 m,10 m),图 14a、14b、14c是高斯波包法计算的相应的三炮散射波,两种方法计算结果的走时、振幅和子波相位特征均一致,只是在观测范围的边界附近子波有差异,如图中箭头标识区域所示,相对FD模拟结果,高斯波包的子波较窄,主频较高.以FD模拟结果为标准,高斯波包在计算大角度散射波时精度比小角度低,这与6.1节的结论一致.

图 12 理论模型速度场(a),光滑的背景速度场(b),相应速度差(c)与散射体(d)Fig. 12 The true velocity model(a),the smoothed one(b),the difference of the above two velocity models(c), and the scatterer(d)

表 3 理论模型参数与观测系统参数Table 3 The parameters of model and geometry

图 13 有限差分法模拟的共炮点道集
三炮数据的震源分别位于(a)(1500 m,10 m),(b)(3000 m,10 m),(c)(4500 m,10 m).
Fig. 13 The common-shot-gathers simulated with FD method
The sources located at(a)(1500 m,10 m),(b)(3000 m,10 m),(c)(4500 m,10 m).

图 14 高斯波包模拟的共炮点道集
三炮数据的震源分别位于(a)(1500 m,10 m),(b)(3000 m,10 m),(c)(4500 m,10 m).
Fig. 14 The common-shot-gathers simulated with Gaussian packet
The sources located at(a)(1500 m,10 m),(b)(3000 m,10 m),(c)(4500 m,10 m).
7 结论

速度模型可分成背景部分和扰动部分(散射体),波场在背景模型中传播遇到散射体时产生散射波场.当散射体在空间的分布是一个Gabor函数时,相应的散射波可用高斯束(频率域)或高斯波包(时间域)描述.这就给散射波的模拟与偏移、层析等应用提供了新的思路——在Gabor变换域以Gabor函数为框架描述散射体,散射波的正演工具为高斯束和高斯波包.

本文在常密度声波假设下推导了利用高斯束和高斯波包计算Gabor散射体产生散射波的具体方法.与刚度系数和密度为模型参数的弹性理论下的散射波模拟方法相比,常密度声波方程中散射波模拟方法的主要不同之处在于散射高斯束(或高斯波包)的个数以及高斯束(或高斯波包)的初始振幅.此外,三维介质和二维介质中的高斯束(或高斯波包)初始振幅也有差异.常密度声波假设下高斯束与高斯波包初始振幅在三维和二维介质中的计算方法是本文的主要理论工作.

在上述理论方法的基础上,本文在时间域利用高斯波包算子实现了一个散射体单元(Gabor函数)产生的散射波的计算.计算结果与FD对比,两者在时间域非常相似,频率域中高斯波包法与FD计算的结果稍有差异,前者是一个单一的窄带信号,后者是若干个窄带信号,但后者的主要能量仍然集中于一个窄带信号中,文中对此现象做了具体分析.

同一高斯波包描述的散射体,背景波场入射方向不同时,高斯波包计算散射波的精度不同,文中实验对此做了精度分析,入射角度越小精度越高.总体误差低于9%,可见本文推导的振幅计算公式精度较高.

多Gabor函数表达的散射体所产生的散射波实验中给出了计算散射波的具体策略.从对水平层散射体的高斯波包散射波模拟实验中,通过分析可知哪些局部散射体对观测数据有影响,哪些局部散射体对观测数据没有影响,这点在偏移和层析等应用中很重要.

针对理论速度模型,模拟了相应的地表观测数据.本文方法与FD的对比说明本文方法可以保证模拟波场的精度.但在大偏移距部分稍有差异,这与前文中 “大角度模拟精度小于小角度”的结论一致.

高斯波包计算散射波的方法主要基于如下几个假设:零阶WKBJ近似、一阶Born近似和入射波无限宽频.这些假设限制了方法的精度.但由于高斯波包方法效率高、应用灵活、易于处理局部数据,所以本文散射波模拟方法在偏移、反演中有一定的应用前景.

致谢 感谢中国石化集团公司科技部及中国石化石油物探技术研究院(南京)对WPI研究组及本项研究的支持.

附录A 公式(8)中变量的物理意义

公式(4)所示背景波场的走时T(′)在Gabor散射体扰动的中心点α的值记为T(α). T(α)在点α附近的一阶和二阶微分简写为

公式(5)所示格林函数的走时为τ(′),′和分别是格林函数的源点和观测点.τ(′)在点′附近的一阶和二阶微分简写成

散射波场在主散射路径上的点为0,点源为α的格林函数走时及其微分在点0处的值τ(0α)、 τ,i(0α)、τ,ij(0α)、τ′ ,i(0α)、τ′,ij(0α)和τ″ ,ij(0α) 分别简写成τ0、τ,i0、τ,ij0、τ,i0′、τ,ij0′和τ0″,ij.点源为α的格林函数走时及其微分在点处的值τ(α)、τ′ ,i(α)和τ′ ,ij(α)分别简记为τα()、τ,iα′()和τ,ijα″().

背景波场走时和格林函数走时之和是单点散射波(散射体是一个点)的走时,写成如下形式:

把θ(′)在Gabor函数的中心点α附近做泰勒展开,有

其中

θα()、θα′,i()和θα″,ij()在主散射路径上的点0附近做泰勒展开,有

附录B 高斯束计算Gabor散射体产生的散射波相位部分

这里首先给出复数相位的计算公式,之后讨论公式中相位函数二阶导数的计算方法.(8)式中主频ω0对应的相位部分利用高斯束计算的公式为

此时0取点在射线上的垂足,MIJ是相位函数在射线中心坐标系q1和q 2方向的Hessian矩阵M中的元素.主频ω0附近的频率ω对应的散射波相位部分计算公式写成

其中MI4是相位函数在空间和频率方向混合二阶导数向量M4的元素,M44则表达了相位函数在频率方向的二阶导数.

已知初始矩阵Mα时,高斯束相位函数在射线中心坐标系q1和q 2方向的Hessian矩阵M

其中

是动力学射线追踪传播矩阵(Červený,2001).高斯束相位函数在垂直射线方向的空间二阶导数初值为

其中,T(q)是入射波走时在射线中心坐标系中点α处对空间微分Hessian矩阵的左上2×2子矩阵,即相位函数在射线中心坐标系q1和q 2方向的二阶导数,射线中心坐标系中完整的Hessian矩阵写成

矩阵(q)可利用笛卡儿坐标系下Hessian矩阵通过Jacobian矩阵h计算得到,其中第a行第b列元素

Jacobian矩阵是笛卡儿坐标系向射线中心坐标系变换的旋转矩阵,

其中第i列向量是笛卡儿坐标系中的第i个坐标轴单位向量在射线中心坐标系中的元素,第a行向量是射线中心坐标系中的第a个坐标轴单位向量在笛卡儿坐标系中的元素. (q)是Gabor函数参数 Kα在射线中心坐标系中的表达形式,

同样可利用坐标转化矩阵计算得到(q)

(B5)式中的 m和m分别为

相位函数的空间一阶导数在频率方向的一阶导数M4和频率方向的二阶导数M44可利用如下公式计算:

其初值为

初值公式中 m和m 用(B11)和(B12)式计算,k(q)为射线中心坐标系中的Gabor参数中的元素,有

可通过坐标转换矩阵 计算

主频率附近其他频率成分对应的高斯束中心射线路径和主频率成分对应的高斯束中心射线路径不同,但仍然可通过(B2)式利用主频率散射成分的中心射线计算,其中等号右端第二项说明了这种现象.更详细的解释参见Klimeš(2012)的工作.

附录C 高斯波包计算Gabor散射体产生的散射波相位部分

Gabor散射体的散射波利用高斯波包计算时,公式(20)中相位部分的计算可通过(C1)式实现

其中矩阵是相位函数对空间坐标的二阶导数,

Mij是相位函数在射线中心坐标系中第i和第j个方向的二阶混合偏导数.根据高斯波包的基础理论,可通过(C5)—(C9)式计算

其中

是辅助变量,

Vi是速度于射线中心坐标系中第i个坐标方向偏导数在点0处的值.相位函数的空间二阶导数在Gabor散射体中心位置的初始值为

把式(B5)、(B15)和(B16)代入式(C12)得到

把式(B15)和(B16)代入式(C13)得到

把式(B16)代入式(C14)得到

至此可利用(C1)式计算出高斯波包的相位部分,其中的矩阵通过公式(C5)—(C9)得到,初值为(C15)—(C17).

参考文献
[1] Beydoun W B, Mendes M. 1989. Elastic ray-Born l2-migration/inversion. Geophysical Journal International, 97(1):151-160.
[2] Bleistein N, Cohen J K. 1979. Direct inversion procedure for Claerbout's equations. Geophysics, 44(6):1034-1040.
[3] Bleistein N, Cohen J K, Hagin F G. 1984. Computational and asymptotic aspects of velocity inversion. Geophysics, 50(8):1253-1265.
[4] Bleistein N. 2008. Mathematics of modeling, migration and inversion with Gaussian beams. http://www.cwp.mines.edu/norm.
[5] Červený V, Popov M M, Pšenčík I. 1982. Computation of wave fields in inhomogeneous media-Gaussian beam approach. Geophysical Journal of the Royal Astronomical Society, 70(1):109-128.
[6] Červený V, Psencik I. 1983. Gaussian beams and paraxial ray approximation in three-dimensional elastic inhomogeneous media. Journal of Geophysics, 53:1-15.
[7] Červený V. 1992. Ray-Born synthetic seismograms for complex structures containing scatters. Journal of Seismic Exploration, 1:191-206.
[8] Červený V. 2001. Seismic Ray Theory. Cambridge:Cambridge University Press.
[9] Cheng G, Coen S. 1984. The relationship between Born inversion and migration for common-midpoint stacked data. Geophysics, 49(12):2117-2131.
[10] Klimeš L. 2004. Gaussian packets in smooth isotropic media. SW3D report 14, 43-54, Charles University.
[11] Klimeš L. 2012. Sensitivity of seismic waves to structure. Studia Geophysica et Geodaetica, 56(2):483-520.
[12] Li H, Feng B, Wang H Z. 2012. The wave-field modeling method with the integration of Gaussian packet. Geophysical Prospecting for Petroleum, 51(4):33-43.
[13] Li H, Wang H Z, Feng B, et al. 2014. Efficient pre-stack depth migration method using characteristic Gaussian packets. Chinese Journal of Geophysics, 57(7):2258-2268.
[14] Ma Y, Hale D. 2013. Wave-equation reflection traveltime inversion with dynamic warping and hybrid waveform inversion. 83rd SEG Annual International Meeting, Expanded Abstracts, 871-876.
[15] Mora P. 1989. Inversion=migration+tomography. Geophysics, 54(12):1575-1586.
[16] Popov M M. 1982. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4(1):85-97.
[17] Ralston J. 1983. Gaussian beams and the propagation of singularities:Studies in partial differential equations. MAA Studies in Mathematics, 1983, 23(1):206-248.
[18] Ren H. 2011. A study of the seismic wave inversion method for imaging in the acoustic media (in Chinese). Shanghai:Tongji University.
[19] Ren H R, Wang H Z, Chen S C. 2013. Least-squares reverse time migration in frequency domain using the adjoint-state method. Journal of Geophysics and Engineering, 10(3):035002, doi:10.1088/1742-2132/10/3/035002.
[20] Schneider W A. 1978. Integral formulation for migration in two and three dimensions. Geophysics, 43(1):49-76.
[21] Stolt R H. 1978. Migration by Fourier transform. Geophysics, 43(1):23-48.
[22] Xu S, Wang D, Chen F, et al. 2012. Full waveform inversion for reflected seismic data. 74th EAGE Annual International Meeting, Expanded Abstracts, W024.
[23] 李辉, 冯波, 王华忠. 2012. 波场模拟的高斯波包叠加方法. 石油物探, 51(4): 33-43.
[24] 李辉, 王华忠, 冯波等. 2014. 特征高斯波包叠前深度偏移方法. 地球物理学报, 57(7): 2258-2268.
[25] 任浩然. 2011. 声介质地震波反演成像方法研究[博士论文]. 上海: 同济大学.