摘要：基于Katchalov提出的矢量時(shí)空域高斯束表達(dá)式，構(gòu)建彈性波介質(zhì)中的矢量時(shí)空域格林函數(shù)；使用Kirchhoff-Helmholtz積分解得到檢波點(diǎn)矢量波場，并利用散度、旋度算子在矢量波場中分離出純縱波和純橫波分量；基于衰減補(bǔ)償原理，對波場進(jìn)行衰減補(bǔ)償，最后使用震源歸一化成像條件得到PP、PS成像結(jié)果。數(shù)值實(shí)驗(yàn)證明了所提方法的正確性和適應(yīng)性。結(jié)果表明，所提方法相比于傳統(tǒng)的高斯束方法，成像精度近似，計(jì)算效率方面有較大優(yōu)勢。

關(guān)鍵詞：高斯束； 時(shí)空域高斯波包； 多分量； 黏彈介質(zhì)

中圖分類號：P 631.4"" 文獻(xiàn)標(biāo)志碼：A

引用格式：陳超，李振春，黃建平. 黏彈介質(zhì)時(shí)空域高斯波包偏移［J］. 中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2024，48（5）：46-58.

CHEN Chao， LI Zhenchun， HUANG Jianping. Time-space" Gaussian packets migration in viscoelastic media ［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2024，48（5）：46-58.

Time-space Gaussian packets migration in viscoelastic media

CHEN Chao1， LI Zhenchun1， HUANG Jianping1，2

（1.School of Geosciences in China University of Petroleum（East China）， Qingdao 266580， China;

2.State Key Laboratory of Deep Oil and Gas， China University of Petroleum（East China）， Qingdao 266580， China）

Abstract：Based on the vector time-space domain Gaussian beam expression proposed by Katchalov， the vector time-space domain Greens function in elastic wave media is constructed. The Kirchhoff-Helmholtz integral is applied to obtain the vector wave field of the receiver point， and the divergence and curl operators are used to separate the pure P-wave and pure S-wave components from the vector wave field. Based on the attenuation compensation principle， the wave field is adjusted and the final PP and PS imaging results are achieved using source-normalized imaging conditions. Numerical experiments confirm the correctness and adaptability of the proposed method. Compared with the traditional Gaussian beam method， the approach maintains similar imaging accuracy while offering significant improvements in computational efficiency.

Keywords： Gaussian beam; time-space Gaussian packets; multi-component; viscoelastic media

傳統(tǒng)的高斯束偏移方法在頻率域中計(jì)算，需要頻繁的傅里葉變換，而且會(huì)出現(xiàn)假頻和邊緣截?cái)鄦栴}，降低了計(jì)算效率和成像精度。而由Katchalov提出的時(shí)空域高斯波包偏移方法具有很多優(yōu)點(diǎn)：在時(shí)間域中計(jì)算，避免頻繁的傅里葉變換，減少假頻和邊緣截?cái)鄦栴}的出現(xiàn)，提高了計(jì)算效率和成像精度。因此時(shí)空域高斯波包偏移方法具有很高的研究價(jià)值。高斯束偏移方法兼具了Kirchhoff偏移方法的計(jì)算效率和波動(dòng)方程偏移的成像精度，同時(shí)具有高效、靈活的特點(diǎn)，對地震波傳播成像有較強(qiáng)的適應(yīng)性［1］。Hill［2］首先提出了高斯束偏移，介紹了高斯束偏移的基本原理，研究了鞍點(diǎn)積分提高計(jì)算效率；Gray［3］提出了真振幅成像條件，對高斯束偏移成像精度有較高的提升；李振春等［4-5］從共角度域和共炮域研究了保幅高斯束偏移方法，提高了成像質(zhì)量；黃建平等［6-7］研究了起伏地表彈性波高斯束偏移成像和保幅高斯束偏移成像，證明高斯束方法對起伏地表有較好的適應(yīng)性。目前，時(shí)空域高斯束偏移有高斯束逆時(shí)偏移和高斯波包法偏移［8］。時(shí)空域高斯波包法在沿中心射線路徑傳播時(shí)能夠描述局部波場特征，對目標(biāo)成像具有一定優(yōu)勢［9］。Katchalov等［10］提出了一種基于時(shí)空高斯束的理論地震記錄計(jì)算新方法，該方法可以克服計(jì)算焦散線附近波場的困難，并且不包含任何傅里葉變換的計(jì)算。Geng等［11-12］提出了一種基于Gabor分解和高斯波包漸近傳播的高斯波包偏移方法，與更精確但更耗時(shí)的基于單向波動(dòng)方程的偏移（如波束波偏移）相比，高斯波包法能夠以更高的效率正確成像復(fù)雜模型的主要結(jié)構(gòu)，特別是在鹽下區(qū)域。李輝等［13-15］基于高斯波包按中心射線出射角度積分存在解析式與Gabor點(diǎn)源函數(shù)產(chǎn)生波場解析式是等價(jià)的條件，提出一種利用高斯波包疊加計(jì)算任意點(diǎn)源函數(shù)產(chǎn)生的波場的新方法，并將其用于疊前深度偏移，證明了高斯波包方法是有效且快速的。之后，L等［16-17］將動(dòng)態(tài)參數(shù)控制、主頻近似的思想與時(shí)空域高斯束結(jié)合，使地震深度成像在淺層和深層方面具有更高的精度，將計(jì)算效率提高到一個(gè)甚至兩個(gè)數(shù)量級。黏滯性主要是由巖石孔隙造成的。張巖等［18］將低場核磁共振成像技術(shù)應(yīng)用到巖石物理研究中，獲取了巖心的局部孔隙結(jié)構(gòu)，實(shí)現(xiàn)了三維空間孔隙的系統(tǒng)性研究。劉學(xué)峰等［19］利用機(jī)器學(xué)習(xí)劃分孔隙和巖石礦物骨架，實(shí)現(xiàn)了對頁巖孔隙精確劃分。巖石孔隙與黏滯性密切相關(guān)，加強(qiáng)對巖石孔隙的研究對黏滯性偏移成像有重要意義。Bickel等［20］將地震波的傳播效應(yīng)使用衰減和頻散的常數(shù)Q（與頻率無關(guān)）模型進(jìn)行近似，并提出了一種時(shí)變逆濾波器消除了地震帶內(nèi)的相位畸變，減少了地震波之間的干擾，是黏滯性偏移成像的基礎(chǔ)。Zhang等［21-22］在非均勻介質(zhì)中進(jìn)行波場外推，發(fā)展了黏滯聲波效應(yīng)補(bǔ)償?shù)哪鏁r(shí)偏移方法。筆者基于時(shí)空域矢量高斯束，構(gòu)建時(shí)空域矢量格林函數(shù)，通過Kirchhoff-Helmholtz積分解得到檢波點(diǎn)矢量波場，并利用散度、旋度算子在矢量波場中分離出純縱波和純橫波分量，基于衰減補(bǔ)償原理，對縱橫波波場進(jìn)行衰減補(bǔ)償，將時(shí)間域高斯波包偏移擴(kuò)展到黏彈介質(zhì)，最后使用震源歸一化成像條件得到PP、PS成像結(jié)果。

1 方法原理

1.1 彈性介質(zhì)時(shí)空域高斯束

根據(jù)Katchalov的研究［10］，彈性介質(zhì)中矢量時(shí)空域高斯束表達(dá)式為

P（t-τ，s，n）=A0（t-τ）ρ0（s）α0（s）Q（s）×expip

θ0（t-τ）-12θ′0（t-τ）P（s）Q（s）n2eP（s），（1）

S（t-τ，s，n）=A0（t-τ）ρ0（s）β0（s）Q（s）×expipθ0（t-τ）-12θ′0（t-τ）P（s）Q（s）n2

eS（s）.（2）

式中，P（t-τ，s，n）和S（t-τ，s，n）分別為純縱波震源和純橫波震源激發(fā)產(chǎn)生的時(shí)空域高斯束；α0（s）和β0（s）為出射點(diǎn)的初始縱波速度和橫波速度；A（t）、θ（t）均為充分光滑的函數(shù)且θ′（t）=dθdt≠0，θ′（t）lt;constlt;0，A0（t-τ）為射線振幅的初始系數(shù)，θ0（t-τ）滿足中心射線的程函方程；（s，n）為射線中心坐標(biāo)系；P（s）和Q（s）為動(dòng)力學(xué)射線參量；ρ0（s）為密度；t為時(shí)間；τ為射線走時(shí)；p為數(shù)值較大的參數(shù)，當(dāng)p→∞，瞬時(shí)頻率ω=-pθ′（t）滿足高頻假設(shè)；eP（s）和eP（s）為極化矢量。eP（s）和eP（s）表達(dá)式為

eP=t+nα0（s）P（s）Q（s）n，

eS=n+nβ0（s）P（s）Q（s）t.（3）

式中，對P波，t為主分量，n為次分量；對S波，n為主分量，t為次分量。根據(jù)Hill給定的初始值［2］

Q1（s0）=1，P1（s0）=0，

Q2（s0）=0，P2（s0）=1.（4）

將式（3）代入式（1）、（2）得

P（t-τ，s，n）=exp［-ω2（t-τ）2］ρ0（s）α0（s）Q（s）×

expiω-（t-τ）+12P（s）Q（s）n2t+nα0（s）P（s）Q（s）n， （5）

S（t-τ，s，n）=exp［-ω2（t-τ）2］ρ0（s）β0（s）Q（s）×

expiω-（t-τ）+12P（s）Q（s）n2n+nβ0（s）P（s）Q（s）t .（6）

時(shí)空域位移矢量可以通過一系列震源出射、不同出射角度的時(shí)空域高斯束積分疊加表示

UP（t-τ，s，n）=∫2π0ΨP（ψ）P（t-τ，s，n）dψ，（7）

US（t-τ，s，n）=∫2π0ΨS（ψ）S（t-τ，s，n）dψ.（8）

式中，ΨP和ΨS為權(quán)因子；ψ為不同角度出射的高斯束的出射角。用式（7）和（8）求得的高頻漸進(jìn)解與均勻介質(zhì)中的位移矢量的解析解對比可求得

ΨP（ψ）=-iρ（s0）α（s0）ρ（s0）ε4π ，（9）

ΨS（ψ）=-iρ（s0）β（s0）ρ（s0）ε4π .（10）

式中，ε為射線參數(shù)，Hill［2］給定了其初始值ε=-iωrl20，ωr為參考頻率；l0為初始束寬。

1.2 點(diǎn)源函數(shù)的Gabor分解及波場構(gòu)建

根據(jù)李輝的研究［13］，點(diǎn)源函數(shù)f（t）可以分解為一系列的Gabor基函數(shù)gk（t-τk） ，即

f（t）=∑nk=1akgk（t-τk）. （11）

式中，ak為疊加權(quán)重系數(shù)；τk為第k條Gabor基函數(shù)的延遲時(shí)間，其可以通過將式（11）兩側(cè)同時(shí)Gabor基函數(shù)gi（t-τk）做內(nèi)積求得

〈gk（t-τk），f（t）〉=∑nk=1ak〈gk（t-τk），gk（t-τk）〉，k=1，2，3，….（12）

通過Gabor基函數(shù)構(gòu)建的震源函數(shù)激發(fā)不同時(shí)間延遲的P波時(shí)空域高斯束，并將其疊加可得到時(shí)空域高斯束波場，結(jié)合式（5）～（10）得到震源波場

uD

（x，xs，t）=Re-iρ（s0）4π∑kakA0（t-τ-τk）∫2π0dψα（s0）ρ（s0）εα（s）ρ（s）Q（s）×expiω-（t-τ-τk）+12P（s）Q（s）n2.（13）

式中，xs為震源位置；x為成像點(diǎn)位置。

假設(shè)S為自由地表，使用Kirchhoff-Helmholtz積分表示反向延拓的彈性波位移場uU

（x，xr，t）的表達(dá)式［23］為

uU（x，xr，t）=-∫t2t1dt∫Sdxr·u（xr，t0）·m·∑T（x，xr）.（14）

式中，xr為接收點(diǎn)位置;m為邊界S的外法向向量;u（xr，t0）為t0時(shí)刻地表檢波點(diǎn)xr處產(chǎn)生的多分量地震波場;∑（x，xr）為格林應(yīng)力張量，其表達(dá)式為

∑lmn=λδlmGkn，k+μ（Gmn，l+Gln，m）.（15）

式中，xi（i=l，m，n，k）為笛卡爾坐標(biāo)分量;δ為Kronecker Delta函數(shù);λ和μ為拉梅系數(shù);Gmn，l=Gxm，xnxl，表示彈性波格林函數(shù)的偏導(dǎo)數(shù)，其高頻近似表達(dá)式為

Gmn，liωplGmn.（16）

其中pl為慢度矢量的xl分量。

將式（7）和（8）中射線中心坐標(biāo)系的UP（t-τ，s，n）和US（t-τ，s，n）轉(zhuǎn)換為笛卡爾坐標(biāo)系UP（x，xr，t0，t）和US（x，xr，t0，t），其表達(dá)式為

UP（x，xr，t0，t）=Re-iρ（s0）4π∫2π0dψ×

α（s0）ρ（s0）εα（s）ρ（s）Q（s）exp［-ω2（t-t0）］×expiω-（t-t0）+12P（s）Q（s）n2×t+nα（s）P（s）Q（s）n. （17）

US（x，xr，t0，t）=Re-iρ（s0）4π∫2π0dψ×

β（s0）ρ（s0）εβ（s）ρ（s）Q（s）exp［-ω2（t-t0）］×

expiω-（t-t0）+12P（s）Q（s）n2×n+nβ（s）P（s）Q（s）t. （18）

則此時(shí)G（x，xr，t）的高頻近似形式為

G（x，xr，t）=κP（xr）［UP（x，xr，t0，t）］T+κS（xr）［US（x，xr，t0，t）］T."" （19）

其中

κP（xr）=［sinψcosψ］T，

κS（xr）=［cosψsinψ］T.

式中，ψ為式（7）和（8）中的出射角；κP（xr）和κS（xr）為接收點(diǎn)xr處的應(yīng)力方向。

在水平地表，即m=（0，0，-1），將式（15）、（17）～（19）代入式（14）中得到檢波點(diǎn)波場

uU（x，xr，t）=iω∫t2t1dt∫SUT（x，xr，t0，t）Φu（xr，t0）dxr.（20）

其中，U（x，xr，t0，t）為式（17）和（18）中的時(shí)空域位移矢量，其表達(dá)式為

U（x，xr，t0，t）=UPx（x，xr，t0，t）USx（x，xr，t0，t）

UPz（x，xr，t0，t）USz（x，xr，t0，t）.（21）

Φ為二階2×2系數(shù)矩陣，其表達(dá)式為

Φ11=ρ（xr）β2（xr）［pPz（xr）sinψ+pPx（xr）cosψ］，

Φ12=ρ（xr）α2（xr）［（1-2

γ2）pPx（xr）sinψ+pPz（xr）cosψ］，

Φ21=ρ（xr）β2（xr）［pSz（xr）cosψ+pSx（xr）sinψ］，

Φ22=ρ（xr）α2（xr）［（1-2γ2）pSx（xr）cosψ+pSz（xr）sinψ］.（22）

式中，pPx和pPz為P波時(shí)空域高斯束慢度的水平和垂直分量；pSx和pSz為P波時(shí)空域高斯束慢度的水平和垂直分量；γ=β（xr）α（xr）為檢波點(diǎn)處的縱橫波速度比。

根據(jù)Yan的研究［24］，使用散度、旋度算子可以在矢量波場中分離出純縱波和純橫波分量，其表達(dá)式為

uPs（x，t）=α2（x）［·uD（x，xs，t）］，（23）

uPr（x，t）=α2（x）［·uU（x，xr，t）］，（24）

uSr（x，t）=-β2（x）×［×uU（x，xr，t）］.（25）

式中，uPs（x，t）、uPr（x，t） 和uSr（x，t）分別為震源純縱波分量、檢波點(diǎn)純縱波分量和檢波點(diǎn)純橫波分量。

1.3 成像條件和極性反轉(zhuǎn)

為了方便成像結(jié)果的對比，本文中使用震源歸一化成像條件，成像公式為

IPP（x）=∑xs∫T2T1uPs（x，t）uPr（x，t）dt∫T2T1uPs（x，t）uPr（x，t）dt ，（26）

IPS（x）=∑xs∫T2T1uPs（x，t）uSr（x，t）dt∫T2T1uPs（x，t）uPr（x，t）dt .（27）

當(dāng)對轉(zhuǎn)換波成像時(shí)，轉(zhuǎn)換波P-S成像可以顯示由于極化方向不同而導(dǎo)致的偏移剖面極性反轉(zhuǎn)。圖1展示了P-S轉(zhuǎn)換波的形成原理［25］。α為P波相對于偏移傾角的入射角。

在圖2中，以逆時(shí)針和向右為正方向的情況下，當(dāng)入射P波到達(dá)O1和O2時(shí)，由于α1lt;0和α2gt;0，在點(diǎn)O1產(chǎn)生的P-S轉(zhuǎn)換波具有負(fù)水平位移分量，而在點(diǎn)O2產(chǎn)生的P-S轉(zhuǎn)換波具有正水平位移分量，進(jìn)而導(dǎo)致在點(diǎn)R1和R2處接收的水平分量地震記錄具有相反的極性。因此根據(jù)P波的正入射角或負(fù)入射角α來校正成像結(jié)果?？紤]極性反轉(zhuǎn)，將式（23）～（25）代入式（26）和（27）得到最終的成像公式

IPP（x）=∑xs∫T2T1α2（x）［·uD（x，xs，t）］α2（x）［·uU（x，xr，t）］dt∫T2T1α2（x）［·uD（x，xs，t）］α2（x）［·uD（x，xs，t）］dt=

∑xs∫T2T1α4（x）［·uD（x，x，t）］［·uU（x，xr，t）］dt∫T2T1α4（x）［·uD（x，xs，t）］［·uD（x，xs，t）］dt .

IPS（x）=∑xs∫T2T1α2（x）［·uD（x，xs，t）］{-β2（x）×［×uU（x，xr，t）］}dt∫T2T1α2（x）［·uD（x，xs，t）］α2（x）［·uD（x，xs，t）］dtsgn（α）=

∑xs∫T2T1-α2β2（x）［·uD（x，xs，t）］×［×uU（x，xr，t）］dt∫T2T1α4（x）［·uD（x，x，t）］［·uD（x，xs，t）］dtsgn（α）.

1.4 黏彈介質(zhì)中的衰減補(bǔ)償

在黏聲介質(zhì)中的復(fù)值速度表達(dá)式為

V（x）=V0（x）1+i2Q-1（x）+1πQ-1（x）ln（ωωr）.

式中，V0（x）為聲波介質(zhì)中的速度；V（x）為黏聲介質(zhì)中的速度；Q為質(zhì)量因子；ω為圓頻率；ωr為參考頻率。

相似的，在黏彈介質(zhì)中的P波和S波復(fù)值速度表達(dá)式為

VP（x）=VP01+i2Q-1P（x）+1πQ-1P（x）lnωωr，

VS（x）=VS01+i2Q-1S（x）+1πQ-1S（x）lnωωr.

式中，VP（x）和VS（x）為黏彈介質(zhì)中的復(fù)速度；VP0和VS0為彈性波介質(zhì)中的速度；QP和QS為品質(zhì)因子。使用復(fù)值速度可以推導(dǎo)出黏彈介質(zhì)中的復(fù)值時(shí)間：

TP（x）=TP0（x）-i2T′P（x）-1πT′P（x）lnωωr，

TS（x）=TS0（x）-i2T′S（x）-1πT′S（x）lnωωr.（28）

其中

T′P（x）=∫ray1VPQPds，

T′S（x）=∫ray1VSQSds.

式中，TP（x）和TS（x）為黏彈介質(zhì)中的復(fù)值時(shí)間;TP0（x）和TS0（x）為彈性波介質(zhì)中的復(fù)值時(shí)間;T′P（x）和T′S（x）為造成衰減和頻散的實(shí)值因子；∫ray

為沿射線路徑的積分。

在震源波場傳播過程中，式（28）右側(cè)第二項(xiàng)隨著射線的時(shí)間向前傳播而產(chǎn)生指數(shù)振幅衰減，反轉(zhuǎn)這些項(xiàng)的符號將通過應(yīng)用指數(shù)增益來補(bǔ)償振幅衰減。式（28）右側(cè)第三項(xiàng)產(chǎn)生相位畸變，其效應(yīng)也通過反轉(zhuǎn)這些項(xiàng)的符號來補(bǔ)償。對于檢波點(diǎn)波場的傳播，也需要對式（28）的右側(cè)二三項(xiàng)反轉(zhuǎn)符號進(jìn)行補(bǔ)償。因此需要對震源波場和檢波點(diǎn)波場由于黏滯性存在造成的振幅能量損失和相位畸變進(jìn)行補(bǔ)償，這種補(bǔ)償是通過反轉(zhuǎn)式（28）右側(cè)第二項(xiàng)和第三項(xiàng)的符號實(shí)現(xiàn)的。

在黏彈介質(zhì)時(shí)空域高斯波包偏移中，需要對震源波場和檢波點(diǎn)波場都進(jìn)行補(bǔ)償，將震源波場和檢波點(diǎn)波場的衰減和復(fù)值走時(shí)聯(lián)系起來。使用時(shí)空域高斯束表征的格林函數(shù)構(gòu)建的震源波場和檢波點(diǎn)波場可以通過將式（28）中帶有TP′（x）、TS′（x）項(xiàng)反號進(jìn)行高斯束復(fù)值走時(shí)校正達(dá)到對地震波場補(bǔ)償?shù)哪康?，最后通過震源歸一化成像條件成像，得到黏彈介質(zhì)時(shí)空域高斯波包偏移方法。

1.5 穩(wěn)定性解決方法

目前補(bǔ)償吸收衰減的偏移方法最主要的困難在于補(bǔ)償振幅衰減和校正相位畸變過程中的不穩(wěn)定問題，穩(wěn)定性問題產(chǎn)生的根本原因在于高頻補(bǔ)償算子的指數(shù)增長造成的不穩(wěn)定現(xiàn)象。為了解決穩(wěn)定性問題，

本文中以P波為例來介紹解決穩(wěn)定性問題的方法。式（28）中，-i2T′P（x）和-1πT′P（x）ln （ωωr）分別表示能量衰減和相位畸變。在弱黏滯性介質(zhì)中（1/Q≤1），黏滯性不會(huì)改變射線路徑，而是通過頻率相關(guān)的耗散函數(shù)影響波形［26］。

D（x）=exp-ω2T′P（x）exp

iωπT′P（x）lnωωr.

其中，exp-ω2T′P（x）造成能量衰減，exp

iωπT′P（x）lnωωr造成相位畸變。為了解決穩(wěn)定性問題，本文中應(yīng)用

一個(gè)平滑最大限制增益函數(shù)［27］代替能量補(bǔ)償因子expω2T′P（x），該增益函數(shù)表達(dá)式為

φ（η）=exp（η），η≤lnG，

G［1-lnG-2.5（lnG）2］+G（1+5lnG）η-

2.5Gη2，ln（G）lt;η≤ln（G）+0.2，

1.1G，ηgt;lnG+0.2.

其中

η=ω2T′P（x）.

式中，G為預(yù)定義增益限制。為了比較能量補(bǔ)償因子exp

-ω2T′P（x）和增益函數(shù)φ（η），本文中設(shè)定G=3 000，圖2展示了兩個(gè)函數(shù)的形態(tài)?？梢钥闯?，當(dāng)η增加時(shí)，φ（η）不會(huì)無限增加，但會(huì)接近常數(shù)，也就是說，φ（η）會(huì)遵循式（28）的指數(shù)增益增長，直至接近預(yù)設(shè)定的最大增益，然后平穩(wěn)地過渡到最大增益。這樣，限制了高頻補(bǔ)償算子的指數(shù)的無限增長，解決了穩(wěn)定性問題。

2 模型測試

2.1 斷層模型

首先使用一個(gè)較為簡單的斷層模型驗(yàn)證本文方法的正確性。如圖3所示，該模型底部高速巖體向上上侵，造成了多個(gè)斷裂系統(tǒng)的出現(xiàn)，有較強(qiáng)的橫向速度變化；該模型淺層Q較小，有利于測試本文方法的補(bǔ)償作用。模型大小為1201×501，水平和垂直網(wǎng)格間距均為10 m。分別使用彈性波和黏彈性高階交錯(cuò)網(wǎng)格有限差分正演模擬方法生成炮記錄。在自由表面激發(fā)主頻為15 Hz的雷克子波構(gòu)建的爆炸P波震源。觀測系統(tǒng)采用移動(dòng)接收，共81炮，炮點(diǎn)間距為100 m，每個(gè)炮點(diǎn)有400個(gè)接收點(diǎn)，接收點(diǎn)間距為10 m。采樣時(shí)間為5.5 s，時(shí)間間隔為0.8 ms。圖4展示了彈性波和黏彈正演模擬方法第41炮水平x分量和垂直z分量。

圖5為彈性波介質(zhì)高斯束偏移結(jié)果、黏彈介質(zhì)高斯束偏移結(jié)果和黏彈介質(zhì)時(shí)空域高斯波包偏移結(jié)果對比?？梢钥闯?，3種方法都能對模型的斷層構(gòu)造進(jìn)行準(zhǔn)確的恢復(fù)，使繞射波收斂，成像位置準(zhǔn)確，成像分辨率較高。通過與參考剖面的對比，傳統(tǒng)高斯束方法和本文方法都能對中深層構(gòu)造進(jìn)行準(zhǔn)確成像，同相軸能量也得到了較好的恢復(fù)。為了更好地對比兩種偏移方法，從圖5中提取單道進(jìn)行能量對比。圖6展示了單道能量對比，可以看出，由于淺層的衰減作用，中深層的能量較弱，但是兩種偏移方法的能量與參考剖面幾乎相同，兩種方法都能夠?qū)φ穹M(jìn)行有效的恢復(fù)。圖7展示了參考剖面、傳統(tǒng)高斯束方法和本文方法的平均波數(shù)譜對比，3條曲線相近，兩種偏移方法的分辨率大致相同，也證明了時(shí)空域高斯波包偏移方法的正確性。黏彈介質(zhì)時(shí)空域高斯波包偏移相比于普通的高斯束偏移方法，其計(jì)算效率方面有較大優(yōu)勢。分別統(tǒng)計(jì)傳統(tǒng)高斯束方法和本文方法的運(yùn)行時(shí)間，兩種方法使用相同的運(yùn)行環(huán)境，傳統(tǒng)高斯束方法的運(yùn)行時(shí)間為2 821 s，本文方法的運(yùn)行時(shí)間為1 519 s，可以看出，本文方法能一定程度地提高計(jì)算效率。綜上，本文方法相比于傳統(tǒng)高斯束方法，成像效果接近，但本文方法的計(jì)算效率較高，證明了本文方法的正確性。

2.2 鹽丘模型

使用一個(gè)復(fù)雜的鹽丘模型測試本文方法的適應(yīng)性。如圖8所示，該模型有較大的鹽丘侵入體，具有高陡的斷層系統(tǒng)；其中，質(zhì)量因子模型頂部的Q值為20，有比較強(qiáng)的衰減作用。模型大小為1 290×300，水平和垂直網(wǎng)格間距均為10 m。仍然分別使用彈性波和黏彈性高階交錯(cuò)網(wǎng)格有限差分正演模擬方法生成炮記錄。在自由表面激發(fā)主頻為15 Hz的雷克子波構(gòu)建的爆炸P波震源。觀測系統(tǒng)采用移動(dòng)接收，共90炮，炮點(diǎn)間距為100 m，每個(gè)炮點(diǎn)有400個(gè)接收點(diǎn)，接收點(diǎn)間距為10 m。采樣時(shí)間為40 s，時(shí)間間隔為0.8 ms。圖9展示了彈性波和黏彈正演模擬方法第40炮水平x分量和垂直z分量，可以看出該炮記錄右側(cè)反射、繞射波場交疊混雜，反映出了鹽丘模型的復(fù)雜性。

圖10為彈性波介質(zhì)高斯束偏移結(jié)果、黏彈介質(zhì)高斯束偏移結(jié)果和黏彈介質(zhì)時(shí)空域高斯波包偏移結(jié)果對比?？梢钥闯?，參考剖面、傳統(tǒng)的高斯束方法和本文方法得到的剖面層位都比較清晰，相比于傳統(tǒng)的高斯束方法，本文提出的時(shí)空域高斯波包方法對鹽丘頂端兩側(cè)更加清晰，分辨率更高（紅色箭頭所示）。從圖10提取單道能量對比如圖11所示，可以看出參考剖面、傳統(tǒng)的高斯束方法和本文方法的能量保持效果近似，本文方法也能夠?qū)δ芰窟M(jìn)行較好的恢復(fù)。圖12為3種方法的10 km處的角道集對比，圖中的同相軸都被拉直，證明了本文方法的成像有效性和正確性，也證明了時(shí)空域高斯波包偏移方法對復(fù)雜模型的適應(yīng)性。同時(shí)，對比了傳統(tǒng)的高斯束方法和本文方法的計(jì)算效率。分別統(tǒng)計(jì)傳統(tǒng)高斯束方法和本文方法的運(yùn)行時(shí)間，兩種方法使用相同的運(yùn)行環(huán)境，傳統(tǒng)高斯束方法的運(yùn)行時(shí)間為1 637 s，本文方法的運(yùn)行時(shí)間為1 125 s，可以看出，本文方法在提高計(jì)算效率方面有一定的優(yōu)勢。為了測試噪聲對本文方法的影響，在合成地震記錄中加入高頻噪聲，圖13展示了加入高頻噪聲后的第40炮的水平x分量和垂直z分量，可以看到整個(gè)炮記錄上有大量的高頻噪聲。

圖14展示了未加入噪聲和加入噪聲后的偏移成像結(jié)果。通過對比，可以看出，加入噪聲后的偏移結(jié)果整個(gè)剖面存在大量的噪聲，但是同相軸基本歸位，層位較為清晰。因此合成地震記錄中的噪聲會(huì)使偏移剖面中存在大量的噪聲，但不會(huì)影響關(guān)鍵層位的成像，說明本文方法對存在噪聲的地震記錄也具有適應(yīng)性。綜上所述，本文方法對含有噪聲的復(fù)雜構(gòu)造有較強(qiáng)的適應(yīng)性，能夠?qū)?fù)雜構(gòu)造進(jìn)行準(zhǔn)確的成像，而且相比于傳統(tǒng)的高斯束方法，計(jì)算效率更高。

3 結(jié)束語

針對傳統(tǒng)高斯束方法頻繁的傅里葉變換造成的計(jì)算效率低以及黏滯性造成的地震波傳播過程中的能量損失和相位畸變的問題，本文

中將時(shí)空域高斯束擴(kuò)展到黏彈介質(zhì)，同時(shí)考慮黏滯性對高斯束復(fù)值走時(shí)的影響，補(bǔ)償了黏滯性造成的波場衰減，發(fā)展了黏彈介質(zhì)時(shí)空域高斯波包偏移方法。該方法能夠處理多分量黏性地震記錄，消除黏滯性造成的負(fù)面影響，充分利用彈性波場的矢量特性，在保持成像精度的同時(shí)，避免了頻繁的傅里葉變換，提高了計(jì)算效率。斷層模型和鹽丘模型的測試表明，本文方法能夠?qū)δＰ瓦M(jìn)行準(zhǔn)確的成像，對復(fù)雜構(gòu)造也具有良好的適應(yīng)性，提高了計(jì)算效率。

參考文獻(xiàn)：

［1］ 黃建平，張晴，張凱，等.格林函數(shù)高斯束逆時(shí)偏移［J］.石油地球物理勘探，2014，49（1）：101-106.

HUANG Jianping， ZHANG Qing， ZHANG Kai， et al. Greens function Gaussian beam inverse time migration［J］.Oil Geophysical Prospecting，2014，49（1）：101-106.

［2］ HILL N R. Gaussian beam migration［J］. Geophysics， 1990，55（11）：1416-1428.

［3］ GRAY S H， BLEISTEIN N. True-amplitude Gaussian-beam migration［J］. Geophysics， 2009，74（2）：S11-S23.

［4］ 李振春，岳玉波，郭朝斌，等.高斯波束共角度保幅深度偏移［J］.石油地球物理勘探，2010（3）：360-365.

LI Zhenchun， YUE Yubo， GUO Chaobin， et al. Amplitude preserved common angle Guassian beam depth migration［J］. Oil Geophysical Prospecting， 2010（3）：360-365.

［5］ 岳玉波，李振春，劉偉，等.保幅炮域高斯波束偏移［J］.中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2011，35（1）：52-55.

YUE Yubo， LI Zhenchun，" LIU Wei， et al. Preserved amplitude shot domain Gaussian beam migration［J］. Journal of China University of Petroleum （Edition of Natural Sciences）， 2011，35（1）：52-55.

［6］ 黃建平，袁茂林，段新意，等.一種解耦的起伏地表彈性波高斯束偏移方法［J］.石油地球物理勘探，2015，50（3）：460-468.

HUANG Jianping， YUAN Maolin， DUN Xinyi， et al. Decoupled elastic Gaussian beam migration for rugged topography［J］. Oil Geophysical Prospecting， 2015，50（3）：460-468.

［7］ 黃建平，楊繼東，李振春，等.基于有效鄰域波場近似的起伏地表保幅高斯束偏移［J］.地球物理學(xué)報(bào)，2016，59（6）：2245-2256.

HUANG Jianping， YANG Jidong， LI Zhenchun， et al. An amplitude-preserved Gaussian beam migration based on wave field approximation in effective vicinity under irregular topographical conditions［J］.Chinese Journal of Geophysics， 2016，59（6）：2245-2256.

［8］ 楊繼東.基于動(dòng)態(tài)參數(shù)控制的地震波束偏移方法研究［D］.青島：中國石油大學(xué)（華東），2016.

YANG Jidong. The research of the seismic beam migration limited by dynamic parameters［D］. Qingdao： China University of Petroleum（East China），2016.

［9］ 呂慶達(dá).基于子波重構(gòu)的時(shí)空域高斯束偏移方法研究［D］.青島：中國石油大學(xué)（華東），2019.

L Qingda. The research of prestack time-space Gaussian beam depth migration based on wavelet reconstruction［D］. Qingdao： China University of Petroleum （East China），2019.

［10］ KACHALOV A P， POPOV M M. Application of the method of summation of Gaussian beams to the calculation of theoretical seismograms［J］. Journal of Soviet Mathematics， 1990，50（4）：1727-1742.

［11］ GENG Y， WU R S， GAO J H. Efficient Gaussian packets representation and seismic imaging［C/OL］. SEG Technical Program Expanded Abstracts，2011［2022-10-05］. https：//library.seg.org/doi/abs/10.1190/1.3627903.

［12］ GENG Y， WU R S， GAO J H. Gabor-frame-based Gaussian packet migration［J］. Geophysical Prospecting， 2014，62（6）：1432-1452.

［13］ 李輝，馮波，王華忠.波場模擬的高斯波包疊加方法［J］.石油物探，2012，51（4）：327-337.

LI Hui， FENG Bo， WANG Huazhong. A new wave field modeling method by using Gaussian Packets［J］. Geophysical Prosepecting for Petroleum， 2012，51（4）：327-337.

［14］ 李輝，王華忠，馮波，等.特征高斯波包疊前深度偏移方法［J］.地球物理學(xué)報(bào)，2014，57（7）：2258-2268.

LI Hui， WANG Huazhong， FENG Bo， et al. Efficient pre-stack depth migration method using characteristic Gaussian packets［J］. Chinese Journal of Geophysics， 2014，57（7）：2258-2268.

［15］ 李輝，王華忠.基于高斯束與高斯波包的Gabor框架散射波模擬方法［J］.地球物理學(xué)報(bào)，2015，58（4）：1317-1332.

LI Hui， WANG Huazhong. A scattering wave modeling method using Gaussian beam and Gaussian packet in the Gabor domain［J］. Chinese Journal of Geophysics， 2015，58（4）：1317-1332.

［16］ L Q， HUANG J， YANG J， et al. An optimized time-space Gaussian beam migration method with dynamic parameter control［J］. Journal of Applied Geophysics， 2019，160：47-56.

［17］ ZHANG D L， HUANG J P， YANG J D， et al. A fast time-space-domain Gaussian beam migration approach using the dominant frequency approximation［J］. Petroleum Science， 2022，19（4）：1555-1565.

［18］ 張巖，黃曉瑩.巖心孔隙結(jié)構(gòu)的空間分辨核磁表征［J］.中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2021，45（4）：49-58.

ZHANG Yan， HUANG Xiaoying. Spatially resolved NMR for rock core analysis ［J］.Journal of China University of Petroleum（Edition of Natural Science）， 2021，45（4）：49-58.

［19］ 劉學(xué)鋒，張曉偉，曾鑫，等.采用機(jī)器學(xué)習(xí)分割算法和掃描電鏡分析頁巖微觀孔隙結(jié)構(gòu)［J］.中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2022，46（1）：23-33.

LIU Xuefeng， ZHANG Xiaowei， ZENG Xin， et al. Pore structure characterization of shales using SEM and machine learning based segmentation method［J］.Journal of China University of Petroleum（ Edition of Natural Science），2022，46（1）：23-33.

［20］ BICKEL S H， NATARAJAN R R. Plane-wave Q deconvolution［J］. Geophysics， 1985，50（9）：1426-1439.

［21］ ZHANG J F， WAPENAAR K. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media［J］. Geophysical Prospecting， 2002，50（6）：629-643.

［22］ ZHANG J F， WU J Z， LI X Y. Compensation for absorption and dispersion in prestack migration：an effective Q approach［J］. Geophysics， 2013，78（1）：S1-S14.

［23］ PAO Y H， VARATHARAJULU V. Huygens principle， radiation conditions， and integral formulas for the scattering of elastic waves［J］. The Journal of the Acoustical Society of America， 1976，59（6）：1361-1371.

［24］ YAN J ， SAVA P. Isotropic angle-domain elastic reverse-time migration［J］. Geophysics， 2008，73（6）：S229-S239.

［25］ DU Q， ZHU Y， BA J. Polarity reversal correction for elastic reverse time migration［J］. Geophysics， 2012，77（2）：S31-S41.

［26］ KEERS H， VASCO D W， JOHNSON L R. Viscoacoustic crosswell imaging using asymptotic waveforms［J］. Geophysics， 2001，66（3）：861-870.

［27］ ZHANG J， WU J， LI X. Compensation for absorption and dispersion inprestack migration： an effective Q approach［J］. Geophysics， 2013，78（1）：S1-S14.

（編輯 修榮榮）