黏滯聲波高斯束疊前深度偏移

吳 娟，陳小宏，白 敏

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

?

黏滯聲波高斯束疊前深度偏移

吳 娟1,2，陳小宏1,2，白 敏1,2

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

高斯束偏移不僅具有接近于波動(dòng)方程偏移的成像精度，而且保留了Kirchhoff 積分法高效、靈活的優(yōu)點(diǎn)，可以對(duì)復(fù)雜介質(zhì)準(zhǔn)確成像。由于實(shí)際地下介質(zhì)具有黏滯性，因此研究黏滯聲波疊前深度偏移具有一定的現(xiàn)實(shí)意義。筆者采用高斯束偏移方法對(duì)地震數(shù)據(jù)進(jìn)行吸收衰減補(bǔ)償。首先給出共炮域高斯束疊前深度偏移原理；然后在此基礎(chǔ)上推導(dǎo)補(bǔ)償吸收衰減的表達(dá)式，校正品質(zhì)因子Q引起的振幅衰減和相位畸變，實(shí)現(xiàn)基于吸收衰減補(bǔ)償?shù)母咚故B前深度偏移；最后用兩層模型和氣云模型對(duì)偏移方法進(jìn)行了測(cè)試。結(jié)果表明，在考慮地下介質(zhì)的黏滯性時(shí)，黏滯聲波高斯束疊前深度偏移比聲波高斯束疊前深度偏移具有更高的成像分辨率。

高斯束；衰減補(bǔ)償；黏滯聲波；格林函數(shù)；疊前深度偏移

0 引言

地下介質(zhì)的黏滯性造成地震波的振幅衰減和相位畸變，尤其是在氣云這樣的高衰減區(qū)。對(duì)這類數(shù)據(jù)偏移時(shí)，氣云和它下方的構(gòu)造不能很好地成像[1]。為了提高偏移成像分辨率，需要補(bǔ)償這類衰減效應(yīng)。

早期的吸收補(bǔ)償是對(duì)地震道做反Q(品質(zhì)因子)濾波[2-4]，但是其基于一維波動(dòng)方程對(duì)地震記錄逐道補(bǔ)償，不能適應(yīng)復(fù)雜構(gòu)造。實(shí)際上，Q值是空間變化的，吸收衰減對(duì)地震信號(hào)的影響與波的傳播路徑有關(guān)[5]。所以，理想的補(bǔ)償方法應(yīng)該根據(jù)其實(shí)際傳播路徑，在偏移成像的過(guò)程中進(jìn)行[6]。過(guò)去的二十幾年中，有許多學(xué)者利用頻率域的單程波方程補(bǔ)償上、下行波的衰減和頻散效應(yīng)[7-10]。概括來(lái)說(shuō)，這些方法通常是把衰減表示為復(fù)相速度的虛部，因此可以靈活處理任意衰減模型。

基于波動(dòng)方程偏移的衰減補(bǔ)償能對(duì)復(fù)雜構(gòu)造進(jìn)行更準(zhǔn)確的成像，但是計(jì)算效率較低。Kirchhoff偏移在補(bǔ)償衰減的過(guò)程中不需要輸出每一個(gè)CIG(共成像點(diǎn))道集，計(jì)算速度快，不足的地方是存在多值走時(shí)問(wèn)題。高斯束偏移[11-17]是近年來(lái)發(fā)展的一種優(yōu)秀的偏移算法，不僅具有接近于波動(dòng)方程偏移的成像精度，而且保留了Kirchhoff 積分法高效靈活的優(yōu)點(diǎn)，能夠有效處理多初至。

本文在高斯束偏移的基礎(chǔ)上，研究基于吸收衰減補(bǔ)償?shù)母咚故B前深度偏移。首先給出共炮域高斯束疊前深度偏移原理；然后推導(dǎo)補(bǔ)償吸收衰減的表達(dá)式，校正Q引起的振幅衰減和相位畸變，實(shí)現(xiàn)衰減補(bǔ)償?shù)母咚故B前深度偏移；最后用數(shù)值算例證明該方法的有效性和優(yōu)越性。

1 方法原理

高斯束方法的核心是格林函數(shù)的合成。利用格林函數(shù)可以得到地震波場(chǎng)，描述地震波的傳播過(guò)程。因此，這里首先介紹格林函數(shù)的合成方法，然后直接給出由格林函數(shù)表示的共炮道集高斯束疊前深度偏移公式。

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

在高斯束方法中，二維格林函數(shù)G2D(x,xs,ω)是通過(guò)一系列由源點(diǎn)出射的具有不同出射角θ的高斯束疊加積分表示的[18]。其中：ω是圓頻率；x=(x,z)和xs=(xs,zs)是地下任意一點(diǎn)和震源點(diǎn)的坐標(biāo)，如圖1所示。

圖1 格林函數(shù)的高斯束疊加合成Fig.1 Green function in terms of Gaussian beams summation

從源點(diǎn)xs處以不同的角度出射中心射線束，每條中心束附近波場(chǎng)值用高斯束方法求取。地下介質(zhì)中x點(diǎn)的波場(chǎng)值由與其臨近的多條高斯束疊加而得，則高斯束表示的格林函數(shù)公式為[12,16]

(1)

式中：ps=(px,pz)為中心射線的參數(shù)矢量，px和pz分別表示射線參數(shù)的水平分量和垂直分量；uGB(x,xs,ps,ω)為高斯束方法求取的波場(chǎng)位移,

(2)

式中，A,tc分別是高斯束的復(fù)值振幅和復(fù)值旅行時(shí)，且

(3)

(4)

式中：s0是射線追蹤的起點(diǎn)；s是射線路徑；v(s)是沿射線的速度；τ(s)是沿著射線的旅行時(shí)；p(s)和q(s)是動(dòng)力學(xué)射線追蹤方程組的復(fù)值解[18]，它們決定了高斯束的能量分布；n表示沿射線法線方向的距離。

高斯束構(gòu)建的格林函數(shù)由多條中心射線附近有限區(qū)域的局部波場(chǎng)疊加得到，利用不同的波束疊加可以解決多值走時(shí)問(wèn)題。采用Hill[11]給定的初始值求解動(dòng)力學(xué)射線追蹤方程組，可以保證高斯束是正則的，因而由高斯束所表示的格林函數(shù)也是處處正則的，避免了復(fù)雜構(gòu)造下的焦散問(wèn)題。

1.2 共炮集高斯束疊前深度偏移

高斯束疊前偏移由上行波場(chǎng)和下行波場(chǎng)的互相關(guān)得到，原理如圖2所示。在震源和束中心處分別以不同的射線參數(shù)ps和pLr出射高斯束進(jìn)行波場(chǎng)計(jì)算。其中：s表示震源位置；Lr表示束中心位置。

圖2 共炮道集高斯束疊前深度偏移原理圖Fig.2 Sketch of common shot Gaussian beam prestack depth migration

利用高斯束積分表征的格林函數(shù)來(lái)描述上、下行波場(chǎng)，則共炮道集高斯束疊前成像公式[19]為

(5)

高斯函數(shù)具有如下性質(zhì)：

(6)

式中：ωr為參考頻率；ΔL為束中心間隔；L0表示ωr下的初始束寬[11,13]。將式(6)代入式(5)，得到

(7)

式中：Ics(x)表示x點(diǎn)的共炮偏移成像值；震源格林函數(shù)G(x,xs,ω)用式(1)表示；束中心Lr有效范圍內(nèi)接收點(diǎn)的格林函數(shù)G(x,xr,ω)用下式表示：

(8)

將這兩個(gè)格林函數(shù)分別代入式(7)可得

(9)

式(9)進(jìn)一步簡(jiǎn)化為

(10)

式(10)即為最終的疊前共炮點(diǎn)道集高斯束偏移公式。其中：pLr=(pLrx,pLrz)是接收點(diǎn)射線參數(shù)；Ds(Lr,pLr,ω)為地震記錄的加窗局部?jī)A斜疊加。

(11)

1.3 衰減介質(zhì)共炮集高斯束疊前深度偏移

圖3為衰減介質(zhì)中的正演與偏移示意圖。圖3a中，由于品質(zhì)因子Q的存在，波從震源到反射點(diǎn)、再到檢波點(diǎn)，沿整個(gè)傳播路徑的累積衰減量為exp(-αLU)exp(-αLD)(α是衰減系數(shù)，LU和LD分別為上行波和下行波的傳播距離)，則衰減后的檢波點(diǎn)波場(chǎng)為

(12)

為了實(shí)現(xiàn)對(duì)檢波點(diǎn)波場(chǎng)的完全補(bǔ)償，檢波點(diǎn)的校正量應(yīng)該為exp(αLU)exp(αLD)。然而在互相關(guān)成像中，傳播距離為L(zhǎng)U(檢波點(diǎn)到反射點(diǎn))的檢波點(diǎn)補(bǔ)償算子為exp(αLU)，補(bǔ)償后的檢波點(diǎn)波場(chǎng)為

(13)

因此，還需要同時(shí)補(bǔ)償下行(震源)波場(chǎng)。沿著射線路徑、傳播距離為L(zhǎng)D(震源到反射點(diǎn))的震源補(bǔ)償算子為exp(αLD)，補(bǔ)償后的震源波場(chǎng)為

(14)

a. 正演；b. 偏移。S(x,t)為參考震源波場(chǎng)；R(x,t)為參考檢波點(diǎn)波場(chǎng)(即無(wú)衰減地震波場(chǎng))。圖3 衰減介質(zhì)中的正演與偏移示意圖Fig.3 Schematic of forward modeling and migration in an attenuating medium

如圖3b所示。這樣補(bǔ)償后的互相關(guān)成像剖面才能達(dá)到與參考偏移剖面相當(dāng)?shù)男Ч?/p>

波在黏聲波介質(zhì)中傳播可認(rèn)為是復(fù)速度的波在聲波介質(zhì)中傳播，復(fù)速度的實(shí)部是聲波介質(zhì)中的速度 ，品質(zhì)因子Q代表衰減。如果衰減小(1/Q≤1)，則復(fù)速度可以表示為[20-21]

(15)

其中：Q不依賴于頻率；速度的虛部引起沿射線路徑的指數(shù)衰減；速度的實(shí)部包含頻散項(xiàng)，確保波動(dòng)方程解的因果性。當(dāng)1/Q≤1時(shí)，對(duì)于一階項(xiàng)，射線路徑保持不變。因此，衰減只通過(guò)依賴于頻率的復(fù)值旅行時(shí)影響波形。復(fù)值旅行時(shí)為

(16)

2 模型算例分析

2.1 兩層模型

筆者設(shè)計(jì)了簡(jiǎn)單的兩層模型。模型大小為3 000 m×3 000 m，如圖4所示。震源采用主頻為20 Hz的Ricker子波，坐標(biāo)為(1 500 m，0 m)，共設(shè)計(jì)了301個(gè)檢波點(diǎn)，布置在地表0～ 3 000 m處，道間距為10 m，接收記錄長(zhǎng)度為2 s，時(shí)間采樣間隔是0.001 s。

圖4 兩層介質(zhì)速度和Q值模型Fig.4 Two layer model with velocity and Q values

a. 震源波場(chǎng)；b. 檢波點(diǎn)波場(chǎng)；c. 互相關(guān)結(jié)果。圖5 參考的震源、檢波點(diǎn)波場(chǎng)以及互相關(guān)成像剖面Fig.5 Reference source wavefield, receiver wavefield and cross-correlation results

a. 震源波場(chǎng)；b. 檢波點(diǎn)波場(chǎng)；c. 互相關(guān)結(jié)果。圖6 未補(bǔ)償?shù)恼鹪?、檢波點(diǎn)波場(chǎng)以及互相關(guān)成像剖面Fig.6 Noncompensated source wavefield, receiver wavefield and cross-correlation results

a. 震源波場(chǎng)；b. 檢波點(diǎn)波場(chǎng)；c. 互相關(guān)結(jié)果。圖7 補(bǔ)償后的震源、檢波點(diǎn)波場(chǎng)以及互相關(guān)成像剖面Fig.7 Compensated source wavefield, receiver wavefield and cross-correlation results

為了更直觀地反映衰減和補(bǔ)償過(guò)程，將分別輸出界面附近某一時(shí)刻的參考、未補(bǔ)償和補(bǔ)償后的震源、檢波點(diǎn)波場(chǎng)快照以及互相關(guān)成像剖面。震源波場(chǎng)由點(diǎn)源出射的高斯束正向延拓得到；檢波點(diǎn)波場(chǎng)由檢波點(diǎn)處束中心局部?jī)A斜疊加后的平面波反向外推得到；互相關(guān)成像剖面為單炮偏移結(jié)果。如圖5，6和7所示，炮點(diǎn)位于(1 500 m，0 m)，圖中虛線為界面位置(深度1 500 m)。

圖5為參考波場(chǎng)，沒(méi)有考慮吸收衰減。圖6是未補(bǔ)償?shù)牟▓?chǎng)：圖6a與圖5a相同；6b是檢波點(diǎn)波場(chǎng)，其輸入為黏聲波地震記錄，與參考波場(chǎng)圖5b相比，存在明顯的振幅衰減和相位畸變(能量較弱，同一時(shí)刻的波場(chǎng)位置下移)。圖7是補(bǔ)償后的結(jié)果：與未補(bǔ)償?shù)?圖6)相比，補(bǔ)償后的震源(圖7a)和檢波點(diǎn)波場(chǎng)(圖7b)振幅均得到加強(qiáng)，相位得到校正(位置上移，圖7b、c恢復(fù)到圖5b、c的位置)。對(duì)比互相關(guān)結(jié)果(圖5c，圖6c，圖7c)可以看出，本文方法在偏移過(guò)程中有效地補(bǔ)償了地震記錄的能量衰減，補(bǔ)償后的成像效果(圖7c)與參考結(jié)果(圖5c)相當(dāng)。圖8是從圖5c，6c和7c中抽取的單道記錄對(duì)比結(jié)果。從圖8中也可以看出，黏聲波高斯束疊前深度偏移有效地補(bǔ)償了地下介質(zhì)對(duì)地震波的吸收。

2.2 氣云模型

在第二個(gè)例子中，我們考慮實(shí)際的衰減模型。圖9是氣云速度模型和相應(yīng)的Q模型(BP模型的一部分[22])。模型的頂部偏中是由于氣云的存在導(dǎo)致的低速和高衰減區(qū)。模型網(wǎng)格數(shù)為398×161，網(wǎng)格大小是10 m×10 m。共40炮，震源采用主頻為15 Hz的Ricker子波，炮間距為100 m，每一炮都是全地表等距(10 m)接收。接收記錄長(zhǎng)度為3 s，時(shí)間采樣間隔為0.001 s。

采用高階交錯(cuò)網(wǎng)格有限差分方法合成聲波和黏聲波地震記錄。如圖10所示，選取一炮的正演結(jié)果，炮點(diǎn)位于1 000 m處。從圖10中可以看出，黏聲波記錄反射波的同相軸能量較弱，且同相軸的波形較寬；表明介質(zhì)的黏滯性對(duì)地震波能量的吸收和衰減作用明顯。

圖11是對(duì)40炮記錄所做的偏移結(jié)果。從圖11中可以看出，補(bǔ)償后的能量強(qiáng)于沒(méi)補(bǔ)償?shù)钠破拭?，接近參考偏移剖面，尤其在深部，并且補(bǔ)償后的同相軸更連續(xù)。

圖12是分別從圖11a，b，c的剖面上抽出的第100道(橫向距離1 000 m處)和第200道(橫向距離2 000 m處)的單道記錄對(duì)比結(jié)果，由圖可知，補(bǔ)償后的振幅增強(qiáng)，相位得到校正。

圖13為圖12對(duì)應(yīng)的振幅譜，從頻譜對(duì)比可以看到，黏滯聲波高斯束疊前深度偏移有效地補(bǔ)償了地震波的能量，尤其是中高頻能量成分得到了加強(qiáng)，即相對(duì)聲波疊前深度偏移結(jié)果，黏滯聲波疊前深度偏移的分辨率要更高些。

圖8 從成像剖面(圖5c，6c和7c)中抽出的單道波形對(duì)比(水平距離1 500 m處)Fig.8 Comparison of the single trace waveforms (at lateral distance 1 500 m) extracted from image results (Figure 5c, 6c and 7c)

圖9 氣云速度和Q模型Fig.9 Gas chimney model with velocity and Q values

a. 聲波；b. 黏聲波。圖10 單炮記錄(炮點(diǎn)位于水平距離1 000 m處)Fig.10 Synthetic data of one shot (at a distance of 1 000 m)

a. 參考(聲波正演，聲波偏移)；b. 未補(bǔ)償(黏聲波正演，聲波偏移)；c. 補(bǔ)償后(黏聲波正演，黏聲波偏移)。圖11 氣云模型偏移疊加剖面Fig.11 Migration results of Gas chimney model

a. 第100道；b. 第200道。圖12 偏移剖面的單道波形對(duì)比Fig.12 Comparison of the single trace waveforms extracted from image results

a. 第100道；b. 第200道。圖13 偏移剖面的振幅譜對(duì)比Fig.13 Comparison of amplitude spectrum from image results

3 結(jié)語(yǔ)

高斯束偏移是一種優(yōu)秀的偏移算法，基于高斯束偏移方法進(jìn)行吸收衰減補(bǔ)償，不僅具有Kirchhoff 偏移方法靈活高效的優(yōu)點(diǎn)，而且具有接近于波動(dòng)方程偏移的成像精度。筆者提出的黏滯聲波高斯束疊前深度偏移可以補(bǔ)償衰減對(duì)成像結(jié)果的影響。實(shí)算表明，新方法能有效提高分辨率，獲得更好的成像剖面，尤其對(duì)于強(qiáng)吸收區(qū)域。

[1] Zhou J,Birdus S,Hung B,et al.Compensating Attenuation Due to Shallow Gas ThroughQTomography andQ-PSDM: A Case Study in Brazil[C]//81st Annual International Meeting. San Antonio: SEG, 2011: 3332-3336.

[2] Bickel S H, Natarajan R R. Plane-WaveQDeconvolution[J]. Geophysics, 1985, 50(9): 1426-1439.

[3] Hargreaves N D, Calvert A J. InverseQFiltering by Fourier Transform[J]. Geophysics, 1991, 56(4): 519-527.

[4] Wang Y H.InverseQFilter for Seismic Resolution Enhancement[J]. Geophysics, 2006, 71(3): 51-60.

[5] 馮晅, 張瑾, 劉財(cái), 等. 基于改進(jìn)的 Kolsky 模型波場(chǎng)延拓公式的縱波Q值, 橫波Q值估計(jì)[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版, 2014, 44(1): 359-368. Feng Xuan, Zhang Jin, Liu Cai, et al. Estimation of P-and S-Waves Quality Factors Based on the Formula of the Wave-Field Continuation in Modified Kolsky Model[J]. Journal of Jilin University: Earth Science Edition, 2014, 44(1): 359-368.

[6] Zhang C H, Ulrych T J. Refocusing Migrated Seismic Images in Absorptive Media[J]. Geophysics, 2010, 75(3): 103-110.

[7] Mittet R, Sollie R, Hokstad K. Prestack Depth Migration with Compensation for Absorption and Dispersion[J]. Geophysics, 1995, 60 (5): 1485-1494.

[8] Zhang J, Wapenaar K. Wavefield Extrapolation and Prestack Depth Migration in Anelastic Inhomogeneous Media[J]. Geophysical Prospecting, 2002, 50(10): 629-643.

[9] Mittet R. A Simple Design Procedure for Depth Extrapolation Operators that Compensate for Absorption and Dispersion[J]. Geophysics, 2007, 72(2): 105-112.

[10] Zhang J, Wu J, Li X. Compensation for Absorption and Dispersion in Prestack Migration: An EffectiveQApproach[J]. Geophysics, 2013, 78(1): 1-14.

[11] Hill N R. Gaussian Beam Migration[J]. Geophysics, 1990, 55: 1416-1428.

[12] Hill N R. Prestack Gaussian-Beam Depth Migration[J]. Geophysics, 2001, 66: 1240-1250.

[13] Hale D. Migration by the Kirchhoff, Slant Stack and Gaussian Beam Methods[R]. Golden: Colorado School of Mines Center for Wave Phenomena, 1992.

[14] Hale D. Computational Aspects of Gaussian Beam Migration[R]. Golden: Colorado School of Mines Center for Wave Phenomena, 1992.

[15] Gray S H. Gaussian Beam Migration of Common-Shot Records[J]. Geophysics, 2005, 70(4): 71-77.

[16] Gray S H, Bleistein N. True-Amplitude Gaussian-Beam Migration[J]. Geophysics, 2009, 74(2): 11-23.

[17] Popov M M, Semtchenok N M, Popov P M, et al. Depth Migration by the Gaussian Beam Summation Method[J]. Geophysics, 2010, 75(2): 81-93.

[18] Cerveny V. A Note on Dynamic Ray Tracing in Ray-Centered Coordinates in Anisotropic Inhomogeneous Media[J]. Studia Geophysica et Geodaetica, 2007, 51: 411-22.

[19] 岳玉波. 復(fù)雜介質(zhì)高斯束偏移成像方法研究[D]. 青島: 中國(guó)石油大學(xué), 2011. Yue Yubo. Study on Gaussian Beam Migration Methods in Complex Medium[D]. Qingdao: China University of Petroleum, 2011.

[20] Xie Y, Notfors C, Sun J, et al. 3D Prestack Beam Migration with Compensation for Frequency Depen-dent Absorption and Dispersion[C]//72nd EAGE Conference & Exhibition Incorporating SPE EUROPEC. Barcelona: [s. n.], 2010.

[21] Keers H, Vasco D W, Johnson L R. Viscoacoustic Crosswell Imaging Using Asymptotic Waveforms[J]. Geophysics, 2001, 66: 861-870.

[22] Zhu T Y, Harris J M, Biondi B.Q-Compensated Reverse-Time Migration[J]. Geophysics, 2014, 79(3): 77-87.

Viscoacoustic Gaussian Beam Prestack Depth Migration

Wu Juan1, 2，Chen Xiaohong1, 2，Bai Min1, 2

1.StateKeyLaboratoryofPetroleumResourcesandProspecting，ChinaUniversityofPetroleum，Beijing102249，China2.NationalEngineeringLaboratoryforOffshoreOilExploration，ChinaUniversityofPetroleum，Beijing102249，China

Gaussian beam migration is an effective and efficient depth migration method with its accuracy comparable to the wave equation migration and its flexibility comparable to Kirchhoff migration. The viscoacoustic prestack depth migration is of practical significance because it considers the viscosity of the subsurface media. Gaussian beam migration is used to perform seismic data compensation for frequency dependent absorption and dispersion. Gaussian beam prestack depth migration algorithms that operate on common-shot gathers are presented to derive expressions about attenuation and compensation which correct amplitude attenuation and phase distortion caused byQin order to realize Gaussian beam prestack depth migration. The numerical modeling suggests that the viscoacoustic Gaussian beam prestack depth migration have a higher imaging resolution than the acoustic Gaussian beam prestack depth when the viscosity of the subsurface is considered.

Gaussian beam；attenuation and compensation；viscoacoustic；Green function；prestack depth migration

10.13278/j.cnki.jjuese.201505302.

2014-12-22

國(guó)家自然科學(xué)基金項(xiàng)目(U1262207)；國(guó)家科技重大專項(xiàng)課題(2011ZX05023-005-005，2011ZX05019-006)

吳娟(1984--)，女，博士研究生，主要從事地震正演與偏移方法研究，E-mail:wujuan-0909@163.com。

10.13278/j.cnki.jjuese.201505302

P631.4

A

吳娟，陳小宏，白敏. 黏滯聲波高斯束疊前深度偏移.吉林大學(xué)學(xué)報(bào):地球科學(xué)版,2015,45(5):1530-1538.

Wu Juan, Chen Xiaohong, Bai Min. Viscoacoustic Gaussian Beam Prestack Depth Migration.Journal of Jilin University:Earth Science Edition,2015,45(5):1530-1538.doi:10.13278/j.cnki.jjuese.201505302.