基于全局笛卡尔坐标系的高斯束地震波场模拟

白　敏　吴　娟*　孙章庆　杨　磊

(①华北水利水电大学资源与环境学院，河南郑州 450046； ②吉林大学地球探测科学与技术学院，吉林长春 130026； ③黄河水利委员会黄河水利科学研究院，河南郑州 450003)

1　引言

地震波场模拟的目的是研究地震波在介质中的传播规律，为采集、处理和解释提供有效的辅助依据。地震波场数值模拟方法主要包括射线类和波动方程类。最常用的射线类方法是射线追踪法，在高频近似条件下，将波动理论简化为射线理论，主要考虑地震波的运动学特征，计算速度快。然而射线类方法不适应焦散区等复杂介质情况。波动方程类方法涉及丰富的波动信息，精度高，但是相比射线类方法，计算效率较低。

高斯束是波动方程集中于中心射线的高频近似解，在与中心射线垂直的局部平面上，其振幅和相位都为高斯函数的形式，即随着距中心射线距离的增加，振幅呈高斯衰减，相位呈抛物线近似特征。在常规的几何射线理论中，振幅在焦散区奇异，而高斯束沿中心射线做复值展开，能正确处理焦散区和多值旅行时问题[1-5]，同时具有较高精度，是常规射线类方法与波动方程类方法之间的一种折衷。

业界也对高斯束正演模拟进行了应用研究，如：吴立明等[14]研究了非均匀复杂介质高斯束正演模拟； 刘学才等[15]将高斯束应用于三维VSP记录的正演模拟； 张汝杰等[16]利用高斯束模拟井间地震； 邓飞等[17]将高斯束正演推广到三维； 孙成禹等[18]将高斯束正演与Zoeppritz方程结合，研究了能量约束下的声波高斯束法地震波场正演。栗学磊[19]基于矢量波场叠加方法探讨了弹性格林函数分解和矢量高斯波包法，并用高斯波包法合成地震记录。刘鹏等[20]把简化的由相速度和群速度及其导数表示的运动学射线追踪和动力学射线追踪系统用于高斯束方法，计算三维TTI介质的格林函数。吴娟[21]基于已有的弹性高斯束理论，引入与品质因子Q有关的复速度和精确的黏弹性Zoeppritz方程合成黏弹性地震记录。

Leung等[28]详细推导了全局笛卡尔坐标系下求解高斯束的数学理论方法，并利用数学模型加以验证。本文首次尝试将其理论应用于全局笛卡尔坐标系的高斯束地震波场模拟，并将模拟结果与射线中心坐标系下的高斯束进行比较，验证了方法的可行性。

2　方法原理

2.1　射线中心坐标系下的高斯束

在射线中心坐标系(图1a)下求解高斯束。地下任一点P(x,z)的波场由距其最近的中心射线上的一点S(xm,zm)的波场值求得，假设S点的切线方向与z轴的夹角为θ，则射线中心坐标系与笛卡尔坐标系之间的坐标变换公式为[23]

(1)

对于射线中心坐标系下的高斯束，求解地下每一个点的波场值都要进行坐标变换。为了避免坐标变换， Leung等[28]推导了在全局笛卡尔坐标系(图1b)下求解高斯束的数学理论方法，下面简要阐述其基本原理。

图1　射线中心坐标系(a)和全局笛卡尔坐标系(b)

2.2　全局笛卡尔坐标系下的高斯束

二维标量波场的Helmholtz方程为

=-δ(x-xs)δ(z-zs)

(2)

式中：ω为频率；v为速度； 坐标(xs,zs)为震源位置。当ω很大时，把射线级数

(3)

代入式(2)，得到程函方程

(4a)

与传输方程

(4b)

程函方程是旅行时τ的函数，传输方程是振幅A的函数。假设射线只沿z轴向下传播，即旅行时τ满足

(5)

(6a)

τ(x,0)=τ0(x)Im(τ0)≥0

(6b)

(6c)

式中τ0和ξ(x)是给定的复值光滑函数，满足以下相容性条件

(7a)

(7b)

式中ξ1(x)是τ0在x方向的初值。给定震源处的初始条件

τ0(xs)=0

(8a)

(8b)

式中：θs为射线初始出射角；θmax为最大出射角。在震源邻域构建旅行时τ的近似函数

τ0(x;xs)=τ0(xs)+ξ1(xs;θs)(x-xs)+

(9)

式中ε>0，一般取值为1。根据Ralston[26]、Tanushev等[27]提出的高斯束理论，令x=X(z)，τ=T(z)，然后引入Hamiltonian算子

(10)

式中射线参数p(z)=τx[z,X(z)]。再构建射线追踪系统

(11a)

(11b)

(11c)

(12a)

(12b)

其中

(13a)

(13b)

(13c)

中心射线邻域内任一点X的复值旅行时为

τ(z,x;xs,θs)=T(z;xs,θs)+p(z)·[x-X(z)]+

(14)

式中：τx=p；τxx=δp/δx=(∂p/∂α)/(∂x/∂α)=B/C。

(15a)

(15b)

(15c)

式中vx、vz分别为v沿x、z方向的偏导数。式(12)可表示为

(16a)

(16b)

其中

(17a)

(17b)

(17c)

沿射线的振幅为

(18)

式中振幅处处非零。基于上面的推导，得到全局笛卡尔坐标系下频率域高斯束的表达式

uGB=A(z;xs,θs)exp[iωτ(z,x;xs,θs)]

(19)

式中τ(z,x;xs,θs)为复值旅行时(式(14))。

3　数值试验

本文主要研究全局笛卡尔坐标系下的高斯束，为验证方法的可行性和有效性，分别通过一条高斯束、格林函数、连续介质的单频波场以及复杂介质的波场进行模拟，并对比射线中心坐标系、全局笛卡尔坐标系的高斯束精度。

3.1　一条高斯束

图2为均匀介质中一条高斯束在全局笛卡尔坐标系、射线中心坐标系的实部与虚部。由图可见，全局笛卡尔坐标系与射线中心坐标系下高斯束的形态(束宽、走时等)是一致的。为验证本文方法的有效性，从图2横向距离1000m处各抽取一道，对全局笛卡尔坐标系与射线中心坐标系下一条高斯束的实部和虚部及其差值进行对比(图3)，结果表明两者的精度很接近，且随着传播距离的增加两者的差值逐渐减小。

3.2　高斯束表征的格林函数

高斯束波场传播与偏移的核心是计算格林函数，因此格林函数的精度决定了波场的精度。在高斯束方法中，格林函数通过一系列由震源(zs,xs)出射且具有不同出射角θ的高斯束的叠加积分所表示，即

(20)

同解析格林函数相比，可用稳相法求得权因子Φ[28]

(21)

图4为均匀介质格林函数在全局笛卡尔坐标系、射线中心坐标系的实部与虚部，图5为图4数据的精度分析结果。由图可见，当传播距离大于200m时，由两种坐标系下的高斯束积分计算的格林函数均与格林函数解析解很接近，且相对误差控制在2%以内。

3.3　连续介质单频波场

利用Sinusoidal模型验证全局笛卡尔坐标系高斯束对焦散区的处理效果，其速度函数为v(z,x)=1+0.2sin(0.5πz)sin[3π(x-0.55)](单位为km/s)。模型有两个焦散区，用本文方法对其进行射线追踪。图6为炮点位于1000、1500m处Sinusoidal模型的射线路径及单频波场。由图可见，利用全局笛卡尔坐标系下的高斯束可以很好地处理焦散问题。

图2　均匀介质中一条高斯束在全局笛卡尔坐标系(a)、射线中心坐标系(b)的实部(左)与虚部(右)

图3　图2数据的精度分析结果(横向距离1000m处)

3.4　复杂模型波场模拟

为检验本文方法对复杂介质的适用性，采用光滑的Marmousi模型进行测试。通过叠加Marmousi模型(图7a)的不同频率的格林函数，然后进行傅里叶逆变换得到两个坐标系下的波场快照(图7b、图7c)。可见，全局笛卡尔坐标系高斯束较好地模拟了波的传播特征，可以适应较复杂的介质(图7b)，并且有望进一步用于地震偏移成像。

3.5　算法精度分析

射线追踪算法决定波场模拟精度。因此，文中主要从射线追踪的角度分析全局笛卡尔坐标系高斯束与射线中心坐标系高斯束波场模拟的精度。图8为两种射线追踪算法的射线路径。由图可见：全局笛卡尔坐标系高斯束算法的射线追踪范围更小(图8b左、图8c左)，对于复杂介质而言，射线的灵活性比射线中心坐标系算法略低(箭头标示区域)。在射线追踪过程中，随着深度的增加，只有当追踪步长越来越小时才能保证波场模拟精度，而全局笛卡尔坐标系高斯束要求等步长追踪。因此，复杂介质的全局笛卡尔坐标系高斯束的精度会降低，横向能量较弱。为此，Leung等[28]、Wu等[32]、Yang等[33]提出了针对全局笛卡尔坐标系的欧式高斯束方法。

图4　均匀介质格林函数在全局笛卡尔坐标系(a)、射线中心坐标系(b)的实部(左)与虚部(右)

图5　两种方法计算的格林函数与解析解的精度对比

图6　炮点位于1000m(a)、1500m(b)处Sinusoidal模型的射线路径及单频波场

图7　Marmousi模型及其波场快照

图8　两种射线追踪算法的射线路径

(a)速度模型； (b)炮点位于横向距离800m处的全局笛卡尔坐标系高斯束算法(左)与射线中心坐标系高斯束算法(右)； (c)炮点位于横向距离1400m处的全局笛卡尔坐标系高斯束算法(左)与射线中心坐标系高斯束算法(右)

4　结束语

本文首次将Leung等[28]提出的全局笛卡尔坐标系高斯束数学理论用于研究地震波传播，理论分析和数值试验结果表明：在均匀介质中，由全局笛卡尔坐标系高斯束、射线中心坐标系高斯束积分所计算的格林函数均很好地近似了解析格林函数；对连续介质模型而言，全局笛卡尔坐标系下的高斯束能很好地处理焦散问题，较好地模拟了复杂介质的波场传播过程，为进一步研究欧式高斯束提供了基础。尚需指出，本文只是对全局笛卡尔坐标系下的高斯束进行了初步研究，全局笛卡尔坐标系下的高斯束同样可用于地震偏移成像。

感谢美国密歇根州立大学钱建良教授提供的代码和帮助。

[3]Popov M M.Ray Theory and Gaussian Beam Method for Geophysicists.Edufba,2002.

[4]刘强,张敏,李振春等.各向异性介质共炮域高斯束偏移.石油地球物理勘探,2016,51(5):930-937.

Liu Qiang,Zhang Min,Li Zhenchun et al.Common-shot domain Gaussian beam migration in anisotropic media.OGP,2016,51(5):930-937.

[5]代福材,黄建平,李振春等.角度域黏声介质高斯束叠前深度偏移方法.石油地球物理勘探,2017,52(2):283-293.

Dai Fucai,Huang Jianping,Li Zhenchun et al.Angle domain prestack Gaussian beam migration for visco-acoustic media.OGP,2017,52(2):283-293.

[9]Nowack R,Aki K.The two-dimensional Gaussian beam synthetic method:testing and application.Journal of Geophysical Research,1984,89:7797-7819.

[10]Popov M M.A method of summation of Gaussian beams in isotropic theory of elasticity.Physics of the Solid Earth,1984,19:699-706.

[11]Raz S.Beam stacking:A generalized preprocessing technique.Geophysics,1987,52(9):1199-1210.

[12]Hanyga A.Gaussian beams in anisotropic elastic media.Geophysical Journal International,1986,85(3):473-504.

[13]Liu Y X.Gaussian-beam Modeling in Heterogeneous Anisotropic Media[D].Colorado School of Mines,2010.

[14]吴立明,许云,乌达巴拉.高斯束射线法在二维非均匀介质复杂构造中的应用.地球物理学报,1995,38(增刊1):144-152.

Wu Liming,Xu Yun,Wudabala.Applied study of Gaussian beam method in 2-D inhomogeneous media and laterally varying layered structures.Chinese Journal of Geophysics,1995,38(S1):144-152.

[15]刘学才,周熙襄,沙椿.三维高斯射线束法合成三分量 VSP 记录.石油地球物理勘探,1995,30(5):669-680.

Liu Xuecai,Zhou Xixiang,Sha Chun.Synthesizing three-component VSP data by 3-D Gaussian beam method.OGP,1995,30(5):669-680.

[16]张汝杰,贺振华,王理.井间地震高斯射线束正演方法.物探化探计算技术,1997,19(2):128-137.

Zhang Rujie,He Zhenhua,Wang Li.Forward modeling by Gaussian beam method in crosswell seismic.Computing Techniques for Geophysical and Geochemical Exploration,1997,19(2):128-137.

[17]邓飞,孙祥娥.三维高斯射线束正演的研究与实现.大庆石油地质与开发,2009,28(1):127-131.

Deng Fei,Sun Xiang’e.Study and realization of 3D Gaussian beam forward modeling.Petroleum Geology and Oilfield Development in Daqing,2009,28(1):127-131.

[18]孙成禹,张文颖,倪长宽等.能量约束下的高斯射线束法地震波场正演.石油地球物理勘探,2011,46(6):856-861.

Sun Chengyu,Zhang Wenying,Ni Changkuan et al.Seismic wave field modeling by energy controlled Gaussian beam method.OGP,2011,46(6):856-861.

[19]栗学磊.基于高斯束的弹性波场计算方法研究[学位论文].吉林长春:吉林大学,2013.

[20]刘鹏,杨长春,王彦飞.三维TTI介质格林函数的高斯束计算方法.地球物理学进展,2013,28(5):2547-2553.

Liu Peng,Yang Changchun,Wang Yanfei.Three-dimensional Green’s function calculation in TTI media using the Gaussian beam method.Progress in Geophysics,2013,28(5):2547-2553.

[21]吴娟.基于衰减补偿的高斯束偏移方法研究[学位论文].北京：中国石油大学(北京),2015.

[22]George T,Virieux J,Madariaga R.Seismic wave synthesis by Gaussian beam summation:A comparison with finite differences.Geophysics,1987,52(8):1065-1073.

[23]Hill N R.Gaussian beam migration.Geophysics,1990,55(11):1416-1428.

[26]Ralston J.Gaussian beams and the propagation of singularities.Mathematics Association of America Studies in Mathematics,1983,23:206-248.

[27]Tanushev N M,Qian J,Ralston J V.Mountain waves and Gaussian beams.Multiscale Modeling & Simulation,2007,6(2):688-709.

[28]Leung S,Qian J,Burridge R.Eulerian Gaussian beams for high-frequency wave propagation.Geophysics,2007,72(5):SM61-SM76.

[29]Gray S and May W.Kirchhoff migration using eikonal equation traveltimes.Geophysics,1994,59(5):810-817.

[30]Qian J and Symes W W.An adaptive finite difference method for traveltime and amplitude.Geophysics,2002,67(1):167-176.

[31]Symes W W and Qian J.A slowness matching Eulerian method for multivalued solutions of eikonal equations.Journal of Scientific Computing.2003,19(1-3):501-526.

[32]Wu H,Yang X.The Eulerian Gaussian beam method for high frequency wave propagation in the reduced momentum space.Wave Motion,2013,50(6):1036-1049.

[33]Yang Xu,Lu Jianfeng,Fomel S.Seismic modeling using the frozen Gaussian approximation.SEG Technical Program Expanded Abstracts,2013,32:52-58.