白　敏　吳　娟*　孫章慶　楊　磊

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

1　引言

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

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

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

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

2　方法原理

2.1　射線中心坐標系下的高斯束

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

(1)

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

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

2.2　全局笛卡爾坐標系下的高斯束

二維標量波場的Helmholtz方程為

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

(2)

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

(3)

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

(4a)

與傳輸方程

(4b)

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

(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ù)

τ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)任一點X的復(fù)值旅行時為

τ(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)，得到全局笛卡爾坐標系下頻率域高斯束的表達式

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

(19)

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

3　數(shù)值試驗

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

3.1　一條高斯束

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

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

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

(20)

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

(21)

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

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

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

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

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

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

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

3.5　算法精度分析

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

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

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

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

圖7　Marmousi模型及其波場快照

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

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

4　結(jié)束語

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

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

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

[4]劉強,張敏,李振春等.各向異性介質(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]吳立明,許云,烏達巴拉.高斯束射線法在二維非均勻介質(zhì)復(fù)雜構(gòu)造中的應(yīng)用.地球物理學(xué)報,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]張汝杰,賀振華,王理.井間地震高斯射線束正演方法.物探化探計算技術(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]鄧飛,孫祥娥.三維高斯射線束正演的研究與實現(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]孫成禹,張文穎,倪長寬等.能量約束下的高斯射線束法地震波場正演.石油地球物理勘探,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é)磊.基于高斯束的彈性波場計算方法研究[學(xué)位論文].吉林長春:吉林大學(xué),2013.

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