白敏， 陳小宏， 吳娟， 陳陽康， 劉國昌， 王恩江

1 中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室， 北京 102249 2 中國石油大學(xué)(北京)海洋石油勘探國家工程實驗室， 北京 102249 3 華北水利水電大學(xué)資源與環(huán)境學(xué)院， 鄭州 450045 4 Bureau of Economic Geology，John A．a(chǎn)nd Katherine G．Jackson School of Geosciences，The University of Texas at Austin，University Station，Box X，Austin，TX 78713-8924，USA

?

基于吸收衰減補償?shù)亩喾至扛咚故鏁r偏移

白敏1,2,3， 陳小宏1,2， 吳娟1,2,3， 陳陽康4， 劉國昌1， 王恩江1

1 中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室， 北京 102249 2 中國石油大學(xué)(北京)海洋石油勘探國家工程實驗室， 北京 102249 3 華北水利水電大學(xué)資源與環(huán)境學(xué)院， 鄭州 450045 4 Bureau of Economic Geology，John A．a(chǎn)nd Katherine G．Jackson School of Geosciences，The University of Texas at Austin，University Station，Box X，Austin，TX 78713-8924，USA

高斯束逆時偏移結(jié)合了射線類偏移的高計算效率和波動方程逆時偏移的高精度，能很好地處理焦散點、大傾角成像問題，并且具有面向目標(biāo)成像的能力．多分量地震資料的偏移技術(shù)可以對地下復(fù)雜構(gòu)造進行更準(zhǔn)確的成像，由于實際地下介質(zhì)具有黏滯性，研究黏彈性疊前逆時偏移具有一定的現(xiàn)實意義．本文采用高斯束逆時偏移方法對多分量地震數(shù)據(jù)進行吸收衰減補償，首先分別給出縱波和轉(zhuǎn)換波共炮域高斯束疊前逆時偏移方法原理，在此基礎(chǔ)上推導(dǎo)補償吸收衰減的表達(dá)式，校正Q引起的振幅衰減和相位畸變，實現(xiàn)基于吸收衰減補償?shù)亩喾至扛咚故B前逆時偏移．?dāng)?shù)值模型的測試結(jié)果顯示，在考慮地下介質(zhì)的黏滯性時，本文方法具有更高的成像分辨率．

衰減補償； 多分量； 高斯束； 逆時偏移； 格林函數(shù)

Multiple-component seismic data contains both PP and PS waves. The PS-wave image complements the traditional PP-wave image，resulting in a more accurate subsurface characterization．Reverse-time migration of multiple-component seismic data can improve the accuracy of imaging subsurface complex geological structures. While viscoelastic prestack reverse-time migration is of practical significance because it considers the viscosity of subsurface media．As a new migration tool，Gaussian beam reverse-time migration (GBRTM) combines the high efficiency and flexibility of Gaussian beam migration with the high accuracy of wave equation reverse-time migration，which can overcome the problems of caustics，handle all arrivals，yield good images of steep flanks，and is easy to extend to target-oriented implementation．However，GBRTM studies have focused on acoustic waves, and multiple-component GBRTM has been little investigated．Besides，it is not clear how the method should be applied for multiple-component seismic data recorded in attenuating media．Therefore，we propose a multiple-component GBRTM to perform seismic data compensation for frequency-dependent absorption and dispersion．We separate multiple-component seismic data into PP- and PS-waves，and migrate by scalar migration methods．The purpose is to provide a new effective method for multiple-component seismic data migration imaging and to compensate the attenuation simultaneously．First，we derive a common-shot gathers GBRTM algorithm of PP and PS waves．Then，the expressions of attenuation equation and the precision analysis of Green function based on Gaussian beams are developed. Finally we present the principle and procedures of compensation，and then propose an attenuation-compensated multiple-component GBRTM．

The migration results of PP and PS waves illustrate that the new method is effective in compensating the amplitude loss and phase shift caused by the anelastic properties of rocks in the field．The migration results have higher amplitudes and more continuous reflectors，especially in deep sections．Comparison of single trace waveforms extracted from migration results shows that the proposed approach effectively compensates the absorption of the subsurface medium．From the amplitude spectra and power spectra，we see the new method effectively compensates the seismic wave energy，and especially enhances the energy of the middle- and high- frequency components．

We propose a attenuation-compensated multiple component method based on GBRTM to compensate the energy and correct the phase in the seismic wave migration．Compared to the Gaussian beam prestack depth migration proposed by Hill，GBRTM is superior in theory because it does not require local slant stack and steepest-descents evaluation．The new attenuation compensation method is an attractive migration algorithm，because it not only has the advantages of the high computational efficiency of ray-basedQ-compensated migration，but also retains the high accuracy of the attenuation compensation method based on wave equation reverse-time migration．We have also demonstrated that the images obtained by the new method can compensate the attenuation and dispersion effects．Numerical results further verify that the proposed approach can effectively improve the resolution and quality of migrate images for both PP- and PS waves，particularly beneath high-attenuation zones．

1 引言

在地震成像早期，偏移主要是為了得到地下地質(zhì)體的構(gòu)造信息．現(xiàn)在獲取地下介質(zhì)真反射系數(shù)已經(jīng)引起了地震勘探界的廣泛關(guān)注，這就要求在偏移延拓和成像中對影響成像振幅的因素做校正處理(Zhang et al．，2002)．然而，固有的吸收因子在常規(guī)偏移方法中經(jīng)常被忽略．因為深層高頻損失，地震波穿過覆蓋層時，介質(zhì)的空間變化引起地震振幅衰減，子波相位畸變，成像分辨率降低，使得異常區(qū)(比如氣層)振幅衰減，頻帶變窄，造成深層識別和解釋的困難，也影響了準(zhǔn)確預(yù)測儲層的能力．因此，需要補償因為這類衰減引起的吸收效應(yīng)．

在疊后Q偏移方面，Wang和Guo (2004)提出了偏移加反Q濾波算法，并給出了穩(wěn)定的偏移算子，但這種算子只適用于一維模型．Wang (2008)將一維模型推廣到二維，提出了適用于速度和Q值隨空間變化的Q偏移算法．在疊前Q偏移方面，Traynin等(2008)利用Kirchhoff 疊前深度Q偏移對氣藏下方的區(qū)域進行振幅和頻帶恢復(fù)；Xie等(2009)發(fā)展了疊前KirchhoffQ偏移，利用射線追蹤計算沿射線的吸收效應(yīng)，并在偏移中對每一個頻帶進行補償，但這些方法都是基于Kirchhoff．由于射線理論的缺陷(比如在焦散區(qū)、陰影區(qū)失效)，該方法不能處理復(fù)雜的地質(zhì)體．為解決這一問題，Xie等(2010)將Kirchhoff疊前深度Q偏移拓展到高斯束．Xiao等(2014)進一步提出共偏移距激光束算法，通過限制束的發(fā)散，使其傳播類似激光，實現(xiàn)疊前激光束Q偏移．對波動方程Q偏移，Mittet(2007)提出了補償振幅衰減和校正相位的策略，并且證明在Q有10%誤差的情況下，Q補償?shù)钠埔廊荒艿玫捷^好的成像效果，相較于未補償?shù)钠拭?，分辨率明顯提高．Zhang和Wapenaar(2002)通過限制延拓步數(shù)和最大偏移角，提出了條件穩(wěn)定的延拓算子，但是該方法只能實現(xiàn)地下介質(zhì)有限傾角和深度的成像．Valenciano等(2011)利用傅里葉有限差分法進行波場延拓，在衰減介質(zhì)中引入一種新的波動方程偏移方法．通過改進的標(biāo)量方程，Deng 和 McMechan(2007)提出時間域Q補償?shù)膹椥阅鏁r偏移，然而，該方法只補償振幅衰減，而忽略了相位校正．基于吸收系數(shù)和頻率滿足近似線性關(guān)系的假設(shè)，Zhang等(2010)推導(dǎo)了常Q模型的波動方程，在疊前逆時偏移中補償?shù)卣饠?shù)據(jù)的衰減效應(yīng)，Suh等(2012)將該方法拓展到各向異性介質(zhì)．Fletcher等(2012)通過兩次波傳播來估算沿傳播路徑的衰減走時，利用該走時對常規(guī)的震源和檢波點波場進行濾波，在偏移前補償Q的影響．基于標(biāo)準(zhǔn)的線性固體模型假設(shè)，Bai 等(2013)推導(dǎo)了時空域的黏滯波動方程，方程中包含偽微分算子項來模擬黏滯性，與Zhang 等(2010)相比，此方法更容易實現(xiàn)．Yan和Liu (2013)采用優(yōu)化的時空域高階有限差分實現(xiàn)黏滯波動方程的疊前逆時偏移，該方程具有高的差分精度．Zhu 等(2014)分析了逆時偏移中Q對震源波場和檢波點波場的影響，在黏滯波動方程波場延拓中對影響振幅和相位的算子進行解耦，實現(xiàn)Q逆時偏移．李振春等(2014)提出了黏滯波動方程的最小二乘逆時偏移，與未補償?shù)慕Y(jié)果相比，能得到更準(zhǔn)確的反射界面及更均衡的振幅值，然而該方法所需計算量較大．

多分量地震勘探能得到PP波和PS波數(shù)據(jù)，與單純的PP波相比，PS數(shù)據(jù)攜帶了目標(biāo)區(qū)域不同的信息，能提供更多的物性參數(shù)，有效改進地震勘探的準(zhǔn)確性．比如在PP波信號弱的區(qū)域，PS波可以增強照明．聯(lián)合PP波和PS波的多分量成像技術(shù)能夠更準(zhǔn)確地反映地下特征．然而，常規(guī)的地震數(shù)據(jù)處理只考慮記錄的P波分量，而忽略了轉(zhuǎn)換PS波．因此，為了充分利用多分量信息，需要更有效的多分量偏移技術(shù)．

在多分量偏移方面，基于彈性波的Kirchhoff偏移方法最先被提出(Kuo and Dai， 1984；Wapenaar et al．，1990；Hokstad, 2000；Du and Hou，2008)．這些方法的關(guān)鍵是計算PP波和PS波的走時，以此求和計算成像振幅．為了克服Kirchhoff成像的弱點，一些學(xué)者提出了彈性單程波動方程偏移(Fisk and McCartor，1991；Wu，1994；Xie and Wu，2005)．由于是雙程波的近似，單程波偏移方法受到成像角度、波的傳播路徑問題的限制，影響成像精度．因此，近年來利用雙程波動方程的多分量逆時偏移方法成為熱點．多分量的逆時偏移方法有兩類：彈性波偏移和標(biāo)量偏移(Hou and Marfurt，2002)．彈性波偏移是直接輸入多分量地震數(shù)據(jù)，利用彈性波動方程構(gòu)建矢量波場(Chang and McMechan，1994；Yan and Sava，2008；Du et al．，2012，2014)．此外，也有些研究者通過把波場分離成PP波和PS波，提出了多分量地震數(shù)據(jù)的標(biāo)量偏移方法(Wang and Nemeth，1997；Jin et al．，1998；Sun and McMechan，2001；Hou and Marfurt，2002；Sun et al．，2004，2006；Han et al．，2013)．

以上的方法中，射線類偏移具有高的效率和靈活性，但是不能處理復(fù)雜的地質(zhì)體．波動方程逆時偏移更精確，但是計算量大．Popov等(2007，2008，2010)結(jié)合高斯束的高效靈活性和波動方程逆時偏移的高精度，提出了高斯束逆時偏移，即以Kirchhoff積分為基礎(chǔ)，利用高斯束的疊加積分計算格林函數(shù)的逆時偏移算法．由于高斯束的疊加積分是嚴(yán)格按照數(shù)學(xué)理論推導(dǎo)的，因此避免了多值走時和不精確的振幅計算，另外可通過控制成像條件的時間窗來消除續(xù)至波，更適合基于目標(biāo)的成像．與Hill(2001)提出的高斯束疊前深度偏移方法相比，高斯束逆時偏移不做局部傾斜疊加和最速下降近似，在理論上要優(yōu)于前者．

目前對于高斯束偏移方法的研究主要集中于縱波，對多分量地震記錄的研究較少(Li and Mao，2015)．另外，高斯束逆時偏移作為一種新的偏移方法，對轉(zhuǎn)換波成像及吸收衰減補償?shù)男Ч胁磺宄虼?，本文在前人研究基礎(chǔ)上提出了基于吸收衰減補償?shù)亩喾至扛咚故鏁r偏移，首先把波場分解為PP波和PS波(李志遠(yuǎn)等，2013)，然后分別基于聲波方程的標(biāo)量方法偏移，目的是提供一個新的有效的多分量數(shù)據(jù)成像技術(shù)，同時補償?shù)卣鹩涗浀奈账p．文中首先給出聲波高斯束逆時偏移中的正向和反向延拓波場表達(dá)式，將聲波反向延拓波場格林函數(shù)拓展到彈性波，得到P波和S波的反向延拓波場，再將其與正向延拓波場互相關(guān)，推導(dǎo)波場分離后的PP波和PS波共炮域高斯束疊前逆時偏移表達(dá)式．在此基礎(chǔ)上引入與Q有關(guān)的補償算子，實現(xiàn)基于吸收衰減補償?shù)亩喾至扛咚故B前逆時偏移．最后用數(shù)值算例證明了本文方法的有效性和優(yōu)越性．

2 多分量高斯束逆時偏移

Popov 等(2010) 首先提出了高斯束逆時偏移的計算方法，并給出詳細(xì)的思路和推導(dǎo)過程．這里，我們只討論關(guān)鍵步驟．考慮標(biāo)量聲波方程，波場U(x,t)滿足：

(1)

其中，x=(x,z)是空間內(nèi)一點，Δ是拉普拉斯算子，v(x)是波傳播速度．

2.1 正向和反向延拓波場的構(gòu)建

f(t)δ(x-xs)； U(D)|t

(2)

式中f(t)是初始波(震源子波f(t)|t<0=0)．利用高斯束計算格林函數(shù)，則正向延拓波場可以表示為：

×GGB(x,xs;ω)，

(3)其中，fF(ω)是初始子波的傅里葉變換，GGB(x,xs;ω)是高斯束表示的格林函數(shù)的漸近式．

對于反傳波場，地下成像區(qū)域內(nèi)一點x0的格林函數(shù)G(x,t;x0,t0)滿足：

δ(t-t0)δ(x-x0); G|t

(4)

其中t0是瞬時時刻，滿足0≤t0≤T．利用Kirchhoff積分，得到x0點反傳波場U(x0,t0)的表達(dá)式：

(5)

式中：? Ω是閉合空間Ω的邊界，U(0)(x,t)是地震記錄，?/? nx是沿Ω外法線方向的導(dǎo)數(shù)，如圖1所示．

圖1 偏移域Ω的參數(shù)化xs和x0分別表示震源和地下成像區(qū)域(虛線內(nèi))任一點，v(x)是速度，nx是Ω在記錄表面的外法向?qū)?shù)，Ω的邊界向外延伸至無窮.Fig.1 Parametrization of the migration domain Ωxs and x0 are positions of the source and a point in the migration domain (refer to the dashed range)，respectively，v(x) is the velocity，and nx is the external normal to domain Ω on the recording surface．The boundary of Ω is extended outward to infinite.

此外，假設(shè)地球的水平面為邊界? Ω的一部分，則? Ω的兩側(cè)對波場的貢獻可以忽略．格林函數(shù)滿足Kirchhoff近似中的邊界條件G|z=0=0，則偏移域內(nèi)一點x0的反向延拓波場為：

(6)

2.2 高斯束積分表示的格林函數(shù)

高斯束逆時偏移中最重要的步驟是構(gòu)建高斯束表示的格林函數(shù)漸近式．格林函數(shù)GGB(x,xs,ω)是通過一系列由源點出射的，具有不同出射角的高斯束疊加積分表示的，如圖2．

圖2 高斯束表示的格林函數(shù)Fig.2 Sketch of Green function in terms of Gaussian beams

對于正向延拓波場，從源點xs處以不同的角度出射中心射線束，每條中心束附近波場值用高斯束方法求取，地下介質(zhì)中一點x的波場值由與其臨近的多條高斯束疊加而得，格林函數(shù)GGB(x,xs,ω)的表達(dá)式為(Hill, 2001，Gray and Bleistein, 2009)：

(7)

其中，uGB(x,xs,p,ω)為高斯束方法求取的波場位移：

(8)

式中，px和pz分別表示射線參數(shù)的水平分量和垂直分量，ω是頻率，A為振幅，T是復(fù)值走時．

反向延拓波場的格林函數(shù)為：

×GGB(x,x0,ω)，

(9)

其中，

(10)

式中，x0和x分別為偏移域內(nèi)任一點和地下介質(zhì)中計算點的位置．

對于P波，有

(11)

對于S波，有

(12)

所以，對成像區(qū)域內(nèi)一點x0，PP波和PS波的走時分別為：

(13)

(14)

式(8)和式(11)、(12)中高斯束的復(fù)值振幅A和走時T分別表示為(Hill,1990,2001)：

(15)

(16)

2.3 成像條件

對原始地震記錄進行波場分離后，采用適用于PP波和PS波的高斯束方法對其進行偏移．本文應(yīng)用互相關(guān)成像條件，則PP波和PS波高斯束逆時偏移公式為：

IPP(x0,xs)=∫dt0U(D)(x0,t0;xs)UPP(x0,t0)，

(17)

IPS(x0,xs)=∫dt0U(D)(x0,t0;xs)UPS(x0,t0)，

(18)

其中，U(D)(x0,t0;xs)和U(x0,t0)的表達(dá)式見方程(3)和(6)．

3 衰減介質(zhì)多分量高斯束逆時偏移

為了得到衰減補償?shù)亩喾至扛咚故鏁r偏移表達(dá)式，本小節(jié)首先給出衰減補償原理，然后分析高斯束在衰減介質(zhì)中傳播的精度．

3.1 衰減和補償原理

圖3 無衰減介質(zhì)與衰減介質(zhì)正演和偏移示意圖. (a)、(c)為無衰減介質(zhì)正演和偏移； (b)、(d)為衰減介質(zhì)正演和偏移Fig.3 Schematic of forward modeling (a) and migration extrapolation (c) in a non-attenuation medium，and forward modeling (b) and migration extrapolation (d) in an attenuating medium

3.2 高斯束在衰減介質(zhì)中的傳播精度

波在黏彈性介質(zhì)中傳播可認(rèn)為是復(fù)速度的波在彈性介質(zhì)中傳播，復(fù)速度的實部分別是彈性介質(zhì)中的速度vP和vS，品質(zhì)因子Q代表衰減．QP和QS分別是縱橫波的品質(zhì)因子，如果衰減小(1/Q?1)，則復(fù)速度為(Keers et al．，2001)：

(19)

(20)

式中ω0是參考頻率，Q不依賴于頻率，速度的虛部引起沿著射線路徑的指數(shù)振幅衰減，速度的實部包含頻散項，確保波動方程解的因果性．

考慮衰減(補償)的PP波和PS波的復(fù)值走時分別為：

(21)

(22)

其中，

(23)

(24)

由于本文是基于波場分離后的標(biāo)量方程計算，因此，下面以聲波和黏聲介質(zhì)為例，分析高斯束在衰減和無衰減情況下的波場傳播精度．

我們設(shè)計了一個均勻介質(zhì)模型，速度為2000m·s-1，模型大小為201×201，dx=dz=10m，震源使用主頻為20Hz的雷克子波，位于(1000m，1000m)處．分別用高斯束積分和頻率域波動方程來模擬無衰減介質(zhì)和衰減介質(zhì)中地震波的傳播，并與解析解對比，分析兩者的精度．

基于Futterman(1962)的理論，利用Valenciano等(2011)推導(dǎo)的頻率域波動方程：

(25)

其中，復(fù)速度

(26)

該式求解采用頻率域有限差分法(Joetal．，1996)．

(27)

在無衰減均勻介質(zhì)中，分別用高斯束積分和頻率域聲波方程合成格林函數(shù)，并將其與解析的格林函數(shù)進行比較(如圖4)．圖4a左是用頻率域有限差分法求解得到的格林函數(shù)(WE)，頻率為20Hz(下同)；圖4a右是用高斯束求解得到的格林函數(shù)(GB)．圖中可以看出：無論實部還是虛部，兩個波場形態(tài)基本一致．圖4b是從圖4a的單頻波場中，抽出的深度z=1000m處的值與解析格林函數(shù)的對比．曲線圖表明：當(dāng)傳播距離大于100m時，圖4a中單頻波場的實部、虛部均與解析解的實部、虛部吻合較好(震源附近的振幅差異是因為射線理論是一種高頻近似方法，屬于遠(yuǎn)場近似，震源是射線理論的奇異點)．這進一步證明了高斯束和頻率域波動方程求解的準(zhǔn)確性，也說明高斯束積分所計算的格林函數(shù)很好地近似了解析的格林函數(shù)．

圖5是在衰減的均勻介質(zhì)中，Q=50，30，10時，分別用高斯束積分和頻率域黏聲波方程合成的格林函數(shù)，以及它們的精度比較(為了更好地對比格林函數(shù)的細(xì)節(jié)，選取局部放大圖)．圖5a，5c，5e分別為兩種方法求得的單頻波場；圖5b，5d，5f分別是抽取的z=1000m處的值所做的精度對比．從兩者實部和虛部的對比可以看出，這兩種方法模擬得到的實部和虛部的振幅和相位均吻合較好，由此也驗證了高斯束表示的格林函數(shù)能較好地描述波在黏聲介質(zhì)中的傳播．

圖4 無衰減介質(zhì)中，高斯束積分和波動方程計算的格林函數(shù)(a)及其與解析格林函數(shù)的精度對比(b)Fig.4 Green function solved by wave equation and Gaussian beam in non-attenuation medium (a) and comparison to analytical Green function (b)

圖5 衰減介質(zhì)中，高斯束積分和波動方程計算的格林函數(shù)及其精度對比(局部放大圖)(a，b) Q=50； (c，d) Q=30； (e，f) Q=10.Fig.5 Green function solved by wave equation and Gaussian beam and accuracy comparison in attenuation medium (partial enlargement)

4 數(shù)值算例

為驗證基于吸收衰減補償?shù)亩喾至扛咚故B前逆時偏移方法的有效性，本文設(shè)計了兩個試驗．

4.1 兩層模型

模型大小為3000 m×3000 m，如圖6所示．震源采用雷克子波，其頻率為20 Hz．震源坐標(biāo)為(1500 m，0 m)，共設(shè)計了301個檢波點，分別布置在地表0 m到3000 m處，道間距為10 m，接收記錄長度為2.0 s，時間采樣間隔是0.001 s．

為了更直觀地反映補償過程，以地下一點(1500 m,1500 m)處的成像為例，分別輸出PP波與PS波震源和反傳波場的原始波場快照以及補償后的波場快照，如圖7所示．其中圖7a和7b為參考和補償后的震源波場，圖7c和7d為參考和補償后地下一點x0的PP波波場．由圖可以看出補償后的震源波場振幅加強，物理意義如圖3d所示，用于補償與exp(-αLD)有關(guān)的吸收．同樣，x0點的波場用以補償與exp(-αLB)有關(guān)的吸收，補償后的x0點波場與地震記錄相關(guān)，實現(xiàn)逆時偏移中波場的反傳過程．由于震源和反傳波場都得到了補償，因而相關(guān)之后x0點的成像值，振幅恢復(fù)，相位得到校正．對地下所有點成像之后的PP波單炮偏移結(jié)果如圖8所示．相似地，圖7e和7f為參考和補償后地下一點x0的PS波波場，對地下所有點成像之后的PS波單炮偏移剖面見圖9．

圖6 兩層模型Fig.6 Two-layer model

圖8和圖9分別是三種情況下(參考，未補償，補償后)的PP波和PS波單炮偏移剖面．從圖中可以看出，基于吸收衰減補償?shù)母咚故B前逆時偏移考慮了地層的吸收效應(yīng)，在偏移過程中對地震記錄的衰減能量進行了補償，補償后的偏移剖面更接近理想的偏移剖面．圖10是從圖8和圖9中抽取的單道波形．通過對比發(fā)現(xiàn)，本文方法有效地補償了地下介質(zhì)對地震波的吸收．圖11為圖8和圖9對應(yīng)的振幅譜和功率譜，從頻譜對比可以看到，補償后的剖面能量得到恢復(fù)，尤其是中高頻能量成分得到了加強．

4.2 氣云模型

在第二個例子中，我們考慮實際的衰減模型．圖12是速度模型和相應(yīng)的Q模型(BP模型的一部分)．模型的頂部偏中是由于氣云的存在導(dǎo)致的低速和高衰減區(qū)．橫波速度和Q值由縱波速度與Q值按比例得到，vP/vS=1.73，QP/QS=1.2，模型網(wǎng)格數(shù)為398×161，網(wǎng)格大小是10 m×10 m．

我們采用彈性和黏彈性方程模擬地震記錄 (Zhu and Carcione，2014)，共40炮，震源采用主頻為15 Hz的雷克子波，炮間距為100 m，每一炮都是全地表等距(10 m)接收．接收記錄長度為3 s，時間采樣間隔為0.001 s．

取其中一炮的記錄，如圖13和圖14所示，炮點位于1000 m處．從圖中可以看出，不管是x分量還是z分量，黏彈性記錄中，反射波的同相軸能量均相對較弱，且同相軸的波形要寬，介質(zhì)的黏滯性對地震波能量的吸收和衰減作用明顯．因為經(jīng)過高衰減區(qū)的反射消失了，所以黏彈性記錄更干凈．

圖15是從圖13和圖14中抽出的第101道(橫向距離1000 m)的波形對比，從圖中也可以看出，介質(zhì)的黏滯性對地震波的衰減作用明顯．圖16是圖13和圖14對應(yīng)的振幅譜和功率譜．

圖7 波到達(dá)界面之后某一時刻的原始波場快照和補償后的波場快照震源波場：(a) 參考，(b) 補償后；地下一點 (1500 m,1500 m)的PP波波場：(c)參考，(d)補償后；PS波波場：(e)參考，(f)補償后.Fig.7 Snapshots of reference wavefield and compensated wavefield at the time when waves arrives at the interfaceSource wavefield： (a) Reference; (b) Compensated. PP-wavefield：(c) Reference; (d) Compensated; PS-wavefield：(e) Reference; (f) Compensated at the point x0 (1500 m，1500 m) in the subsurface．

圖8 PP波單炮偏移(a) 參考; (b) 未補償; (c) 補償后.Fig.8 PP-wave migration results of single shot(a) Reference; (b) Not compensated; (c) Compensated.

圖9 PS波單炮偏移(a) 參考; (b) 未補償; (c) 補償后.Fig.9 PS-wave migrated results of single shot(a) Reference; (b) Not compensated; (c) Compensated.

圖10 從圖8和圖9中抽出的單道波形(a) PP波; (b) PS波.Fig.10 Comparison of single-trace waveforms extracted from Figs.8 and 9(a) PP-waves; (b) PS-waves.

圖11 圖8和圖9對應(yīng)的振幅譜和功率譜 (a)和(c)為PP波; (b)和(d)為PS波.Fig.11 Amplitude spectra and power spectra of Figs.8 and 9(a) and (c) PP waves; (b) and (d) PS waves.

圖12 氣云模型(a) 速度模型; (b) Q模型.Fig.12 Gas chimney model(a) Velocity model; (b) Q model.

圖17為PP波和PS波40炮的偏移結(jié)果，其中圖17a、17b、17c分別為PP波參考、未補償和補償后的偏移剖面；圖17d、17e、17f分別為PS波參考、未補償和補償后的偏移剖面．從圖中可以看出，無論是縱波還是轉(zhuǎn)換波，新方法均對地震記錄的衰減能量進行了補償，對相位進行了校正，補償后的能量均強于沒補償?shù)钠破拭?，尤其在深部，并且補償后的同相軸更連續(xù)．

圖18為圖17對應(yīng)的振幅譜和功率譜，從頻譜對比可以看到，頻率成分(尤其是中高頻)能量得到了加強，即相對未補償?shù)钠破拭妫a償后的剖面分辨率更高．

5 結(jié)論

本文提出的基于吸收衰減補償?shù)亩喾至扛咚故B前逆時偏移可以有效補償衰減對多分量成像結(jié)果的影響．基于衰減補償?shù)亩喾至扛咚故鏁r偏移是一種優(yōu)秀的偏移算法，因為它不僅具有常規(guī)射線類Q偏移方法高效的優(yōu)點，而且具有接近于波動方程逆時Q偏移的成像精度．?dāng)?shù)值實驗中，對PP波和PS波的準(zhǔn)確成像表明了新方法能有效地補償多分量地震數(shù)據(jù)的衰減效應(yīng)，提高分辨率，改進成像效果，尤其對于強吸收區(qū)域．本文是在波場分離之后分別對PP和PS波做標(biāo)量偏移，處理過程中沒有考慮波傳播的矢量特征．因此，為了得到更好的轉(zhuǎn)換波成像結(jié)果，作者下一步將研究基于矢量波場延拓的黏彈性高斯束逆時偏移方法．

圖13 x分量：(a) 彈性; (b) 黏彈性Fig.13 x-component：(a) Elastic; (b) Viscoelastic

圖14 z分量：(a) 彈性；(b) 黏彈性Fig.14 z-component：(a) Elastic; (b) Viscoelastic

圖15 從圖13和圖14中抽取的第101道波形(橫向距離1000 m)(a) x分量；(b) z分量.Fig.15 Comparison of 101st-trace wave forms extracted from Figs.13 and 14 (at distance of 1000 m)(a) x-component; (b) z-component.

致謝 感謝朱鐵源博士提供氣云速度和Q值模型；第一作者感謝國家留學(xué)基金委在聯(lián)合培養(yǎng)期間給予的資助．感謝兩位審稿專家提出的寶貴修改意見和本文編輯細(xì)致而負(fù)責(zé)的工作，從而使得論文更加完善.

附錄A 射線追蹤

高斯束是波動方程沿中心射線的高頻近似解，解的過程包括運動學(xué)射線追蹤(計算射線路徑和走時)和動力學(xué)射線追蹤(計算中心射線周圍的能量值)．

對于光滑的二維速度模型(P波或S波)，運動學(xué)射線追蹤公式為：

圖16 圖13、14對應(yīng)的振幅譜和功率譜(a)和(c)分別為x分量; (b)和(d)為z分量.Fig.16 Amplitude spectra and power spectra of Figs.13 and 14(a) and (c) x component; (b) and (d) z component.

圖17 PP波和PS波偏移剖面PP波：(a) 參考; (b) 未補償; (c) 補償后; PS波：(d) 參考; (e) 未補償; (f) 補償后.Fig.17 PP-wave and PS-wave migration resultsPP waves：(a) Reference; (b) Not compensated; (c) Compensated; PS waves：(d) Reference; (e) Not compensated; (f) Compensated．

圖18 圖17對應(yīng)的振幅譜和功率譜PP波：(a) 振幅譜; (c) 功率譜； PS波：(b) 振幅譜; (d) 功率譜.Fig.18 Amplitude spectra and power spectra of Fig.17PP waves: (a) Amplitude spectra； (c) Power spectra; PS waves：(b) Amplitude spectra； (d) Power spectra.

(A1)

動力學(xué)射線追蹤可以表示為：

(A2)

(A3)

式中，p(s)和q(s)是復(fù)值動力學(xué)參數(shù)，它們決定了高斯束的波前曲率和束寬．

Abramowitz M，Stegun I A．1965. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover New York.

Bai J Y, Chen G Q, Yingst D, et al. 2013. Attenuation compensation in viscoacoustic reverse time migration.∥ 83rd Annual International Meeting, SEG, Expanded Abstracts, 3825-3830.

Bi L F, Qin N, Yang X D, et al. 2015. Gauss beam reverse time migration method for elastic multiple wave.GeophysicalProspectingforPetroleum(in Chinese), 54(1): 64-70.

Chang W F, McMechan G A. 1994. 3-D elastic prestack reverse-time depth migration.Geophysics, 59(4): 597-609.

Deng F, McMechan G A. 2007. True-amplitude prestack depth migration.Geophysics, 72(3): S155-S166.

Du Q Z, Hou B. 2008. Elastic Kirchhoff migration of vectorial wave-fields.AppliedGeophysics, 5(4): 284-293.

Du Q Z, Zhu Y T, Ba J. 2012. Polarity reversal correction for elastic reverse time migration.Geophysics, 77(2): S31-S41.

Du Q Z, Gong X F, Zhang M Q, et al. 2014. 3D PS-wave imaging with elastic reverse-time migration.Geophysics, 79(5): S173-S184. Fisk M D, McCartor G D. 1991. The phase screen method for vector elastic waves.JournalofGeophysicalResearch:SolidEarth(1978—2012), 96(B4): 5985-6010.

Fletcher R P, Nichols D, Cavalca M. 2012. Wavepath-consistent effectiveQestimation forQ-compensated reverse-time migration.∥ 74th EAGE Conference and Exhibition, Extended Abstracts. A020. Futterman W I. 1962. Dispersive body waves.JournalofGeophysicalResearch, 67(13): 5279-5291. Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration.Geophysics, 74(2): S11-S23.

Han J G, Wang Y, Lu J. 2013. Multi-component Gaussian beam prestack depth migration.JournalofGeophysicsandEngineering, 10(5): 055008. Hill N R. 1990. Gaussian beam migration.Geophysics, 55(11): 1416-1428. Hill N R. 2001. Prestack Gaussian-beam depth migration.Geophysics, 66(4): 1240-1250. Hokstad K. 2000. Multicomponent Kirchhoff migration.Geophysics, 65(3): 861-873. Hou A, Marfurt K J. 2002. Multicomponent prestack depth migration by scalar wavefield extrapolation.Geophysics, 67(6): 1886-1894.

Jin S W, Wu R S, Xie X B, et al. 1998. Wave equation-based decomposition and imaging for multicomponent seismic data.JournalofSeismicExploration, 7(2): 145-158. Jo C H, Shin C, Suh J H. 1996. An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator.Geophysics, 61(2): 529-537. Keers H, Vasco D W, Johnson L R. 2001. Viscoacoustic crosswell imaging using asymptotic waveforms.Geophysics, 66(3): 861-870. Kuo J T, Dai T F. 1984. Kirchhoff elastic wave migration for the case of noncoincident source and receiver.Geophysics, 49(8): 1223-1238.

Li X L, Mao W J. 2015. Multicomponent Gaussian beam migration in elastic medium.∥ 77th EAGE Conference and Exhibition, Extended Abstracts.

Li Z C, Guo Z B, Tian K. 2014. Least-squares reverse time migration in visco-acoustic medium.ChineseJ.Geophys. (in Chinese), 57(1): 214-228, doi: 10.6038/cjg20140118.

Li Z Y, Liang G H, Gu B L. 2013. Improved method of separating P- and S-waves using divergence and curl.ChineseJ.Geophys. (in Chinese), 56(6): 2012-2022, doi: 10.6038/cjg20130622.

Mittet R. 2007. A simple design procedure for depth extrapolation operators that compensate for absorption and dispersion.Geophysics, 72(2): S105-S112.

Popov M M, Semtchenok N M, Popov P M, et al. 2007. Reverse time migration with Gaussian beams and its application to a few synthetic data sets.∥ 76th EAGE Conference and Exhibition, Extended Abstracts, 2165-2169.

Popov M M, Semtchenok N M, Popov P M, et al. 2008. Reverse time migration with Gaussian beams and velocity analysis applications.∥ 70th EAGE Conference and Exhibition, Extended Abstracts, F048Popov M M, Semtchenok N M, Popov P M, et al. 2010. Depth migration by the Gaussian beam summation method.Geophysics, 75(2): S81-S93. Suh S, Yoon K, Cai J, et al. 2012. Compensating visco-acoustic effects in anisotropic reverse-time migration.∥ 82nd Annual International Meeting, SEG, Expanded Abstracts, 1-5.

Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data.Geophysics, 66(5): 1519-1527. Sun R, McMechan G A, Hsiao H H, et al. 2004. Separating P-and S-waves in prestack 3D elastic seismograms using divergence and curl.Geophysics, 69(1): 286-297.

Sun R, McMechan G A, Lee C S, et al. 2006. Prestack scalar reverse-time depth migration of 3D elastic seismic data.Geophysics, 71(5): S199-S207.

Traynin P, Liu J, Reilly J M. 2008. Amplitude and bandwidth recovery beneath gas zones using Kirchhoff prestack depthQ-migration.∥ 78th Annual International Meeting, SEG, Expanded Abstracts, 2412-2416. Valenciano A A, Chemingui N, Whitmore D, et al. 2011. Wave equation migration with attenuation and anisotropy compensation.∥ 81st Annual International Meeting, SEG, Expanded Abstracts, 232-236. Wang Y, Nemeth T. 1997. Multi-component separation of PP and SS by a least-squares migration method: Synthetic and field data tests.∥ 67th Annual International Meeting, SEG, Expanded Abstracts, 1222-1225.

Wang Y H, Guo J. 2004. Seismic migration with inverseQfiltering.GeophysicalResearchLetters, 31(21): L21608.

Wang Y H. 2008. Inverse-Qfiltered migration.Geophysics, 73(1):

S1-S6.

Wapenaar C P A, Herrmann P, Verschuur D J, et al. 1990. Decomposition of multicomponent seismic data into primary P- and S-wave responses.GeophysicalProspecting, 38(6): 633-661. Wu R S. 1994. Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method.JournalofGeophysicalResearch:SolidEarth(1978—2012), 99(B1): 751-766. Xiao X, Hao F, Egger C, et al. 2014. Final laser-beamQ-migration.∥ 84th Annual International Meeting, SEG, Expanded Abstracts, 3872-3876.

Xie X B, Wu R S. 2005. Multicomponent prestack depth migration using the elastic screen method.Geophysics, 70(1): S30-S37.

Xie Y, Xin K F, Sun J, et al. 2009. 3D prestack depth migration with compensation for frequency dependent absorption and dispersion.∥ 79th Annual International Meeting, SEG, Expanded Abstracts, 2919-2923.

Xie Y, Notfors C, Sun J, et al. 2010. 3D prestack beam migration with compensation for frequency dependent absorption and dispersion.∥ 72nd EAGE Conference and Exhibition, Extended Abstracts.

Yan H Y, Liu Y. 2013. Visco-acoustic prestack reverse-time migration based on the time-space domain adaptive high-order finite-difference method.GeophysicalProspecting, 61(5): 941-954.

Yan J, Sava P. 2008. Isotropic angle-domain elastic reverse-time migration.Geophysics, 73(6): S229-S239.

Zhang J F, Wapenaar K. 2002. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media.GeophysicalProspecting, 50(6): 629-643.

Zhang Y, Karazincir M, Notfors C, et al. 2002. Amplitude preserving v(z) pre-stack Kirchhoff migration, demigration and modeling.∥ 64th EAGE Conference and Exhibition, Extended Abstracts, B011.

Zhang Y, Zhang P, Zhang H Z. 2010. Compensating for visco-acoustic effects in reverse-time migration.∥ 80th Annual International Meeting, SEG, Expanded Abstracts, 3160-3164.

Zhu T Y, Harris J M, Biondi B. 2014.Q-compensated reverse-time migration.Geophysics, 79(3): S77-S87.

Zhu T, Carcione J M. 2014. Theory and modelling of constant-QP- and S-waves using fractional spatial derivatives.GeophysicalJournalInternational, 196(3): 1787-1795.

附中文參考文獻

畢麗飛, 秦寧, 楊曉東等. 2015. 彈性多波高斯束逆時偏移方法. 石油物探, 54(1): 64-70.

李振春, 郭振波, 田坤. 2014. 黏聲介質(zhì)最小平方逆時偏移. 地球物理學(xué)報, 57(1): 214-228, doi: 10.6038/cjg20140118.

李志遠(yuǎn), 梁光河, 谷丙洛. 2013. 基于散度和旋度縱橫波分離方法的改進. 地球物理學(xué)報, 56(6): 2012-2022, doi: 10.6038/cjg20130622.

(本文編輯 何燕)

Multiple-component Gaussian beam reverse-time migration based on attenuation compensation

BAI Min1,2,3， CHEN Xiao-Hong1,2， WU Juan1,2,3， CHEN Yang-Kang4， LIU Guo-Chang1， WANG En-Jiang1

1StateKeyLaboratoryofPetroleumResourcesandProspecting，ChinaUniversityofPetroleum，Beijing102249，China2NationalEngineeringLaboratoryforOffshoreOilExploration，ChinaUniversityofPetroleum，Beijing102249，China3SchoolofResourcesandEnvironment,NorthChinaUniversityofWaterResourcesandElectricPower,Zhengzhou450045,China4BureauofEconomicGeology，JohnA.andKatherineG.JacksonSchoolofGeosciences，TheUniversityofTexasatAustin，UniversityStation，BoxX，Austin，TX78713-8924，USA

Anelastic properties of subsurface media can cause amplitude loss and phase distortion of seismic waves，especially in high-attenuation areas such as the gas chimneys as observed in several oil and gas fields．In migration of such data sets，we usually obtain poor seismic images of the structure within and below high-attenuation gas-filled reservoirs．To improve the resolution of the migration image，we must deal with these attenuation effects．

Attenuation compensation； Multiple-component； Gaussian beam； Reverse-time migration； Green function

10.6038/cjg20160921.

海洋石油勘探國家工程實驗室“斜纜采集地震數(shù)據(jù)分析與處理技術(shù)研究”課題，國家自然科學(xué)基金項目(U1262207，41404099)，河南省重點科技攻關(guān)項目(152102210111)聯(lián)合資助．

白敏，男，1986年生，博士，講師，主要從事地震波傳播及偏移成像方面的研究.E-mail:bmwj103@163.com

10.6038/cjg20160921

P631

2015-06-02，2015-10-12收修定稿

白敏，陳小宏，吳娟等. 2016. 基于吸收衰減補償?shù)亩喾至扛咚故鏁r偏移. 地球物理學(xué)報，59(9):3379-3393,

Bai M， Chen X H， Wu J, et al. 2016. Multiple-component Gaussian beam reverse-time migration based on attenuation compensation.ChineseJ.Geophys. (in Chinese),59(9):3379-3393,doi:10.6038/cjg20160921.