孫輝, 岳玉波, 李猛

1 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院， 成都 611756 2 西南交通大學(xué)高速鐵路線路工程教育部重點(diǎn)實(shí)驗室， 成都 610031 3 東方地球物理公司物探技術(shù)研究中心， 涿州 072751

0 引言

地震波場數(shù)值模擬是勘探地震學(xué)重要組成部分，可以揭示地震波的傳播規(guī)律，有效指導(dǎo)地震資料的采集、處理和解釋，同時還是多種地震反演方法的基礎(chǔ)(Long et al., 2019; Wang and Zhan, 2020).傳統(tǒng)地震正演方法通常忽略介質(zhì)黏滯性對地震波傳播的影響，但這一性質(zhì)會顯著影響地震波的動力學(xué)特征，包括能量衰減、頻帶變窄、相位畸變等(Walcott, 1970; White, 1975；Zong et al., 2020).因此研究快速、高精度黏聲介質(zhì)地震正演方法具有重要意義.

地震波場數(shù)值模擬方法主要分為兩大類.一類是基于波動方程的數(shù)值解法，這類方法雖然計算效率不高但具有較高的計算精度，且可以計算全波場，因此眾多學(xué)者針對這一類方法在黏聲介質(zhì)中的應(yīng)用展開了研究(Arntsen et al., 1998; Operto et al., 2007; Abubakar and Habashy, 2013; 吳玉等，2017；Li et al., 2019; Wang et al., 2020).另一類方法是射線類地震波場數(shù)值模擬方法，這一類方法基于射線理論，具有較高的計算效率，但傳統(tǒng)射線類正演方法存在計算陰影區(qū)，在焦散區(qū)無法獲得準(zhǔn)確的振幅并且難以處理多至走時問題(Jaramillo and Bleistein, 1999; Duquet et al., 2000; Fawcett, 2001; Chapman, 2004; Shekar and Tsvankin, 2014).

雖然射線Born正演相對于傳統(tǒng)射線類地震正演方法的計算精度更高(Bleistein et al., 2001; 符立耘，2010)，但由于該方法采用射線類方法計算格林函數(shù)，仍然無法從根本上解決射線類正演方法所存在的問題.作為解決辦法，可以使用基于高斯束疊加的格林函數(shù)代入Born正演公式中.由高斯束疊加獲取的格林函數(shù)可以處理多路徑及焦散問題，其有效性已在相關(guān)偏移成像方法中得到了驗證(Hill, 1990, 2001; Gray and Bleistein, 2009; Liu and Palacharla, 2010; 曹文俊等；2013；楊繼東等，2016；韓建光等，2017；楊晶和黃建平，2018；Sun et al.,2020).Shekar和Tsvankin(2014)將基于高斯束疊加的格林函數(shù)和Kirchhoff散射理論相結(jié)合提出了一種針對各向異性衰減介質(zhì)的地震波正演方法.Huang等(2016)以及岳玉波等(2019)將由高斯束疊加獲取的格林函數(shù)融入了Born散射理論中，提出了一種針對聲波介質(zhì)的地震波正演方法，并取得了良好的數(shù)值模擬結(jié)果.

本文提出了一種針對黏聲介質(zhì)的二維高斯波束Born正演方法，在Born正演理論的框架下采用由高斯束疊加獲得的格林函數(shù)，在高斯束傳播算子中考慮了介質(zhì)黏滯性信息，并將其融入了wavelet-bank平面波合成方法中.在接下來的章節(jié)中，首先推導(dǎo)出基于高斯束的Born正演公式，接著介紹正演方法的實(shí)現(xiàn)策略，最后通過兩個數(shù)值模擬對本文提出方法的計算精度和計算效率進(jìn)行驗證.

1 基本原理

1.1 2D黏聲介質(zhì)中基于高斯束的格林函數(shù)計算

(1)

其中G表示格林函數(shù)；x為散射點(diǎn)；x′為源點(diǎn)；p=(px,pz)為慢度，px和pz分別為慢度在水平方向和垂直方向上的分量；ω為角頻率.高斯波束uGB(x,x′,p,ω)可以進(jìn)一步表示為:

uGB(x,x′,p,ω)=

(2)

其中s、n分別為射線中心坐標(biāo)系中的坐標(biāo)；p(s)、q(s)為動力學(xué)射線參數(shù)，可以通過求解動力學(xué)射線方程組獲得；τ(s)為沿著射線路徑的走時.

令：

(3)

式(1)中格林函數(shù)可以表示為：

(4)

其中：

(5)

AGB(x,x′)和τGB(x,x′)分別為復(fù)振幅和復(fù)走時，下標(biāo)R和I表示復(fù)數(shù)的實(shí)部和虛部.

黏聲介質(zhì)有多種模型，本文采用的是常Q模型，此時介質(zhì)黏滯性強(qiáng)弱由與頻率無關(guān)的品質(zhì)因子Q決定.地震波在黏聲介質(zhì)中的傳播可以看作為以一個復(fù)速度在聲波介質(zhì)中傳播，當(dāng)滿足條件Q?1,此時的復(fù)速度V(x)可以表示為(Aki and Richards，2002)：

(6)

其中ωr為參考頻率.

(7)

其中：

(8)

2D黏聲介質(zhì)中基于高斯束的Born正演公式

Born散射理論可以通過省略高階項實(shí)現(xiàn)針對一次散射波的正演，此時的Born正演公式為(Moser，2012)：

u(xr,xs,ω)=

(9)

其中xs和xr分別表示源點(diǎn)和接收道的位置；F(ω)為震源函數(shù)；x為散射點(diǎn)；D為散射點(diǎn)的集合；c0(x)為背景速度；c1(x)為擾動速度；G(x,xs,ω)和G(x,xr,ω)分別為源點(diǎn)和接收點(diǎn)的格林函數(shù)，它們的表達(dá)式為：

(10)

將式(10)代入式(9)可得黏聲介質(zhì)中基于高斯束的Born正演公式：

u(xr,xs,ω)=

(11)

-iωpLx(xr-L)]，

(12)

其中pLx和pLz為中心道初始慢度在水平和垂直方向上的參數(shù).

結(jié)合公式(13)：

(13)

式(11)可以進(jìn)一步轉(zhuǎn)化為：

×U(L,xs,psx,pLx,ω)exp[-iωpLx(xr-L)

(14)

其中ΔL為窗間隔；ωr為參考頻率；w0為初始波束寬度；U(L,xs,psx,pLx,ω)的表達(dá)式為：

U(L,xs,psx,pLx,ω)=

(15)

式(15)的含義是將反射率信息映射為局部平面波，接著通過式(14)將局部平面波轉(zhuǎn)化為窗內(nèi)接收道的地震記錄.式(15)為整個計算的核心也決定了正演算法的計算效率，下一小節(jié)將著重討論式(15)的計算方法.

1.2 基于wavelet-bank的時間域?qū)崿F(xiàn)策略

直接按照式(15)在頻率域中對U(L,xs,psx,pLx,ω)進(jìn)行計算效率并不高，文中將給出高效的時間域計算方法.將式(5)和式(7)代入到式(15)中可得：

U(L,xs,psx,pLx,ω)=

(16)

對U(L,xs,psx,pLx,ω)進(jìn)行傅里葉反變換可將其轉(zhuǎn)換到時間域：

(17)

(18)

首先針對τI和τQ建立規(guī)則的序列：

(19)

(20)

(21)

隨后通過雙線性插值以及高斯束的振幅和走時將散射點(diǎn)的反射率信息映射到Rjk(t)和Ijk(t)中：

(22)

(23)

整個正演算法的實(shí)現(xiàn)流程如圖1所示.

2 模型試算

在本節(jié)中，將使用一個水平層狀黏聲模型和一個復(fù)雜黏聲模型對所提出新正演方法的計算能力進(jìn)行驗證.

2.1 水平層狀黏聲模型計算

圖2為水平層狀黏聲模型的的速度場分布和反射率信息.模型橫向和縱向的網(wǎng)格點(diǎn)數(shù)均為301；橫向和縱向網(wǎng)格間距分別為20 m和10 m.從上到下各層的聲波速度分別為2000 m·s-1、2500 m·s-1、3000 m·s-1；Q值分別為60、50以及40.炮點(diǎn)位于模型頂部中間位置，接收道以炮點(diǎn)為中心向兩側(cè)均勻分布201道，道間距為20 m，正演使用的子波為雷克子波，主頻為20 Hz，時間采樣間隔為4 ms，時間采樣點(diǎn)數(shù)為750.

圖3為水平層狀黏聲模型采用不同方法獲得的正演結(jié)果，圖3a為聲波高斯束Born正演結(jié)果，圖3b為黏聲高斯束Born正演結(jié)果.圖4進(jìn)一步展示了兩種方法在零偏移距處的波形對比.通過對比可以發(fā)現(xiàn)：沒有考慮介質(zhì)黏滯性的聲波正演結(jié)果中與地層對應(yīng)的同相軸能量更強(qiáng)，黏聲介質(zhì)高斯束Born正演方法可以很好的將介質(zhì)黏滯性對地震波場的影響反映到正演結(jié)果中.在測試中，由于考慮介質(zhì)黏滯性的相關(guān)計算更加耗時，黏聲正演結(jié)果總共耗時2.2 s，比聲波正演多消耗了約50%的時間.

圖1 黏聲介質(zhì)高斯束Born正演實(shí)現(xiàn)流程Fig.1 Flow chart of the Born forward modeling for visco-acoustic media using Gaussian beam

2.2 各向異性復(fù)雜模型計算

圖5為復(fù)雜黏聲模型的的速度場分布和反射率分布.模型橫向網(wǎng)格點(diǎn)數(shù)為250；縱向網(wǎng)格點(diǎn)數(shù)為750；橫向和縱向網(wǎng)格間距分別為10 m和5 m.模型為常Q模型，Q的取值為50.炮點(diǎn)位于橫向1000 m處，接收道以炮點(diǎn)為中心向兩側(cè)均勻分布50道，道間距為15 m，正演使用的子波為雷克子波，主頻為15 Hz，時間采樣間隔為4 ms，時間采樣點(diǎn)數(shù)為750.

圖6為復(fù)雜黏聲模型采用不同方法獲得的正演結(jié)果，圖6a為黏聲有限差分方法正演結(jié)果，圖6b為黏聲高斯束Born正演結(jié)果.通過對比可以發(fā)現(xiàn)：有限差分方法可以模擬直達(dá)波等全波場，黏聲高斯束Born正演方法只針對一次散射波進(jìn)行模擬.兩種方法得到的一次散射波場結(jié)果幾乎相同.就計算時間而言，黏聲有限差分正演方法共耗時101.4 s，黏聲高斯束Born正演方法共耗時1.3 s，可以看出本文提出的正演方法在計算效率上具有很大的優(yōu)勢.

3 結(jié)論

本文提出了一種針對二維黏聲介質(zhì)一次散射波場模擬的高斯束Born正演方法.該方法通過疊加不同初射方向的高斯束獲取格林函數(shù)，克服了傳統(tǒng)射線類方法計算格林函數(shù)的缺點(diǎn)，保證了格林函數(shù)以及波場模擬的精度.同時基于黏聲高斯束所包含的信息，新的正演方法采用wavelet-bank方法快速有效的將反射率信息映射成為了窗內(nèi)的局部平面波，保證了正演算法的計算效率.通過兩個數(shù)值模型測試表明，本文提出的正演方法具有良好的計算精度，并且計算效率相較于黏聲介質(zhì)有限差分地震波正演要高很多.

圖2 水平層狀黏聲模型(a) 聲波速度分布； (b) 反射率分布.Fig.2 Layered visco-acoustic model(a) Acoustic velocity distribution； (b) Reflectivity distribution.

圖3 水平層狀黏聲模型不同方法正演結(jié)果(a) 聲波高斯束Born正演結(jié)果； (b) 黏聲高斯束Born正演結(jié)果.Fig.3 Modeling results with different modeling methods for layered visco-acoustic model(a) Result of acoustic Gaussian beam Born modeling； (b) Result of visco-acoustic Gaussian beam Born modeling.

圖4 零偏移距不同正演方法波形對比虛線為圖3a中零偏移距波形，實(shí)線為圖3b中零偏移距波形.Fig.4 Waveform comparison of zero-offset traces between different modeling methodsDashed line marks the waveform in Fig.3a, solid line is the waveform in Fig.3b.

圖5 復(fù)雜黏聲模型(a) 聲波速度分布； (b) 反射率分布.Fig.5 Complex visco-acoustic model(a) Acoustic velocity distribution； (b) Reflectivity distribution.

圖6 復(fù)雜黏聲模型不同方法正演結(jié)果(a) 黏聲有限差分正演結(jié)果； (b) 黏聲高斯束Born正演結(jié)果.Fig.6 Forward modeling results with different methods for layered visco-acoustic model(a) Result of visco-acoustic finite-difference modeling；(b) Result of visco-acoustic Gaussian beam Born modeling.