已知背景模型与其中传播的背景波场,同时给出散射体(模型扰动部分),如何计算散射体产生的散射波场是一个有应用价值的工作.散射波的正、反演与偏移成像有密切的联系.Cheng和Coen(1984)详细推导并分析了常速介质中散射波反演与F-K偏移(Stolt,1978)的关系,Bleistein等(1984)指出WKBJ近似下散射波反演的一次迭代与Kirchhoff偏移(Schneider,1978)是等价的.最小二乘偏移即从散射波反演中发展而来,其中的正演过程即为散射波模拟(Beydoun and Mendes, 1989;任浩然,2011).反射全波形反演(RFWI)的研究方兴未艾(Xu et al., 2012;Ma and Hale, 2013),其中关键步骤是利用反偏移计算目标函数的梯度,反偏移在背景介质中利用模型高波数成分(即散射体)模拟散射波.事实上,Mora(1989)在1980年代末就已经分析了反演可分为“偏移”和“层析”两部分,其中层析部分的模拟过程为散射波模拟.从上述分析可知,在勘探地震领域,数据处理的核心部分——偏移与反演均与散射波的正演有着密不可分的关系.所以散射波正演理论研究对勘探地震的发展从基础上存在显著的推动作用.鉴于散射波模拟有重要的理论和应用价值,本文致力于在Gabor域描述散射体时散射波模拟方法的研究.
在网格中描述散射体的空间分布,每一个散射体单元是一个网格点,在此基础上,传统散射波应用(最小二乘偏移等)中需要每个网格点的散射波,实际计算时应用格林函数计算每个检波点观测到的散射波(任浩然,2011).计算格林函数的方法有射线法(Bleistein and Cohen, 1979;Bleistein et al.,1984;Červený,1992)、波动方程法(任浩然,2011;Ren et al., 2013)和高斯束/高斯波包法(Červený and Psencik,1983;李辉等,2012)等.射线方法效率较高,但受到射线理论的约束,无法处理正演时的焦散及阴影区问题(;Červený,1992),且难以计算波场中的子波相位信息;波动方程的计算效率特别低;高斯束/高斯波包计算格林函数时需要对不同初始角度的高斯束/高斯波包进行积分,体现不出高斯束/高斯波包的效率优势.总之,在传统网格描述散射体的基础上,通过格林函数计算散射波的方法存在效率低等诸多问题,同时由于网格描述散射体时冗余性较高,进一步降低了散射波模拟的效率.
在Gabor域描述散射体时,每一个散射体单元变成一个Gabor函数,一个散射体单元散射出的散射波不再是格林函数.在一阶Born近似和零阶WKBJ近似假设下,Klimeš(2012)推导出Gabor函数散射体产生的散射波场可以用高斯束(频率域)或高斯波包(时间域)计算,在各向异性弹性介质中(模型参数为刚度系数和密度),一个Gabor函数表达的散射体单元所产生的散射波可利用最多5个高斯波包(时间域)或高斯束(频率域)描述.高斯束( ;Červený,et al, 1982;Popov,1982)和高斯波包(Ralston,1983;Klimeš 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的计算写成如下形式:



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

′点的振幅和走时.三维介质中格林函数的零阶WKBJ形式为(Bleistein,2008)

,
′)是从点
′到
的走时.
![]() |
图 3 二维空间中的Gabor函数示意图Fig. 3 A Gabor function in 2-dimensions |
将散射体的Gabor函数表达式(3),背景波场与格林函数的零阶WKBJ近似(4)和(5)代入一阶Born近似下的散射波场积分公式(1)中,得到Gabor散射体的散射波积分形式(Klimeš, 2012)


0,Δ
是
0到主散射射线附近点的矢量.公式(8)中变量在Klimeš(2012)的文章中有详细的解释,为了保持本文理论的完整性,我们把上述变量的物理意义表述在附录A中.
0位于主散射射线上,当Δ
=
时,散射波us(
) 的振幅最大,即式(8)中振幅衰减最小,此时有



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

α垂直时,(12)式不存在关于1/ω的非平凡解,此时不产生散射波,否则(12)式存在唯一关于1/ω的非平凡解,此时散射波的主频为

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

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





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散射体产生的散射波场在时间域的表达形式,得到时间域计算散射波场的方式——高斯波包.高斯波包计算散射波场写成振幅和复数相位的表达形式:

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

α,t)是高斯波包初始振幅,Q1和Q 2是公式(B4)所示动力学射线追踪传播矩阵的元素,MGPα是高斯波包相位函数在射线中心坐标系中q方向的二阶导数初值,详见(C12)式.高斯波包(20)式中振幅部分(21)式的计算可通过高斯束的频率积分得到,积分结果的振幅在散射点的表达式即为高斯波包的初始振幅,具体形式为

α,ω0)是频率为ω0的高斯束初始振幅,M44α是高斯束相位函数在频率方向的二阶导数,参见附录B中的(B16)式.结合高斯束初始振幅(19)、(B16)和(22)式,高斯波包的初始振幅可表达如下:

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



至此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)



本节在二维介质中对上述理论进行数值实验.有限差分法(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 |
针对一个理论模型,利用高斯波包模拟散射波的理论,结合公式(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). |
速度模型可分成背景部分和扰动部分(散射体),波场在背景模型中传播遇到散射体时产生散射波场.当散射体在空间的分布是一个Gabor函数时,相应的散射波可用高斯束(频率域)或高斯波包(时间域)描述.这就给散射波的模拟与偏移、层析等应用提供了新的思路——在Gabor变换域以Gabor函数为框架描述散射体,散射波的正演工具为高斯束和高斯波包.
本文在常密度声波假设下推导了利用高斯束和高斯波包计算Gabor散射体产生散射波的具体方法.与刚度系数和密度为模型参数的弹性理论下的散射波模拟方法相比,常密度声波方程中散射波模拟方法的主要不同之处在于散射高斯束(或高斯波包)的个数以及高斯束(或高斯波包)的初始振幅.此外,三维介质和二维介质中的高斯束(或高斯波包)初始振幅也有差异.常密度声波假设下高斯束与高斯波包初始振幅在三维和二维介质中的计算方法是本文的主要理论工作.
在上述理论方法的基础上,本文在时间域利用高斯波包算子实现了一个散射体单元(Gabor函数)产生的散射波的计算.计算结果与FD对比,两者在时间域非常相似,频率域中高斯波包法与FD计算的结果稍有差异,前者是一个单一的窄带信号,后者是若干个窄带信号,但后者的主要能量仍然集中于一个窄带信号中,文中对此现象做了具体分析.
同一高斯波包描述的散射体,背景波场入射方向不同时,高斯波包计算散射波的精度不同,文中实验对此做了精度分析,入射角度越小精度越高.总体误差低于9%,可见本文推导的振幅计算公式精度较高.
多Gabor函数表达的散射体所产生的散射波实验中给出了计算散射波的具体策略.从对水平层散射体的高斯波包散射波模拟实验中,通过分析可知哪些局部散射体对观测数据有影响,哪些局部散射体对观测数据没有影响,这点在偏移和层析等应用中很重要.
针对理论速度模型,模拟了相应的地表观测数据.本文方法与FD的对比说明本文方法可以保证模拟波场的精度.但在大偏移距部分稍有差异,这与前文中 “大角度模拟精度小于小角度”的结论一致.
高斯波包计算散射波的方法主要基于如下几个假设:零阶WKBJ近似、一阶Born近似和入射波无限宽频.这些假设限制了方法的精度.但由于高斯波包方法效率高、应用灵活、易于处理局部数据,所以本文散射波模拟方法在偏移、反演中有一定的应用前景.
致谢 感谢中国石化集团公司科技部及中国石化石油物探技术研究院(南京)对WPI研究组及本项研究的支持. 附录A 公式(8)中变量的物理意义
公式(4)所示背景波场的走时T(
′)在Gabor散射体扰动的中心点
α的值记为T(
α). T(
α)在点
α附近的一阶和二阶微分简写为

,
′),
′和
分别是格林函数的源点和观测点.τ(
,
′)在点
和
′附近的一阶和二阶微分简写成

散射波场在主散射路径上的点为
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附近做泰勒展开,有

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

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

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



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

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

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

(q)是Gabor函数参数 Kα在射线中心坐标系中的表达形式,

计算得到
(q)



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


中的元素,有

计算

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

是相位函数对空间坐标的二阶导数,





0处的值.相位函数的空间二阶导数在Gabor散射体中心位置的初始值为




至此可利用(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. 声介质地震波反演成像方法研究[博士论文]. 上海: 同济大学. |
2015, Vol. 58















