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

白敏， 陳小宏， 吳娟， 陳陽(yáng)康， 劉國(guó)昌， 王恩江

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 2 中國(guó)石油大學(xué)(北京)海洋石油勘探國(guó)家工程實(shí)驗(yàn)室， 北京 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

?

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

白敏1,2,3， 陳小宏1,2， 吳娟1,2,3， 陳陽(yáng)康4， 劉國(guó)昌1， 王恩江1

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

衰減補(bǔ)償； 多分量； 高斯束； 逆時(shí)偏移； 格林函數(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)注，這就要求在偏移延拓和成像中對(duì)影響成像振幅的因素做校正處理(Zhang et al．，2002)．然而，固有的吸收因子在常規(guī)偏移方法中經(jīng)常被忽略．因?yàn)樯顚痈哳l損失，地震波穿過(guò)覆蓋層時(shí)，介質(zhì)的空間變化引起地震振幅衰減，子波相位畸變，成像分辨率降低，使得異常區(qū)(比如氣層)振幅衰減，頻帶變窄，造成深層識(shí)別和解釋的困難，也影響了準(zhǔn)確預(yù)測(cè)儲(chǔ)層的能力．因此，需要補(bǔ)償因?yàn)檫@類衰減引起的吸收效應(yīng)．

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

多分量地震勘探能得到PP波和PS波數(shù)據(jù)，與單純的PP波相比，PS數(shù)據(jù)攜帶了目標(biāo)區(qū)域不同的信息，能提供更多的物性參數(shù)，有效改進(jìn)地震勘探的準(zhǔn)確性．比如在PP波信號(hào)弱的區(qū)域，PS波可以增強(qiáng)照明．聯(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)鍵是計(jì)算PP波和PS波的走時(shí)，以此求和計(jì)算成像振幅．為了克服Kirchhoff成像的弱點(diǎn)，一些學(xué)者提出了彈性單程波動(dòng)方程偏移(Fisk and McCartor，1991；Wu，1994；Xie and Wu，2005)．由于是雙程波的近似，單程波偏移方法受到成像角度、波的傳播路徑問(wèn)題的限制，影響成像精度．因此，近年來(lái)利用雙程波動(dòng)方程的多分量逆時(shí)偏移方法成為熱點(diǎn)．多分量的逆時(shí)偏移方法有兩類：彈性波偏移和標(biāo)量偏移(Hou and Marfurt，2002)．彈性波偏移是直接輸入多分量地震數(shù)據(jù)，利用彈性波動(dòng)方程構(gòu)建矢量波場(chǎng)(Chang and McMechan，1994；Yan and Sava，2008；Du et al．，2012，2014)．此外，也有些研究者通過(guò)把波場(chǎng)分離成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ì)體．波動(dòng)方程逆時(shí)偏移更精確，但是計(jì)算量大．Popov等(2007，2008，2010)結(jié)合高斯束的高效靈活性和波動(dòng)方程逆時(shí)偏移的高精度，提出了高斯束逆時(shí)偏移，即以Kirchhoff積分為基礎(chǔ)，利用高斯束的疊加積分計(jì)算格林函數(shù)的逆時(shí)偏移算法．由于高斯束的疊加積分是嚴(yán)格按照數(shù)學(xué)理論推導(dǎo)的，因此避免了多值走時(shí)和不精確的振幅計(jì)算，另外可通過(guò)控制成像條件的時(shí)間窗來(lái)消除續(xù)至波，更適合基于目標(biāo)的成像．與Hill(2001)提出的高斯束疊前深度偏移方法相比，高斯束逆時(shí)偏移不做局部?jī)A斜疊加和最速下降近似，在理論上要優(yōu)于前者．

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

2 多分量高斯束逆時(shí)偏移

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

(1)

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

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

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

(2)

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

×GGB(x,xs;ω)，

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

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

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

(4)

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

(5)

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

圖1 偏移域Ω的參數(shù)化xs和x0分別表示震源和地下成像區(qū)域(虛線內(nèi))任一點(diǎn)，v(x)是速度，nx是Ω在記錄表面的外法向?qū)?shù)，Ω的邊界向外延伸至無(wú)窮.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è)對(duì)波場(chǎng)的貢獻(xiàn)可以忽略．格林函數(shù)滿足Kirchhoff近似中的邊界條件G|z=0=0，則偏移域內(nèi)一點(diǎn)x0的反向延拓波場(chǎng)為：

(6)

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

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

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

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

(7)

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

(8)

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

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

×GGB(x,x0,ω)，

(9)

其中，

(10)

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

對(duì)于P波，有

(11)

對(duì)于S波，有

(12)

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

(13)

(14)

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

(15)

(16)

2.3 成像條件

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

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á)式見(jiàn)方程(3)和(6)．

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

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

3.1 衰減和補(bǔ)償原理

圖3 無(wú)衰減介質(zhì)與衰減介質(zhì)正演和偏移示意圖. (a)、(c)為無(wú)衰減介質(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ù)速度的實(shí)部分別是彈性介質(zhì)中的速度vP和vS，品質(zhì)因子Q代表衰減．QP和QS分別是縱橫波的品質(zhì)因子，如果衰減小(1/Q?1)，則復(fù)速度為(Keers et al．，2001)：

(19)

(20)

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

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

(21)

(22)

其中，

(23)

(24)

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

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

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

(25)

其中，復(fù)速度

(26)

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

(27)

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

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

圖4 無(wú)衰減介質(zhì)中，高斯束積分和波動(dòng)方程計(jì)算的格林函數(shù)(a)及其與解析格林函數(shù)的精度對(duì)比(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ì)中，高斯束積分和波動(dòng)方程計(jì)算的格林函數(shù)及其精度對(duì)比(局部放大圖)(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ù)值算例

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

4.1 兩層模型

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

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

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

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

4.2 氣云模型

在第二個(gè)例子中，我們考慮實(shí)際的衰減模型．圖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)接收．接收記錄長(zhǎng)度為3 s，時(shí)間采樣間隔為0.001 s．

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

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

圖7 波到達(dá)界面之后某一時(shí)刻的原始波場(chǎng)快照和補(bǔ)償后的波場(chǎng)快照震源波場(chǎng)：(a) 參考，(b) 補(bǔ)償后；地下一點(diǎn) (1500 m,1500 m)的PP波波場(chǎng)：(c)參考，(d)補(bǔ)償后；PS波波場(chǎng)：(e)參考，(f)補(bǔ)償后.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) 未補(bǔ)償; (c) 補(bǔ)償后.Fig.8 PP-wave migration results of single shot(a) Reference; (b) Not compensated; (c) Compensated.

圖9 PS波單炮偏移(a) 參考; (b) 未補(bǔ)償; (c) 補(bǔ)償后.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對(duì)應(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波參考、未補(bǔ)償和補(bǔ)償后的偏移剖面；圖17d、17e、17f分別為PS波參考、未補(bǔ)償和補(bǔ)償后的偏移剖面．從圖中可以看出，無(wú)論是縱波還是轉(zhuǎn)換波，新方法均對(duì)地震記錄的衰減能量進(jìn)行了補(bǔ)償，對(duì)相位進(jìn)行了校正，補(bǔ)償后的能量均強(qiáng)于沒(méi)補(bǔ)償?shù)钠破拭妫绕湓谏畈?，并且補(bǔ)償后的同相軸更連續(xù)．

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

5 結(jié)論

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

圖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值模型；第一作者感謝國(guó)家留學(xué)基金委在聯(lián)合培養(yǎng)期間給予的資助．感謝兩位審稿專家提出的寶貴修改意見(jiàn)和本文編輯細(xì)致而負(fù)責(zé)的工作，從而使得論文更加完善.

附錄A 射線追蹤

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

對(duì)于光滑的二維速度模型(P波或S波)，運(yùn)動(dòng)學(xué)射線追蹤公式為：

圖16 圖13、14對(duì)應(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) 未補(bǔ)償; (c) 補(bǔ)償后; PS波：(d) 參考; (e) 未補(bǔ)償; (f) 補(bǔ)償后.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對(duì)應(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)

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

(A2)

(A3)

式中，p(s)和q(s)是復(fù)值動(dòng)力學(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.

附中文參考文獻(xiàn)

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

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

李志遠(yuǎn), 梁光河, 谷丙洛. 2013. 基于散度和旋度縱橫波分離方法的改進(jìn). 地球物理學(xué)報(bào), 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.

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

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

10.6038/cjg20160921

P631

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

白敏，陳小宏，吳娟等. 2016. 基于吸收衰減補(bǔ)償?shù)亩喾至扛咚故鏁r(shí)偏移. 地球物理學(xué)報(bào)，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.