基于Lagrangian Gaussian光束的數(shù)值方法

高景璐,郝永樂,孟品超

(1.吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012;2.長春理工大學(xué) 理學(xué)院,長春 130022)

1 基于Green函數(shù)的Gaussian光束

在地震波模擬和高頻狀態(tài)遷移中,Gaussian光束累加方法應(yīng)用廣泛[1-10].對(duì)于標(biāo)量波場(chǎng)U(x,ω),考慮如下Helmholtz波動(dòng)方程:

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

(1)

其中:x∈d(d=2,3);ω為頻率;n(x)表示在x點(diǎn)的波速;xe為源點(diǎn)坐標(biāo).

對(duì)于U(x,ω),用標(biāo)準(zhǔn)幾何光學(xué)高頻率(即ω較大)擬設(shè):

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

(2)

將式(2)代入式(1),并令ω-2和ω-1為0,可得關(guān)于傳播時(shí)間τ和振幅A且有相應(yīng)初始條件的程函方程和遷移方程:

(τ)2=n2,

(3)

(4)

當(dāng)ω足夠大時(shí),可以視Hellmholtz算子為橢圓型微分算子:

(5)

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

對(duì)于給定的光線γ,可以用文獻(xiàn)[11-12]中的梁理論構(gòu)建Gaussian射束,假設(shè)γ(s)={x(s)=x(s;xe,pe)},可得{p(s;xe,pe)}=τ(x(s)).對(duì)于如下射線追蹤系統(tǒng)：

由梁理論[11-12],2τ的Hessian陣滿足如下Riccati方程:

(10)

其中:

(11)

下面用遷移方程(4)得到光線上的振幅函數(shù).由

方程(4)可寫成

(12)

由方程(4)和上述結(jié)果可見,振幅永遠(yuǎn)不會(huì)無窮大,因此它也不會(huì)發(fā)散.

(13)

即定義了相位函數(shù)的一個(gè)全局逼近.類似地,可以定義O中沿中心線的振幅:

(14)

因此,單一光束的解可定義為

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

(15)

先對(duì)Helmholtz方程用Green函數(shù),再對(duì)pe進(jìn)行疊加,得

(16)

其中α(ω)為響應(yīng)的某個(gè)常數(shù)[13].

2 2D波動(dòng)方程的實(shí)現(xiàn)

下面考慮2D Helmholtz波動(dòng)方程Lagrangian Gaussian光束的實(shí)現(xiàn)(d=2,x=(x,y),xe=(xe,ye)):

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

(17)

2.1 預(yù)處理階段

對(duì)于程函方程(3),取p=(p1,p2)=(ncosθ,nsinθ)(用θ代替p1,p2),通過式(6)～(9),射線追蹤系統(tǒng)可化為

(18)

其中φe∈[0,2π]為相應(yīng)于pe=((p1)e,(p2)e)的θ的初值.

進(jìn)而,Riccati方程可化為

(19)

其中

D=-n2n-n·

為了保證非線性Riccati方程[12]有全局解,需對(duì)Hessian陣M(0)做如下約束:

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

取

易見M(0)滿足條件1);由條件2)等價(jià)于

(20)

選取

(21)

則其實(shí)數(shù)部分也滿足條件2);對(duì)于條件3),有

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

(22)

若取

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

(23)

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

從而遷移方程可寫成：

(24)

預(yù)處理階段步驟如下:

1) 離散計(jì)算區(qū)域,指標(biāo)xi,yj,φm,sn分別為:

當(dāng)n=1時(shí)初始化所有函數(shù):

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

其中xs(n,m)和ys(n,m)是光線上點(diǎn)的坐標(biāo).M(n,m)的初始值在式(21),(23)中都有說明.

2) 求解Liouville方程(18),(24).對(duì)于每個(gè)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)的計(jì)算效率由兩部分組成:

1) 對(duì)于一個(gè)給定的光線,需找到一個(gè)合適的鄰域,使得其滿足對(duì)鄰域中每個(gè)點(diǎn)得到的振幅無限小;

2) 對(duì)于1)中找到的給定點(diǎn)的鄰域,需要在過其光線上找到離其最近的一個(gè)點(diǎn).

對(duì)于1),文獻(xiàn)[5]已給出了很多方法,本文給出一個(gè)更有效的新方法.圖1顯示了在光線的鄰域內(nèi),只有一個(gè)集中域存在明顯的振幅.表明要得到一個(gè)Gaussian射束,不需要計(jì)算所有的網(wǎng)格點(diǎn).

(25)

圖1 波導(dǎo)管Gaussian射束Fig.1 Waveguide Gaussian beam

圖2 一小段光線γ的鄰域Fig.2 Neighborhood of a segment of ray γ

圖方法射線γ的鄰域Fig.3 Neighborhood of ray γ method

圖4 分割粗糙時(shí)射線γ的鄰域Fig.4 Neighborhood of ray γ under coarse segmentation

下面討論鄰域的遍歷方法,不失一般性,考慮射線γ斜率小于1的情況.

2.2.1 區(qū)域On,m的遍歷算法

2) 對(duì)每個(gè)當(dāng)前點(diǎn)yj,用y=yj分割clsn-1和clsn可得min_x和max_x,然后遍歷xi從min_x到max_x.

2.2.2 尋找最近點(diǎn)的算法

2) 在當(dāng)前區(qū)域On,m,已知csn-1和csn,運(yùn)用插值公式

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

可以求出s(xi,yj),則對(duì)應(yīng)于s(xi,yj)的點(diǎn)即為要尋找的點(diǎn).

注3由于Δs相對(duì)于網(wǎng)格上的點(diǎn)較小,取coarse_N作為Δcs,同時(shí)要求Δcs的兩個(gè)端點(diǎn)斜率相差不大(以coarse_N=10為例).在上述條件下,On,m即可視為一個(gè)矩形.根據(jù)該算法,只需要計(jì)算α1和s(xi,yj),避免了局域性的搜索,如圖6所示.

圖5 鄰域的遍歷算法Fig.5 Neighborhood traversal algorithms

圖6 在γ上尋找最近點(diǎn)的方法Fig.6 Finding the nearest point on γ

2.2.3 后置處理過程

2) 對(duì)于(xi,yj),找到γm上關(guān)于s(xi,yj)的最近點(diǎn);

3) 對(duì)于固定的點(diǎn)(xi,yj)和對(duì)應(yīng)γm上的最近點(diǎn),用式(13)～(15)得到一個(gè)Gaussian射束;

4) 將所有的Gaussian射束組合即構(gòu)成一個(gè)完整的波域.

3 數(shù)值實(shí)例

3.1 常值模型

當(dāng)v(x,y)=1時(shí),漸近射線理論(ART)的解即為Hankel方程的漸近展開,它也是一點(diǎn)處Hankel方程的精確解.可以用這個(gè)解檢驗(yàn)Eulerian Gaussian射束Eulerian凍結(jié)Gaussian射束的數(shù)值解.下面用這兩種方法考察ω=16π時(shí)的情況.

3.2 波導(dǎo)模型

當(dāng)速度函數(shù)為v(x,y)=3-2.5exp{-2x2}時(shí),射線追蹤顯示在逼近源位置上該模型會(huì)發(fā)生焦散.傳統(tǒng)的射線理論和Gaussian光束理論都預(yù)測(cè)了在焦散處有無窮次震動(dòng).因此可以用該模型檢驗(yàn)Gaussian光束理論的正確性.

圖7和圖8分別為常值模型和波導(dǎo)模型的射線;圖9和圖10分別為ω=16π時(shí)常值模型和波導(dǎo)模型的數(shù)值結(jié)果.

圖7 源為(xe,ye)=(0,0)時(shí)常值模型的射線Fig.7 Ray of constant model at source (xe,ye)=(0,0)

圖8 源為(xe,ye)=(0,0)時(shí)波導(dǎo)模型的射線Fig.8 Ray of waveguide model at source (xe,ye)=(0,0)

(A) 真正的波場(chǎng)部分;(B) 在x=0處的切片;(C) 在y=0處的切片;(D) 在y=0.5處的切片.圖9 ω=16π時(shí)常值模型的數(shù)值結(jié)果Fig.9 Numerical results for constant model when ω=16π

(A) 真正的波場(chǎng)部分;(B) 在x=0處的切片;(C) 在y=0處的切片;(D) 在y=0.5處的切片.圖10 ω=16π時(shí)波導(dǎo)模型的數(shù)值結(jié)果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.