基于全局笛卡爾坐標系的高斯束地震波場模擬

白　敏　吳　娟*　孫章慶　楊　磊

(①華北水利水電大學資源與環境學院，河南鄭州 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.