基于全局笛卡爾坐標(biāo)系的高斯束地震波場(chǎng)模擬

白　敏　吳　娟*　孫章慶　楊　磊

(①華北水利水電大學(xué)資源與環(huán)境學(xué)院，河南鄭州 450046； ②吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春 130026； ③黃河水利委員會(huì)黃河水利科學(xué)研究院，河南鄭州 450003)

1　引言

地震波場(chǎng)模擬的目的是研究地震波在介質(zhì)中的傳播規(guī)律，為采集、處理和解釋提供有效的輔助依據(jù)。地震波場(chǎng)數(shù)值模擬方法主要包括射線類和波動(dòng)方程類。最常用的射線類方法是射線追蹤法，在高頻近似條件下，將波動(dòng)理論簡(jiǎn)化為射線理論，主要考慮地震波的運(yùn)動(dòng)學(xué)特征，計(jì)算速度快。然而射線類方法不適應(yīng)焦散區(qū)等復(fù)雜介質(zhì)情況。波動(dòng)方程類方法涉及豐富的波動(dòng)信息，精度高，但是相比射線類方法，計(jì)算效率較低。

高斯束是波動(dòng)方程集中于中心射線的高頻近似解，在與中心射線垂直的局部平面上，其振幅和相位都為高斯函數(shù)的形式，即隨著距中心射線距離的增加，振幅呈高斯衰減，相位呈拋物線近似特征。在常規(guī)的幾何射線理論中，振幅在焦散區(qū)奇異，而高斯束沿中心射線做復(fù)值展開，能正確處理焦散區(qū)和多值旅行時(shí)問(wèn)題[1-5]，同時(shí)具有較高精度，是常規(guī)射線類方法與波動(dòng)方程類方法之間的一種折衷。

業(yè)界也對(duì)高斯束正演模擬進(jìn)行了應(yīng)用研究，如：吳立明等[14]研究了非均勻復(fù)雜介質(zhì)高斯束正演模擬； 劉學(xué)才等[15]將高斯束應(yīng)用于三維VSP記錄的正演模擬； 張汝杰等[16]利用高斯束模擬井間地震； 鄧飛等[17]將高斯束正演推廣到三維； 孫成禹等[18]將高斯束正演與Zoeppritz方程結(jié)合，研究了能量約束下的聲波高斯束法地震波場(chǎng)正演。栗學(xué)磊[19]基于矢量波場(chǎng)疊加方法探討了彈性格林函數(shù)分解和矢量高斯波包法，并用高斯波包法合成地震記錄。劉鵬等[20]把簡(jiǎn)化的由相速度和群速度及其導(dǎo)數(shù)表示的運(yùn)動(dòng)學(xué)射線追蹤和動(dòng)力學(xué)射線追蹤系統(tǒng)用于高斯束方法，計(jì)算三維TTI介質(zhì)的格林函數(shù)。吳娟[21]基于已有的彈性高斯束理論，引入與品質(zhì)因子Q有關(guān)的復(fù)速度和精確的黏彈性Zoeppritz方程合成黏彈性地震記錄。

Leung等[28]詳細(xì)推導(dǎo)了全局笛卡爾坐標(biāo)系下求解高斯束的數(shù)學(xué)理論方法，并利用數(shù)學(xué)模型加以驗(yàn)證。本文首次嘗試將其理論應(yīng)用于全局笛卡爾坐標(biāo)系的高斯束地震波場(chǎng)模擬，并將模擬結(jié)果與射線中心坐標(biāo)系下的高斯束進(jìn)行比較，驗(yàn)證了方法的可行性。

2　方法原理

2.1　射線中心坐標(biāo)系下的高斯束

在射線中心坐標(biāo)系(圖1a)下求解高斯束。地下任一點(diǎn)P(x,z)的波場(chǎng)由距其最近的中心射線上的一點(diǎn)S(xm,zm)的波場(chǎng)值求得，假設(shè)S點(diǎn)的切線方向與z軸的夾角為θ，則射線中心坐標(biāo)系與笛卡爾坐標(biāo)系之間的坐標(biāo)變換公式為[23]

(1)

對(duì)于射線中心坐標(biāo)系下的高斯束，求解地下每一個(gè)點(diǎn)的波場(chǎng)值都要進(jìn)行坐標(biāo)變換。為了避免坐標(biāo)變換， Leung等[28]推導(dǎo)了在全局笛卡爾坐標(biāo)系(圖1b)下求解高斯束的數(shù)學(xué)理論方法，下面簡(jiǎn)要闡述其基本原理。

圖1　射線中心坐標(biāo)系(a)和全局笛卡爾坐標(biāo)系(b)

2.2　全局笛卡爾坐標(biāo)系下的高斯束

二維標(biāo)量波場(chǎng)的Helmholtz方程為

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

(2)

式中：ω為頻率；v為速度； 坐標(biāo)(xs,zs)為震源位置。當(dāng)ω很大時(shí)，把射線級(jí)數(shù)

(3)

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

(4a)

與傳輸方程

(4b)

程函方程是旅行時(shí)τ的函數(shù)，傳輸方程是振幅A的函數(shù)。假設(shè)射線只沿z軸向下傳播，即旅行時(shí)τ滿足

(5)

(6a)

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

(6b)

(6c)

式中τ0和ξ(x)是給定的復(fù)值光滑函數(shù)，滿足以下相容性條件

(7a)

(7b)

式中ξ1(x)是τ0在x方向的初值。給定震源處的初始條件

τ0(xs)=0

(8a)

(8b)

式中：θs為射線初始出射角；θmax為最大出射角。在震源鄰域構(gòu)建旅行時(shí)τ的近似函數(shù)

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

(9)

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

(10)

式中射線參數(shù)p(z)=τx[z,X(z)]。再構(gòu)建射線追蹤系統(tǒng)

(11a)

(11b)

(11c)

(12a)

(12b)

其中

(13a)

(13b)

(13c)

中心射線鄰域內(nèi)任一點(diǎn)X的復(fù)值旅行時(shí)為

τ(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方向的偏導(dǎo)數(shù)。式(12)可表示為

(16a)

(16b)

其中

(17a)

(17b)

(17c)

沿射線的振幅為

(18)

式中振幅處處非零?；谏厦娴耐茖?dǎo)，得到全局笛卡爾坐標(biāo)系下頻率域高斯束的表達(dá)式

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

(19)

式中τ(z,x;xs,θs)為復(fù)值旅行時(shí)(式(14))。

3　數(shù)值試驗(yàn)

本文主要研究全局笛卡爾坐標(biāo)系下的高斯束，為驗(yàn)證方法的可行性和有效性，分別通過(guò)一條高斯束、格林函數(shù)、連續(xù)介質(zhì)的單頻波場(chǎng)以及復(fù)雜介質(zhì)的波場(chǎng)進(jìn)行模擬，并對(duì)比射線中心坐標(biāo)系、全局笛卡爾坐標(biāo)系的高斯束精度。

3.1　一條高斯束

圖2為均勻介質(zhì)中一條高斯束在全局笛卡爾坐標(biāo)系、射線中心坐標(biāo)系的實(shí)部與虛部。由圖可見，全局笛卡爾坐標(biāo)系與射線中心坐標(biāo)系下高斯束的形態(tài)(束寬、走時(shí)等)是一致的。為驗(yàn)證本文方法的有效性，從圖2橫向距離1000m處各抽取一道，對(duì)全局笛卡爾坐標(biāo)系與射線中心坐標(biāo)系下一條高斯束的實(shí)部和虛部及其差值進(jìn)行對(duì)比(圖3)，結(jié)果表明兩者的精度很接近，且隨著傳播距離的增加兩者的差值逐漸減小。

3.2　高斯束表征的格林函數(shù)

高斯束波場(chǎng)傳播與偏移的核心是計(jì)算格林函數(shù)，因此格林函數(shù)的精度決定了波場(chǎng)的精度。在高斯束方法中，格林函數(shù)通過(guò)一系列由震源(zs,xs)出射且具有不同出射角θ的高斯束的疊加積分所表示，即

(20)

同解析格林函數(shù)相比，可用穩(wěn)相法求得權(quán)因子Φ[28]

(21)

圖4為均勻介質(zhì)格林函數(shù)在全局笛卡爾坐標(biāo)系、射線中心坐標(biāo)系的實(shí)部與虛部，圖5為圖4數(shù)據(jù)的精度分析結(jié)果。由圖可見，當(dāng)傳播距離大于200m時(shí)，由兩種坐標(biāo)系下的高斯束積分計(jì)算的格林函數(shù)均與格林函數(shù)解析解很接近，且相對(duì)誤差控制在2%以內(nèi)。

3.3　連續(xù)介質(zhì)單頻波場(chǎng)

利用Sinusoidal模型驗(yàn)證全局笛卡爾坐標(biāo)系高斯束對(duì)焦散區(qū)的處理效果，其速度函數(shù)為v(z,x)=1+0.2sin(0.5πz)sin[3π(x-0.55)](單位為km/s)。模型有兩個(gè)焦散區(qū)，用本文方法對(duì)其進(jìn)行射線追蹤。圖6為炮點(diǎn)位于1000、1500m處Sinusoidal模型的射線路徑及單頻波場(chǎng)。由圖可見，利用全局笛卡爾坐標(biāo)系下的高斯束可以很好地處理焦散問(wèn)題。

圖2　均勻介質(zhì)中一條高斯束在全局笛卡爾坐標(biāo)系(a)、射線中心坐標(biāo)系(b)的實(shí)部(左)與虛部(右)

圖3　圖2數(shù)據(jù)的精度分析結(jié)果(橫向距離1000m處)

3.4　復(fù)雜模型波場(chǎng)模擬

為檢驗(yàn)本文方法對(duì)復(fù)雜介質(zhì)的適用性，采用光滑的Marmousi模型進(jìn)行測(cè)試。通過(guò)疊加Marmousi模型(圖7a)的不同頻率的格林函數(shù)，然后進(jìn)行傅里葉逆變換得到兩個(gè)坐標(biāo)系下的波場(chǎng)快照(圖7b、圖7c)?？梢?，全局笛卡爾坐標(biāo)系高斯束較好地模擬了波的傳播特征，可以適應(yīng)較復(fù)雜的介質(zhì)(圖7b)，并且有望進(jìn)一步用于地震偏移成像。

3.5　算法精度分析

射線追蹤算法決定波場(chǎng)模擬精度。因此，文中主要從射線追蹤的角度分析全局笛卡爾坐標(biāo)系高斯束與射線中心坐標(biāo)系高斯束波場(chǎng)模擬的精度。圖8為兩種射線追蹤算法的射線路徑。由圖可見：全局笛卡爾坐標(biāo)系高斯束算法的射線追蹤范圍更小(圖8b左、圖8c左)，對(duì)于復(fù)雜介質(zhì)而言，射線的靈活性比射線中心坐標(biāo)系算法略低(箭頭標(biāo)示區(qū)域)。在射線追蹤過(guò)程中，隨著深度的增加，只有當(dāng)追蹤步長(zhǎng)越來(lái)越小時(shí)才能保證波場(chǎng)模擬精度，而全局笛卡爾坐標(biāo)系高斯束要求等步長(zhǎng)追蹤。因此，復(fù)雜介質(zhì)的全局笛卡爾坐標(biāo)系高斯束的精度會(huì)降低，橫向能量較弱。為此，Leung等[28]、Wu等[32]、Yang等[33]提出了針對(duì)全局笛卡爾坐標(biāo)系的歐式高斯束方法。

圖4　均勻介質(zhì)格林函數(shù)在全局笛卡爾坐標(biāo)系(a)、射線中心坐標(biāo)系(b)的實(shí)部(左)與虛部(右)

圖5　兩種方法計(jì)算的格林函數(shù)與解析解的精度對(duì)比

圖6　炮點(diǎn)位于1000m(a)、1500m(b)處Sinusoidal模型的射線路徑及單頻波場(chǎng)

圖7　Marmousi模型及其波場(chǎng)快照

圖8　兩種射線追蹤算法的射線路徑

(a)速度模型； (b)炮點(diǎn)位于橫向距離800m處的全局笛卡爾坐標(biāo)系高斯束算法(左)與射線中心坐標(biāo)系高斯束算法(右)； (c)炮點(diǎn)位于橫向距離1400m處的全局笛卡爾坐標(biāo)系高斯束算法(左)與射線中心坐標(biāo)系高斯束算法(右)

4　結(jié)束語(yǔ)

本文首次將Leung等[28]提出的全局笛卡爾坐標(biāo)系高斯束數(shù)學(xué)理論用于研究地震波傳播，理論分析和數(shù)值試驗(yàn)結(jié)果表明：在均勻介質(zhì)中，由全局笛卡爾坐標(biāo)系高斯束、射線中心坐標(biāo)系高斯束積分所計(jì)算的格林函數(shù)均很好地近似了解析格林函數(shù)；對(duì)連續(xù)介質(zhì)模型而言，全局笛卡爾坐標(biāo)系下的高斯束能很好地處理焦散問(wèn)題，較好地模擬了復(fù)雜介質(zhì)的波場(chǎng)傳播過(guò)程，為進(jìn)一步研究歐式高斯束提供了基礎(chǔ)。尚需指出，本文只是對(duì)全局笛卡爾坐標(biāo)系下的高斯束進(jìn)行了初步研究，全局笛卡爾坐標(biāo)系下的高斯束同樣可用于地震偏移成像。

感謝美國(guó)密歇根州立大學(xué)錢建良教授提供的代碼和幫助。

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

[4]劉強(qiáng),張敏,李振春等.各向異性介質(zhì)共炮域高斯束偏移.石油地球物理勘探,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]代福材,黃建平,李振春等.角度域黏聲介質(zhì)高斯束疊前深度偏移方法.石油地球物理勘探,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]吳立明,許云,烏達(dá)巴拉.高斯束射線法在二維非均勻介質(zhì)復(fù)雜構(gòu)造中的應(yīng)用.地球物理學(xué)報(bào),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]劉學(xué)才,周熙襄,沙椿.三維高斯射線束法合成三分量 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]張汝杰,賀振華,王理.井間地震高斯射線束正演方法.物探化探計(jì)算技術(shù),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]鄧飛,孫祥娥.三維高斯射線束正演的研究與實(shí)現(xiàn).大慶石油地質(zhì)與開發(fā),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]孫成禹,張文穎,倪長(zhǎng)寬等.能量約束下的高斯射線束法地震波場(chǎng)正演.石油地球物理勘探,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]栗學(xué)磊.基于高斯束的彈性波場(chǎng)計(jì)算方法研究[學(xué)位論文].吉林長(zhǎng)春:吉林大學(xué),2013.

[20]劉鵬,楊長(zhǎng)春,王彥飛.三維TTI介質(zhì)格林函數(shù)的高斯束計(jì)算方法.地球物理學(xué)進(jìn)展,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]吳娟.基于衰減補(bǔ)償?shù)母咚故品椒ㄑ芯縖學(xué)位論文].北京：中國(guó)石油大學(xué)(北京),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.