韓冰凱　顧漢明*　劉少勇

(①中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院，湖北武漢 430074； ②地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074)

1　引言

經(jīng)典射線(束)傳播算子在地震波偏移成像中應(yīng)用廣泛，其高頻近似假設(shè)的成立條件是介質(zhì)在第一菲涅耳帶內(nèi)足夠光滑[1]，導(dǎo)致射線類傳播算子在復(fù)雜介質(zhì)中的應(yīng)用受到限制?；陔p程波波動(dòng)方程的逆時(shí)偏移，能夠準(zhǔn)確描述復(fù)雜介質(zhì)中帶限波場的傳播，但計(jì)算效率偏低，且輸出角度域共成像點(diǎn)道集需要額外計(jì)算量[2,3]。在射線理論框架下，發(fā)展更準(zhǔn)確地傳播帶限地震波場的傳播算子進(jìn)行偏移成像和速度分析，有重要的理論意義和應(yīng)用潛力。

Foreman[4]推導(dǎo)了與頻率相關(guān)的射線追蹤系統(tǒng)，相比于傳統(tǒng)的射線追蹤系統(tǒng)，其計(jì)算復(fù)雜度顯著增加。Protasov等[5]分析了與頻率相關(guān)的射線追蹤系統(tǒng)的算法結(jié)構(gòu)，并與傳統(tǒng)射線追蹤、時(shí)間域有限差分求解波動(dòng)方程的波場進(jìn)行了對比，結(jié)果表明對應(yīng)震源子波頻帶特征的帶限射線能夠改善射線陰影區(qū)成像。Biondi[6]從波動(dòng)方程出發(fā)，推導(dǎo)出與頻率相關(guān)的程函方程，通過構(gòu)建非線性偏微分方程求解與頻率相關(guān)的程函方程，但其低頻部分的求解過程需利用已計(jì)算的高頻部分的結(jié)果，需額外的計(jì)算量。Lomax[7]通過對垂直于射線平面波長范圍內(nèi)的速度場平滑，以加權(quán)平均的速度進(jìn)行射線追蹤。Zelt等[8]在有限差分求解程函方程過程中，對差分節(jié)點(diǎn)的速度用波長范圍內(nèi)的加權(quán)平均速度代替，并指出該平滑過程等價(jià)于波長范圍內(nèi)的速度預(yù)平滑，在不改變傳統(tǒng)有限差分求解程函方程的前提下高效計(jì)算有限頻旅行時(shí)場，并發(fā)展了基于該速度預(yù)平滑算子的有限頻層析反演算法。Chen等[9]通過與有限差分波動(dòng)方程正演模擬的波場和波形記錄對比，分析了波長內(nèi)速度平滑算子對有限頻旅行時(shí)場的影響。Protasov等[10]提出一種與頻率有關(guān)的射線追蹤算法，指出帶限射線通過速度界面時(shí)，其方向?qū)?yīng)第一菲涅耳帶內(nèi)地震波場傳播的最大能量方向，即帶限射線的傳播遵循帶限Snell定律。Yarman等[11]從Kirchhoff邊界積分出發(fā)，推導(dǎo)了帶限Snell定律的表達(dá)式，并將帶限射線追蹤應(yīng)用于Kirchhoff積分偏移中。Yarman等的算法不需要對速度模型進(jìn)行平滑，但需要已知明確的速度界面信息，這在一定程度上限制了此方法的應(yīng)用。劉少勇等[12]基于帶限射線追蹤構(gòu)建了各向同性介質(zhì)的帶限局部平面波傳播算子，并應(yīng)用于射線束偏移。

地下介質(zhì)的各向異性普遍存在，且介質(zhì)各向異性對成像聚焦和成像點(diǎn)位置的影響不容忽略[13]。隨著油氣勘探技術(shù)的發(fā)展，在地震資料處理過程中考慮各向異性效應(yīng)已經(jīng)逐漸從前沿技術(shù)探索轉(zhuǎn)變成為常規(guī)處理流程?；诓煌母飨虍愋缘卣鸩▊鞑ニ阕?，可以發(fā)展出不同的各向異性介質(zhì)成像算法?；陔p程波方程的逆時(shí)偏移技術(shù)可以擴(kuò)展到適應(yīng)橫各向同性(TI)介質(zhì)[14,15]，但是其計(jì)算效率低。通過引入TI介質(zhì)射線追蹤[16,17]或適用于TI介質(zhì)的波前類旅行時(shí)計(jì)算方法，如動(dòng)態(tài)規(guī)劃法[18]和有限差分求解程函方程法[19,20]，Kirchhoff積分疊前深度偏移可以方便地推廣到TI介質(zhì)[21,22]。經(jīng)典的高斯束成像[23-26]也可以方便地?cái)U(kuò)展到TI介質(zhì)中[27-29]。各向異性介質(zhì)中射線類偏移算法同樣也受限于高頻近似假設(shè)，存在射線照明的陰影區(qū)和焦散現(xiàn)象等[30]。本文基于TI介質(zhì)射線追蹤系統(tǒng)，在局部平面波近似下構(gòu)建TI介質(zhì)帶限射線追蹤算法；結(jié)合旁軸近似，基于帶限中心射線構(gòu)建帶限射線束傳播算子，并應(yīng)用于TI介質(zhì)的射線束偏移。

2　方法原理

2.1　帶限射線追蹤

基于高頻近似的射線追蹤是假設(shè)地震波的能量沿寬度為零的射線傳播，高頻射線也稱作費(fèi)馬射線。帶限射線理論描述地震波能量沿有寬度的波路徑傳播，為一類胖射線[31]。如圖1所示，對于波長尺度的速度異常體(圖1a)，波傳播路徑(紅色箭頭)受到影響，而經(jīng)過異常體附近的高頻射線路徑(黑色箭頭)并沒有改變；對于遠(yuǎn)小于一個(gè)波長的速度異常體(圖1b)，波傳播路徑?jīng)]有發(fā)生明顯改變，但高頻射線路徑發(fā)生明顯偏轉(zhuǎn)[7]。因此，高頻射線難以準(zhǔn)確描述帶限波場傳播的特征。

從聲波波動(dòng)方程出發(fā)，含震源介質(zhì)透過界面Γ的透射波場[32]表達(dá)為

uI(s,x,f)?nG2(x,r,f)]T(x)e-2πftdsΓ(x)

(1)

圖1　高頻射線與帶限波路徑對比示意圖[7]

uT(s,r,t)≈4?Γ{A(x,s,r)T(x)cos(2πfc[τ(x,s,r)-

t] sinc(2πfb[τ(x,s,r)-t])}dsΓ(x)

(2)

(3)

對于帶寬為B的帶限震源，在x0附近，限定第一菲涅耳帶內(nèi)的Kirchhoff積分由u(s,r,x0)=uT[s,r,τ(x0,s,r)]給出[11]，即

sinc(2πB[Δτ(x,x0,s,r)])dsΓ(x)

(4)

式中： FZx0表示x0附近由式(3)限定的第一菲涅耳帶； Δτ(x,x0,s,r)=τ(x,s,r)-τ(x0,s,r)。帶限Snell定律定義為使式(4)中透射波場u(s,r,x0)振幅取得最大值的波傳播方向，即x0處第一菲涅耳帶FZx0內(nèi)射線入射方向p1(x)和出射方向p2(x)滿足如下關(guān)系式

(5)

式中 ‖p2‖2=1。在局部平面波的假設(shè)下，若A(x,s,r)和T(x)為常數(shù)，對于界面上的任意一點(diǎn)x∈FZx0，當(dāng)p2(x0)滿足下式時(shí)，式(4)有最大值

(6)

式中v1和v2分別表示界面兩側(cè)的速度。式(6)的物理含義為，在第一菲涅耳帶FZx0內(nèi)，存在等效的局部平面波，能夠描述式(4)表征的穿過界面Γ的帶限透射波場，透射波場的傳播方向p2對應(yīng)頻帶B的帶限波場傳播的最大能量方向。由式(5)表征的帶限Snell定律控制帶限射線追蹤系統(tǒng)，帶限射線傳播可由式(6)表達(dá)的等效局部平面波實(shí)現(xiàn)。

2.2　TI介質(zhì)帶限射線束傳播與偏移

在各向異性介質(zhì)中，以Thomsen各向異性參數(shù)表征的qP波相速度是相角的函數(shù)[33]。Alkhalifa[17]通過qS波垂直各向同性面的相速度為零的假設(shè)，得到TI介質(zhì)的各向異性相速度表達(dá)式。Jiang等[34]分析了qP波相速度對彈性剛度系數(shù)C13的敏感性，在|δ|?1的假設(shè)下，對C13以線性近似推導(dǎo)出VTI介質(zhì)近似qP波相速度公式

(7)

(8)

上述射線控制方程可用于描述TI介質(zhì)中局部平面波的傳播。如圖2a所示，在x0處第一菲涅耳帶FZx0內(nèi)構(gòu)建局部平面波，以灰色矩形表示的局部平面波由一定密度且方向相同的各向異性高頻射線構(gòu)成，對應(yīng)射線參數(shù)(pin,pout)的單條射線，在x∈Γx0的傳播遵循式(8)。進(jìn)一步，由透射射線參數(shù)pout加權(quán)得到等效射線參數(shù)p2。數(shù)學(xué)上，該過程對應(yīng)懲罰區(qū)間為第一菲涅耳帶的最優(yōu)化問題，可在最小二乘框架下求解。由此，根據(jù)帶限Snell定律求解各向異性帶限射線的傳播方向，對應(yīng)地描述各向異性介質(zhì)中帶限波場的傳播方向。如圖2b所示，由xs發(fā)出的射線束，經(jīng)過速度界面Γ后，對于不同方向的格林函數(shù)的波場能量不同。基于帶限中心射線的射線束(紅色)對應(yīng)最大能量方向，描述帶限波場的傳播方向。

下面以該帶限射線為中心射線，結(jié)合旁軸射線束構(gòu)建各向異性帶限射線束傳播算子。將中心射線旅行時(shí)在射線束寬度內(nèi)展開

(9)

式中：x0為中心射線坐標(biāo)；xΩ為射線束寬度內(nèi)的坐標(biāo)； Δs為xΩ在中心射線上的投影與x0間的距離；r為xΩ到中心射線的距離；mr為與射線束波前曲率相關(guān)的動(dòng)力學(xué)參數(shù)

(10)

圖2　帶限射線追蹤(a)及傳播示意圖(b)

地震波偏移成像可表達(dá)為地震波傳播算子加上成像條件，基于各向異性帶限射線束傳播算子可構(gòu)建各向異性帶限射線束偏移算法。將傳統(tǒng)的Kirchhoff積分偏移成像條件應(yīng)用于共炮集數(shù)據(jù)的射線束偏移方程中

pr(xr0)]dΩ|τ=τs+τr

(11)

式中：x為地下成像點(diǎn)；Ω包含炮點(diǎn)xs和檢波點(diǎn)中心點(diǎn)xr0的成像孔徑；ps,pr為炮點(diǎn)和檢波點(diǎn)中心的初始射線參數(shù)；w為成像點(diǎn)的加權(quán)函數(shù)[36]；d為時(shí)間域表征的局部平面波；τ為旅行時(shí)。該共炮集成像方程可在時(shí)間空間域高效實(shí)現(xiàn)。

3　數(shù)值實(shí)驗(yàn)

為驗(yàn)證TI介質(zhì)帶限射線束傳播算子的有效性和正確性，本文選取四類簡單模型和各向異性SEG/Hess模型進(jìn)行數(shù)值實(shí)驗(yàn)。

首先，選取四類簡單各向異性模型進(jìn)行帶限射線照明分析，圖3為四類簡單模型高頻射線與帶限射線照明分析，模型網(wǎng)格為10m×10m。對于同一模型設(shè)置相同震源坐標(biāo)，分別進(jìn)行TI介質(zhì)的射線追蹤、帶限射線追蹤和有限差分求解波動(dòng)方程正演模擬，對比相同時(shí)刻的射線旅行時(shí)波前與有限差分計(jì)算波場快照。對比高頻射線和帶限射線的射線路徑與波前面可知，帶限射線較高頻射線在模型中分布更均勻，提高了傳統(tǒng)射線陰影區(qū)的射線覆蓋率，能夠提供更有效的照明。

然后，選擇各向異性SEG/Hess模型對TI介質(zhì)帶限射線束傳播算子及其對應(yīng)的偏移成像方法進(jìn)行照明分析和偏移脈沖響應(yīng)測試。如圖4所示，SEG/Hess模型中存在高速鹽丘體和三個(gè)各向異性透鏡體，同時(shí)模型包含陡傾斷裂構(gòu)造，模型網(wǎng)格為5m×5m。在x=8.0km處設(shè)置震源，分別進(jìn)行射線追蹤得到高頻射線路徑和帶限射線路徑(圖5a和圖5b)。對比高頻射線和帶限射線路徑可知，帶限射線在鹽丘下方區(qū)域“①”、鹽丘側(cè)方各向異性透鏡體區(qū)域“②”、模型斷層區(qū)域“③”有更均勻的射線分布。對應(yīng)地，以帶限射線為中心射線構(gòu)建的帶限射線束進(jìn)一步提高了鹽下射線陰影區(qū)和各向異性透鏡體的照明。相比于傳統(tǒng)射線束偏移脈沖響應(yīng)(圖5c)，帶限射線束偏移脈沖響應(yīng)(圖5d)的能量連續(xù)性更好。

圖6a和6b展示了TI介質(zhì)傳統(tǒng)射線束偏移和TI介質(zhì)帶限射線束偏移結(jié)果，圖6c和6d為圖6a及圖6b中矩形框內(nèi)區(qū)域的放大展示。 對比箭頭所示的區(qū)域，相比于TI介質(zhì)傳統(tǒng)射線束偏移，TI介質(zhì)帶限射線束偏移結(jié)果中對于復(fù)雜鹽丘邊界和陡傾斷裂構(gòu)造成像效果更聚焦，同相軸能量分布更均勻，對于各向異性透鏡體(A、B和C)及其相鄰鹽丘側(cè)邊界的刻畫更準(zhǔn)確。圖7a和7b分別為兩種射線束方法5處的角度域成像道集(x=2.75， 6.95， 8.40，9.50，12.40km)，對比圖中紅色箭頭標(biāo)注的區(qū)域，TI介質(zhì)帶限局部平面波傳播算子能夠?yàn)閷?yīng)的射線束偏移提供更有效的角度域照明，TI介質(zhì)帶限射線束偏移可輸出質(zhì)量更高的角度域共成像點(diǎn)道集。

圖3　四類簡單模型高頻射線與帶限模型照明分析

(a)向斜各向異性模型； (b)圖a的高頻射線路徑及波前面； (c)圖a的帶限射線路徑及波前面； (d)背斜各向異性模型； (e)圖d的高頻射線路徑及波前面； (f)圖d的帶限射線路徑及波前面； (g)崎嶇震蕩各向異性模型； (h)圖g的高頻射線路徑及波前面； (i)圖g的帶限射線路徑及波前面； (j)楔形各向異性模型； (k)圖j的高頻射線路徑及波前面； (l)圖j的帶限射線路徑及波前面

各類各向異性模型中兩種介質(zhì)參數(shù)：vP0為2.5km/s和4.5km/s，ε為0.15和0.05，δ為0.06和0.02。模型速度界面以藍(lán)色實(shí)線表示； 射線以紅色實(shí)線表示； 射線范圍-65°～65°，間隔2.5°； 0.8s旅行時(shí)波前面以綠色實(shí)心圓點(diǎn)表示；背景為有限差分計(jì)算的震源主頻20Hz的0.8s波場快照

圖4　SEG/Hess模型各向異性參數(shù) (a)vP0； (b)ε； (c)δ

4　結(jié)論

在射線理論框架下，基于局部平面波假設(shè)，提出了適用于各向異性介質(zhì)的帶限射線束傳播算子，分析了該算子對鹽丘等復(fù)雜構(gòu)造和各向異性結(jié)構(gòu)的照明，并發(fā)展了基于該算子的帶限射線束偏移方法。數(shù)值實(shí)驗(yàn)表明，各向異性帶限射線束能夠改善鹽下、斷裂等復(fù)雜各向異性構(gòu)造的角度域照明，提高偏移剖面和角度域共成像點(diǎn)道集的質(zhì)量。本文發(fā)展的TI介質(zhì)帶限波傳播算子和偏移成像方法為后續(xù)基于該算子的層析反演、偏移速度分析等奠定了基礎(chǔ)。

圖5　TI介質(zhì)射線束照明與偏移脈沖響應(yīng)對比

(a)TI介質(zhì)高頻射線路徑及旅行時(shí)波前面； (b)TI介質(zhì)帶限射線路徑及旅行時(shí)波前面； (c)TI介質(zhì)傳統(tǒng)射線束偏移脈沖響應(yīng)； (d)TI介質(zhì)帶限射線束偏移脈沖響應(yīng)

炮點(diǎn)坐標(biāo)(8km,0)，射線路徑以紅色實(shí)線表示； 初始射線角度范圍-60°～60°，角度間隔2.5°； 2.4s旅行時(shí)波前面以綠色實(shí)心圓點(diǎn)表示；背景為有限差分正演模擬的震源主頻為20Hz的2.4s波場快照。偏移脈沖響應(yīng)為記錄于x=8km，記錄時(shí)間為4.353s的20Hz雷克子波，對應(yīng)于坐標(biāo)(8.0km，5.825km)的反射點(diǎn)感謝勘探地球物理學(xué)家協(xié)會(huì)(SEG)和Amerada Hess石油公司提供公開模型數(shù)據(jù)。

圖6　TI介質(zhì)射線束偏移結(jié)果對比

[1]Kravtsov Y A,Orlov Y I.Geometrical Optics of Inhomogeneous Media.Springer-Verlag,1990.

[2]Yoon K,Marfurt K J.Reverse-time migration using the Poynting vector.Exploration Geophysics,2006,37(1)：102-107.

[3]Hu J,McMechan G A,Guan H.Comparison of methods for extracting ADCIGs from RTM.Geophysics,2014,79(3)：S89-S103.

[4]Foreman T.A Frequency Dependent Ray Theory[D].The University of Texas at Austin,1987.

[5]Protasov M,Gadylshin K.Exact frequency dependent rays on the basis of Helmholtz solver.SEG Technical Program Expanded Abstracts,2015,34：3739-3743.

[6]Biondi B.Solving the frequency dependent eikonal equation.SEG Technical Program Expanded Abstracts,1992,11：1315-1319.

[7]Lomax A.The wavelength-smoothing method for approximating broad-band wave propagation through complicated velocity structures.Geophysical Journal International,1994,117(2)：313-334.

[8]Zelt C A,Chen J.Frequency-dependent traveltime tomography for near-surface seismic refraction data.Geophysical Journal International,2016,207(1)：72-88.

[9]Chen J,Zelt C A.Comparison of full wavefield syn-thetics with frequency-dependent traveltimes calculated using wavelength-dependent velocity smoothing.Journal of Environmental and Engineering Geophy-sics,2017,22(2)：133-141.

[10]Protasov M I,Yarman C E,Nichols D et al.Frequency-dependent ray-tracing through rugose interfaces.SEG Technical Program Expanded Abstracts,2011,30：2992-2996.

[11]Yarman C E,Cheng X,Osypov K et al.Band-limited ray tracing.Geophysical Prospecting,2013,61(6)：1194-1205.

[12]劉少勇,韓冰凱,顧漢明等.帶限射線束傳播算子及偏移方法.石油地球物理勘探,2017,52(5)：948-955.

Liu Shaoyong,Han Bingkai,Gu Hanming et al.Study on band-limited beam propagator and migration method.OGP,2017,52(5)：948-955.

[13]Isaac J H,Lawton D C.Image mispositioning due to dipping TI media:A physical seismic modeling study.Geophysics,1999,64(4)：1230-1238.

[14]Fowler P J,Du X,Fletcher R P.Coupled equations for reverse time migration in transversely isotropic media.Geophysics,2010,75(1)：S11-S22.

[15]Delaney S J,O’Brien G S,Short R et al.Tilted transverse isotropic reverse time migration with angle gathers:Implementation and efficiency.Geophysics,2016,81(6)：S419-S432.

[17]Alkhalifah T.An acoustic wave equation for aniso-tropic media.Geophysics,2000,65(4)：1239-1250.

[18]Kumar D,Sen M K,Ferguson R J.Traveltime calculation and prestack depth migration in tilted transversely isotropic media.Geophysics,2004,69(1)：37-44.

[19]Qian J,Symes W.Finite-difference quasi-P traveltime for anisotropic media.Geophysics,2002,67(1)：147-155.

[20]Sun J,Sun Z,Han F.A finite difference scheme for solving the eikonal equation including surface topography.Geophysics,2011,76(4)：T53-T63.

[21]Vanelle C,Gajewski D.True-amplitude Kirchhoff depth migration in anisotropic media:The traveltime-based approach.Geophysics,2013,78(5)：WC33-WC39.

[21]Liu S,Wang H,Yang Q et al.Traveltime computation and imaging from rugged topography in 3D TTI media.Journal of Geophysics and Engineering,2014,11(1)：1-9.

[23]Hill N R.Gaussian beam migration.Geophysics,1990,55(11)：1416-1428.

[24]Hill N R.Prestack Gaussian-beam depth migration.Geophysics,2001,66(4)：1240-1250.

[25]鄧飛,劉超穎,趙波等.高斯射線束正演與偏移.石油地球物理勘探,2009,44(3)：265-269.

Deng Fei,Liu Chaoying,Zhao Bo et al.Gaussian beams forward simulation and migration.OGP,2009,44(3)：265-269.

[26]李振春,岳玉波,郭朝斌等.高斯波束共角度保幅深度偏移.石油地球物理勘探,2010,43(3):360-365.

Li Zhenchun,Yue Yubo,Guo Chaobin et al.Gaussian beam common angle preserved-amplitude migration.OGP,2010,45(3)：360-365.

[27]Zhu T,Gray S H,Wang D.Prestack Gaussian-beam depth migration in anisotropic media.Geophysics,2007,72(3)：S133-S138.

[28]韓建光,王赟,張曉波等.VTI介質(zhì)高斯束疊前深度偏移.石油地球物理勘探,2015,50(2)：267-273.

Han Jianguang, Wang Yun, Zhang Xiaobo et al.Gaussian beam prestack depth migration in VTI media.OGP,2015,50(2)：267-273.

[29]劉強(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,51(5)：930-937.

[31]Woodward M J.Wave-equation tomography.Geophysics,1992,57(1)：15-26.

[32]Chapman C H.Fundamentals of Seismic Wave Propagation.Cambridge University Press,2004.

[33]Thomsen L.Weak elastic anisotropy.Geophysics,1986,51(10)：1954-1966.

[34]Jiang B,Wang H Z,Liu S Y.A phase-velocity based ray tracing system for anisotropic media.76th EAGE Conference & Exhibition,Amsterdam,2014,We G105 15.

[35]Liu S,Wang H,Feng B.The characteristic wave decomposition and imaging in VTI media.Journal of Applied Geophysics,2015,115(13)：51-58.

[36]Hu C,Stoffa P L.Slowness-driven gaussian-beam prestack depth migration for low-fold seismic data.Geophysics,2009,74(6)：WCA35-WCA45.