基于Lagrangian Gaussian光束的数值方法

高景璐,郝永乐,孟品超

(1.吉林大学 数学学院,长春 130012;2.长春理工大学 理学院,长春 130022)

1 基于Green函数的Gaussian光束

在地震波模拟和高频状态迁移中,Gaussian光束累加方法应用广泛[1-10].对于标量波场U(x,ω),考虑如下Helmholtz波动方程:

2U(x,ω)+ω2n2(x)U(x,ω)=-δ(x-xe),

(1)

其中:x∈d(d=2,3);ω为频率;n(x)表示在x点的波速;xe为源点坐标.

对于U(x,ω),用标准几何光学高频率(即ω较大)拟设:

U(x,ω)=[A(x)+O(1/ω)]exp{iωτ(x)}.

(2)

将式(2)代入式(1),并令ω-2和ω-1为0,可得关于传播时间τ和振幅A且有相应初始条件的程函方程和迁移方程:

(τ)2=n2,

(3)

(4)

当ω足够大时,可以视Hellmholtz算子为椭圆型微分算子:

(5)

其中:p=(p1,p2,…,pn);x=(x1,x2,…,xn).

对于给定的光线γ,可以用文献[11-12]中的梁理论构建Gaussian射束,假设γ(s)={x(s)=x(s;xe,pe)},可得{p(s;xe,pe)}=τ(x(s)).对于如下射线追踪系统：

由梁理论[11-12],2τ的Hessian阵满足如下Riccati方程:

(10)

其中:

(11)

下面用迁移方程(4)得到光线上的振幅函数.由

方程(4)可写成

(12)

由方程(4)和上述结果可见,振幅永远不会无穷大,因此它也不会发散.

(13)

即定义了相位函数的一个全局逼近.类似地,可以定义O中沿中心线的振幅:

(14)

因此,单一光束的解可定义为

u(x;xe,pe)=A(x;xe,pe)exp{iωτ(x;xe,pe)}.

(15)

先对Helmholtz方程用Green函数,再对pe进行叠加,得

(16)

其中α(ω)为响应的某个常数[13].

2 2D波动方程的实现

下面考虑2D Helmholtz波动方程Lagrangian Gaussian光束的实现(d=2,x=(x,y),xe=(xe,ye)):

2U(x,y,ω)+ω2n2(x,y)U(x,y,ω) =-δ(x-xe)δ(y-ye),

(17)

2.1 预处理阶段

对于程函方程(3),取p=(p1,p2)=(ncosθ,nsinθ)(用θ代替p1,p2),通过式(6)～(9),射线追踪系统可化为

(18)

其中φe∈[0,2π]为相应于pe=((p1)e,(p2)e)的θ的初值.

进而,Riccati方程可化为

(19)

其中

D=-n2n-n·

为了保证非线性Riccati方程[12]有全局解,需对Hessian阵M(0)做如下约束:

1)M(0)=MT(0);

取

易见M(0)满足条件1);由条件2)等价于

(20)

选取

(21)

则其实数部分也满足条件2);对于条件3),有

a2sin2φe-2b2sinφecosφe+c2cos2φe≥0,

(22)

若取

(a2,b2,c2)=(sin2φe,-sinφecosφe,cos2φe),

(23)

则方程(22)和(20)的第二部分都成立.

从而迁移方程可写成：

(24)

预处理阶段步骤如下:

1) 离散计算区域,指标xi,yj,φm,sn分别为:

当n=1时初始化所有函数:

xs(n,m)=xe,ys(n,m)=ye,θ(n,m)=φm,τ(n,m)=0,A(n,m)=1,

其中xs(n,m)和ys(n,m)是光线上点的坐标.M(n,m)的初始值在式(21),(23)中都有说明.

2) 求解Liouville方程(18),(24).对于每个m=1,2,…,Nφ和n=1,2,…,Ns,定义xs(n,m),ys(n,m),θ(n,m),τ(n,m),M(n,m),A(n,m).

2.2 后置处理阶段

Lagrangian Gaussion光束(LGB)的计算效率由两部分组成:

1) 对于一个给定的光线,需找到一个合适的邻域,使得其满足对邻域中每个点得到的振幅无限小;

2) 对于1)中找到的给定点的邻域,需要在过其光线上找到离其最近的一个点.

对于1),文献[5]已给出了很多方法,本文给出一个更有效的新方法.图1显示了在光线的邻域内,只有一个集中域存在明显的振幅.表明要得到一个Gaussian射束,不需要计算所有的网格点.

(25)

图1 波导管Gaussian射束Fig.1 Waveguide Gaussian beam

图2 一小段光线γ的邻域Fig.2 Neighborhood of a segment of ray γ

图方法射线γ的邻域Fig.3 Neighborhood of ray γ method

图4 分割粗糙时射线γ的邻域Fig.4 Neighborhood of ray γ under coarse segmentation

下面讨论邻域的遍历方法,不失一般性,考虑射线γ斜率小于1的情况.

2.2.1 区域On,m的遍历算法

2) 对每个当前点yj,用y=yj分割clsn-1和clsn可得min_x和max_x,然后遍历xi从min_x到max_x.

2.2.2 寻找最近点的算法

2) 在当前区域On,m,已知csn-1和csn,运用插值公式

s(xi,yj)=csn-1+α1(csn-csn-1)

可以求出s(xi,yj),则对应于s(xi,yj)的点即为要寻找的点.

注3由于Δs相对于网格上的点较小,取coarse_N作为Δcs,同时要求Δcs的两个端点斜率相差不大(以coarse_N=10为例).在上述条件下,On,m即可视为一个矩形.根据该算法,只需要计算α1和s(xi,yj),避免了局域性的搜索,如图6所示.

图5 邻域的遍历算法Fig.5 Neighborhood traversal algorithms

图6 在γ上寻找最近点的方法Fig.6 Finding the nearest point on γ

2.2.3 后置处理过程

2) 对于(xi,yj),找到γm上关于s(xi,yj)的最近点;

3) 对于固定的点(xi,yj)和对应γm上的最近点,用式(13)～(15)得到一个Gaussian射束;

4) 将所有的Gaussian射束组合即构成一个完整的波域.

3 数值实例

3.1 常值模型

当v(x,y)=1时,渐近射线理论(ART)的解即为Hankel方程的渐近展开,它也是一点处Hankel方程的精确解.可以用这个解检验Eulerian Gaussian射束Eulerian冻结Gaussian射束的数值解.下面用这两种方法考察ω=16π时的情况.

3.2 波导模型

当速度函数为v(x,y)=3-2.5exp{-2x2}时,射线追踪显示在逼近源位置上该模型会发生焦散.传统的射线理论和Gaussian光束理论都预测了在焦散处有无穷次震动.因此可以用该模型检验Gaussian光束理论的正确性.

图7和图8分别为常值模型和波导模型的射线;图9和图10分别为ω=16π时常值模型和波导模型的数值结果.

图7 源为(xe,ye)=(0,0)时常值模型的射线Fig.7 Ray of constant model at source (xe,ye)=(0,0)

图8 源为(xe,ye)=(0,0)时波导模型的射线Fig.8 Ray of waveguide model at source (xe,ye)=(0,0)

(A) 真正的波场部分;(B) 在x=0处的切片;(C) 在y=0处的切片;(D) 在y=0.5处的切片.图9 ω=16π时常值模型的数值结果Fig.9 Numerical results for constant model when ω=16π

(A) 真正的波场部分;(B) 在x=0处的切片;(C) 在y=0处的切片;(D) 在y=0.5处的切片.图10 ω=16π时波导模型的数值结果Fig.10 Numerical results for waveguide model when ω=16π

[1] Alkhalifah T.Gaussian Beam Depth Migration for Anisotropic Media [J].Geophysics,1995,60(5):1474-1484.

[2] Babich V M,Buldyrev V S.Asymptotic Methods in Short-Wave Length Diffraction Problems [M].Moscow:Nauka,1972.

[4] Gray S.Gaussian Beam Migration of Common Shot Records [J].Geophysics,2005,70:133-136.

[5] Hill N R.Gaussian Beam Migration [J].Geophysics,1990,55(11):1416-1428.

[6] Hill N.Prestack Gaussian-Beam Depth Migration [J].Geophysics,2001,66:1240-1250.

[7] Norris A V,White B S,Schrieffer J R.Gaussian Wave Packets in Inhomogeneous Media with Curved Interfaces [J].Proc Roy Soc London,Ser:A,1987,412:93-123.

[8] White B S,Norris A N,Bayliss A,et al.Some Remarks on the Gaussian Beam Summation Method [J].Geophysical Journal of the Royal Astronomical Society,1987,89(2):579-636.

[9] LUO Song-ting,QIAN Jian-liang.Factored Singularities and High-Order Lax-Friedrichs Sweeping Schemes for Point-Source Traveltimes and Amplitudes [J].J Comput Phys,2011,230(12):4742-4755.

[10] LEUNG Shing-yu,QIAN Jian-liang.The Backward Phase Flow and FBI-Transform-Based Eulerian Gaussian Beams for the Schrödinger Equation [J].J Comput Phys,2010,229(23):8888-8917.

[11] Ralston J.Gaussian Beams and the Propagation of Singularities,Studies in Partial Differential Equations [J].MAA Studies in Mathematics,1982,23:206-248.

[12] Tanushev N M,QIAN Jian-liang,Ralston J.Mountain Waves and Gaussian Beams [J].Multiscale Model Simul,2007,6(2):688-709.

[13] LEUNG Shing-yu,QIAN Jian-liang,Burridge R.Eulerian Gaussian Beams for High Frequency Wave Propagation [J].Geophysics,2007,72(5):61-76.