基于反射波動方程的疊前地震反射數(shù)據(jù)波阻抗相對變化成像研究

陳生昌,周華敏

(1.浙江大學(xué)地球科學(xué)學(xué)院,浙江杭州310027;2.長江水利委員會長江科學(xué)院,湖北武漢430010)

基于反射波動方程的疊前地震反射數(shù)據(jù)波阻抗相對變化成像研究

陳生昌1,周華敏2

(1.浙江大學(xué)地球科學(xué)學(xué)院,浙江杭州310027;2.長江水利委員會長江科學(xué)院,湖北武漢430010)

波阻抗是地震數(shù)據(jù)巖性處理解釋中的重要參數(shù)。提出了一種基于反射波動方程的疊前地震反射數(shù)據(jù)深度域波阻抗相對變化成像方法。從變密度的聲波方程出發(fā),首先利用高頻近似,推導(dǎo)出基于波阻抗相對變化的一次反射波場的線性傳播方程,然后基于該傳播方程,利用線性反演理論推導(dǎo)出獲得深度域波阻抗相對變化近似估計的成像公式。受反射地震數(shù)據(jù)頻帶范圍的限制和波形線性反演問題近似求解成像方法的限制,該波阻抗相對變化成像方法所獲得的成像結(jié)果的保真性和分辨率會存在不足,但其計算量與逆時偏移的計算量基本相當(dāng)。合成數(shù)據(jù)的模型試驗結(jié)果驗證了波阻抗相對變化成像方法的有效性。

疊前地震反射數(shù)據(jù);高頻近似;反射波動方程;線性反演;波阻抗相對變化;成像

波阻抗是聯(lián)系地質(zhì)與反射地震數(shù)據(jù)的一個重要的物性參數(shù),在油氣地震勘探開發(fā)的數(shù)據(jù)處理解釋中,特別是在巖性油氣藏勘探開發(fā)中,發(fā)揮了十分重要的作用[1-2]。因此,波阻抗反演在地震數(shù)據(jù)處理中具有特殊的地位。

自20世紀(jì)70年代初LAVERGNE[3]首先提出波阻抗反演方法以來,有關(guān)波阻抗反演方法技術(shù)的研究有了很大的發(fā)展,提出了各種不同的波阻抗反演方法。如疊前地震數(shù)據(jù)波阻抗反演和疊后地震數(shù)據(jù)波阻抗反演。在疊前地震數(shù)據(jù)波阻抗反演方法中,有基于Zoeppritz方程近似的彈性波阻抗反演方法[4-5]和全波形反演方法[6-8]。在疊后地震數(shù)據(jù)波阻抗反演方法中,有直接反演方法[9]、基于模型反演方法[10-13]和非線性反演方法[14-15]。疊前地震數(shù)據(jù)波阻抗反演盡管計算復(fù)雜、計算量大,而且對數(shù)據(jù)的要求也高,但可以得到相對多的地層信息,特別是基于全波形的波阻抗反演,不僅可以充分利用地震數(shù)據(jù)的波形信息,還可以適應(yīng)地下介質(zhì)的橫向變化。因此疊前地震數(shù)據(jù)波阻抗反演是波阻抗反演的發(fā)展趨勢。疊后地震數(shù)據(jù)波阻抗反演是當(dāng)前生產(chǎn)中比較常用的方法技術(shù),計算效率高,但由于地震數(shù)據(jù)疊加處理方法技術(shù)本身存在的不足,疊后地震數(shù)據(jù)波阻抗反演在復(fù)雜地區(qū)的應(yīng)用受到限制。

反射地震波是地下波阻抗變化產(chǎn)生的,ZHANG等[16]提出了一種利用真振幅逆時偏移估計地下波阻抗擾動的方法,并認(rèn)為在角度域共成像點道集中的小角度成像波場信息反映了波阻抗的擾動,為利用波動方程偏移方法技術(shù)進(jìn)行地下波阻抗成像提供了思路,但是他們的地下波阻抗擾動估計方法是通過對偏移方法中成像公式的改造得到,而不是根據(jù)與波阻抗有關(guān)的反射波場傳播方程得到?；诓▌臃匠痰哪鏁r偏移對地下介質(zhì)模型空間變化具有很好的適應(yīng)性[17-22]。與其它偏移成像方法一樣,當(dāng)前的逆時偏移也是基于CLAERBOUT[23]所提出的偏移成像原理,并認(rèn)為對偏移成像結(jié)果進(jìn)行積分就可以實現(xiàn)對波阻抗的成像[24]。CLAERBOUT偏移成像的核心思想是其提出的“時間一致性成像原理和成像公式”,我們認(rèn)為其成像原理僅是一個有關(guān)地震波傳播旅行時的概念性描述,成像公式缺乏嚴(yán)格的數(shù)學(xué)物理推導(dǎo),也不滿足描述地震波傳播的數(shù)學(xué)物理方程。陳生昌等[25]將地震數(shù)據(jù)偏移成像定義為近似求解地震數(shù)據(jù)的波形線性反演問題,由地震波所滿足的波動方程,利用高頻近似給出了反射面上反射率的定義,進(jìn)而推導(dǎo)出了基于反射率分布的反射波傳播方程,然后在反射波傳播方程的基礎(chǔ)上應(yīng)用線性反演的方法對地震數(shù)據(jù)偏移成像進(jìn)行了重新推導(dǎo),建立了對地下反射率進(jìn)行成像的偏移成像公式,為基于地震波場傳播方程的偏移成像方法研究與應(yīng)用提供了堅實的理論基礎(chǔ)。

根據(jù)文獻(xiàn)[25]對偏移成像的重新定義與推導(dǎo),本文在變密度聲波方程的基礎(chǔ)上,首先應(yīng)用地震波傳播的散射理論和地震波方法理論中的高頻近似,推導(dǎo)出基于波阻抗相對變化的一次反射地震波線性傳播方程,然后再應(yīng)用線性反演的方法推導(dǎo)出波阻抗相對變化的成像公式,最后將提出的波阻抗成像方法應(yīng)用于簡單的層狀介質(zhì)模型和Marmousi模型的合成地震數(shù)據(jù),均取得了理想的結(jié)果。

1 方法原理

我們采用變密度聲波方程描述地震波傳播:

(1)

(2)

將方程(1)和方程(2)變換到頻率域,分別得到:

式中:ω表示頻率;F(ω)為f(t)的傅里葉變換;G(x,ω;xs)和P(x,ω;xs)分別為g(x,t;xs)和p(x,t;xs)的傅里葉變換。

1.1 基于波阻抗相對變化的反射波動方程

為了構(gòu)建反射地震數(shù)據(jù)的波阻抗相對變化成像方法,首先要推導(dǎo)出基于波阻抗相對變化的反射波動方程。

令v0(x)和ρ0(x)分別為v(x)和ρ(x)的背景速度模型和背景密度模型,則它們的相對擾動分別為:

(5)

將(5)式代入方程(3)可得到由相對擾動引起的散射波場傳播方程:

(6)

式中:Ps(x,ω;xs)為速度和密度相對擾動產(chǎn)生的散射波場;P(x,ω;xs)表示總波場,有P(x,ω;xs)=Pi(x,ω;xs)+Ps(x,ω;xs),其中,Pi(x,ω;xs)表示背景速度和密度模型下的波場,也稱為背景波場或入射波場,有:

(7)

假定v0(x)和ρ0(x)為光滑的背景模型,則方程(7)中的Pi(x,ω;xs)僅包含直達(dá)波和透射波,而不包含散射波和反射波。

在一階Born近似下,對于方程(6)有:

(8)

將方程(8)返回到時間域,并利用形如方程(2)背景模型中的Green函數(shù)將其寫為積分形式,有:

(9)

式中:v0(y)和ρ0(y)為光滑的背景模型;g(x,t;y)為背景模型中描述散射波傳播的Green函數(shù);“*”代表時間褶積。經(jīng)推導(dǎo)與簡化[7],方程(9)可寫為:

(10)

將方程(10)寫為頻率域形式,有:

(11)

(12)

進(jìn)一步化簡,有:

(13)

式中:Pr(x,ω;xs)為高頻近似下由地下速度、密度相對擾動產(chǎn)生的反射波場。

令:

(14)

根據(jù)波阻抗定義I(y)=v(y)ρ(y),對于背景模型有I0(y)=v0(y)ρ0(y),波阻抗擾動有δI(y)=δ·[v(y)ρ(y)]=ρ(y)δv(y)+v(y)δρ(y)。對于法向入射α=0,有:

(15)

式中:Ir(y,α=0)為法向入射時的波阻抗相對擾動。由(14)式可知,Ir(y,α)是一個隨半開角α變化的波阻抗相對擾動。隨著α的增大,Ir(y,α)由阻抗型相對擾動逐漸轉(zhuǎn)變?yōu)樗俣刃拖鄬_動。當(dāng)α=π/2時,有:

(16)

將(14)式代入方程(13),有:

(17)

將(17)式寫為時間域形式,有:

(18)

方程(17)和方程(18)分別為高頻近似條件下基于波阻抗相對變化的反射波動方程的頻率域形式和時間域形式。

1.2 基于反射波動方程的深度域波阻抗相對變化成像

1.1節(jié)在高頻近似條件下得到的基于波阻抗相對變化的反射波動方程(17)、(18)是利用反射波地震數(shù)據(jù)對地下波阻抗相對變化進(jìn)行反演和成像的基礎(chǔ)。我們將成像定義為基于正演算子的伴隨算子的近似反演(或稱為粗糙反演)。相對于基于逆算子的反演,利用成像得到的解的估計通常會存在保幅性差、分辨率低的問題。

利用背景介質(zhì)模型下的Green函數(shù),方程(7)中的入射波場可表示為:

(19)

在光滑背景模型的假定下,根據(jù)方程(17),對于炮點xs激發(fā),檢波點xg接收的反射波地震數(shù)據(jù)Dr(xg,ω;xs)可表示為:

(20)

根據(jù)陳生昌等[26]提出的時間二階積分波場的全波形反演方法,方程(20)可寫為:

(21)

令:

(22)

得到:

(23)

將線性積分方程(23)寫為矩陣形式,有:

(24)

(25)

式中:W-g代表W的最小二乘廣義逆矩陣;W*為W的伴隨矩陣。

將(25)式求得的m和方程(19)得到的入射波場Pi(x,ω;xs)代入(22)式,可得到波阻抗相對變化的估計:

(26)

式中:Re表示取實部運算。

考慮到(25)式計算的穩(wěn)定性和逆矩陣計算的復(fù)雜性,我們用伴隨矩陣W*代替(25)式中的最小二乘廣義逆矩陣W-g:

(27)

(27)式的運算實質(zhì)上是對波場數(shù)據(jù)u的伴隨傳播(反向傳播)。

再考慮(26)式中波場相除運算的不穩(wěn)定性問題,并結(jié)合(27)式將(26)式改寫為:

(28)

式中:“*”代表復(fù)共軛運算。

(19)式、(27)式和(28)式為波阻抗相對變化在角度α上的近似求解,實質(zhì)上實現(xiàn)了對Ir(x,α)的成像。將上述對Ir(x,α)的成像運算轉(zhuǎn)換到時間域。

1) 地下入射波場的計算(對應(yīng)公式(19)所表示的運算):

(29)

2) 地下時間二階積分反射波場的伴隨(逆時)重建(對應(yīng)公式(27)所表示的運算):

(30)

3) 對波阻抗相對變化Ir(x,α)的成像(對應(yīng)公式(28)所表示的運算):

(31)

(32)

對于不考慮地下介質(zhì)密度變化的標(biāo)量波動方程,則可得到基于速度相對變化vr(x,α)的反射波動方程,有vr(x,α)=av(x)cos2α=2δv(x)cos2α/v0(x),相應(yīng)地得到對速度相對變化的成像,對應(yīng)的計算公式如下。

1) 地下入射波場的計算:

(33)

2) 地下時間二階積分反射波場的伴隨(逆時)重建:

(34)

3) 對速度相對變化vr(x,α)的成像:

(35)

2 數(shù)值試驗

為了驗證上述基于反射波動方程提出的疊前地震反射數(shù)據(jù)波阻抗相對變化成像方法的正確性與有效性,我們采用兩個二維模型的合成數(shù)據(jù)進(jìn)行試驗。由于在模型試驗中沒有考慮密度的變化,因此試驗得到的波阻抗相對變化成像結(jié)果實際是速度相對變化的成像結(jié)果。

2.1 簡單模型試驗

圖1顯示了簡單層狀模型的速度分布,在速度值為3000m/s的均勻背景中包含一個低速層和一個高速層,上部低速層的速度值為2500m/s,下部高速層的速度值為3500m/s。在該速度模型上進(jìn)行地震數(shù)據(jù)數(shù)值模擬,得到試驗所需的合成地震反射數(shù)據(jù)。將本文所提出的波阻抗相對變化成像方法應(yīng)用于簡單模型的合成地震反射數(shù)據(jù),得到圖2所示的速度相對變化成像結(jié)果。圖3為圖2紅框中速度相對變化成像結(jié)果的波形放大圖。

對比圖1和圖2可以看出,波阻抗相對變化成像方法所得到的波阻抗相對變化成像結(jié)果中的低阻層和高阻層與模型中的低速層和高速層在空間位置和相對幅值變化方面都對應(yīng)很好,驗證了本文所提出的波阻抗相對變化成像方法的正確性。

2.2 Marmousi模型試驗

Marmousi模型地震反射數(shù)據(jù)是檢驗波動方程疊前深度偏移成像算法的標(biāo)準(zhǔn)模型數(shù)據(jù)。Marmousi模型數(shù)據(jù)共有240炮,每炮96個接收道,炮間距25m,道間距25m,記錄長度3s,時間采樣率4ms。

圖1 簡單模型的速度分布

圖2 簡單模型的速度相對變化成像結(jié)果

圖3 圖2紅框中的波形放大顯示

圖4顯示了Marmousi模型的速度相對變化。圖5 為利用本文提出的波阻抗相對變化成像方法得到的速度相對變化成像結(jié)果。為了對比,我們還利用逆時偏移方法對Marmousi模型數(shù)據(jù)進(jìn)行了構(gòu)造偏移成像,圖6為由逆時偏移得到的構(gòu)造成像結(jié)果。

對比圖4,圖5和圖6可以看出,由波阻抗相對變化偏移成像方法得到的速度相對變化成像結(jié)果能更清楚地反映地下巖性的變化,而由逆時偏移得到的構(gòu)造成像結(jié)果主要反映地下構(gòu)造界面的幾何結(jié)構(gòu)。

圖4 Marmousi模型的速度相對變化分布

圖5 Marmousi模型的速度相對變化成像結(jié)果

圖6 Marmousi模型的逆時偏移成像結(jié)果

3 結(jié)論

地震反射數(shù)據(jù)成像是一種基于反射波動傳播方程的波形線性反演問題的近似求解。建立基于波阻抗相對變化的反射波動方程,通過線性反演方法提出的利用逆時傳播的深度域波阻抗相對變化成像方法是對地震數(shù)據(jù)反射率成像方法的推廣,可為進(jìn)一步的巖性成像方法研究打下基礎(chǔ)。本文方法的計算量與當(dāng)前常用的逆時偏移成像方法的計算量基本相當(dāng),僅增加了計算量非常少的反射地震數(shù)據(jù)時間二階積分計算。與當(dāng)前對反射面空間位置和反射率變化進(jìn)行成像的偏移成像方法不同,本文提出的波阻抗相對變化成像方法可對地下波阻抗的相對變化進(jìn)行成像,可為地震數(shù)據(jù)的巖性處理解釋提供基礎(chǔ)數(shù)據(jù)。

[1] 李慶忠.走向精確勘探的道路——高分辨地震勘探系統(tǒng)工程剖析[M].北京:石油工業(yè)出版社,1993:1-196 LI Q Z.The way to obtain a better resolution in seismic prospecting:a systematical analysis of high resolution seismic exploration[M].Beijing:Petroleum Industry Press,1993:1-196

[2] 李慶忠.巖性油氣藏地震勘探若干問題討論(Ⅰ)[J].巖性油氣藏,2008,20(2):1-6 LI Q Z.Discussion on seismic exploration of lithologic reservoirs(Ⅰ)[J].Lithologic Reservoirs,2008,20(2):1-6

[3] LAVERGNE M.Pseudo-diagraphics de vitesse en offshore profound[J].Geophysical Prospecting,1975,23(4):695-711

[4] CONNOLLY P.Elastic impedance[J].The Leading Edge,1999,18(4):438-452

[5] WHITECOMBE N.Elastic impedance normalization[J].Geophysics,2002,67(1):60-62

[6] LAILLY P.The seismic inverse problem as a sequence of before stack migrations[J].Expanded Abstracts of Conference on Inverse Scattering,Theory and Application,Society for Industrial and Applied Mathematics,1983:206-220

[7] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[8] VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[9] LINDSETH O.Synthetic sonic logs:a process for stratigraphic interpretation[J].Geophysics,1979,44(1):3-26

[10] COOKE D A,SCHNEIDER W A.Generalized linear inversion of reflection seismic data[J].Geophysics,1983,48(6):665-676

[11] OLDENBURG D W,SCHEUER T,LEVY S.Recovery of the acoustic impedance from reflection seismograms[J].Geophysics,1983,48(10):1318-1337

[12] 楊文采.地震道的非線性混沌反演I:理論和數(shù)值試驗[J].地球物理學(xué)報,1993,36(2):222-232 YANG W C.Nonlinear chaotic inversion of seismic traces I:theory and numerical experiments[J].Chinese Journal of Geophysics,1993,36(2):222-232

[13] 楊文采.地震道的非線性混沌反演Ⅱ:關(guān)于Lyapunov指數(shù)和吸引子[J].地球物理學(xué)報,1993,36(3):376-387 YANG W C.Nonlinear chaotic inversion of seismic traces Ⅱ:about Lyapunov exponent and attractor[J].Chinese Journal of Geophysics,1993,36(3):376-387

[14] ROGERS S T,FANG J H,KARR C L,et al.Determination of lithology from well logs using a neural network[J].AAPG Bulletin,1992,76(5):731-739

[15] 楊立強,宋海濱,郝天珧.基于BP神經(jīng)網(wǎng)絡(luò)的波阻抗反演及應(yīng)用[J].地球物理學(xué)進(jìn)展,2005,20(1):34-37 YANG L Q,SONG H B,HAO T Y.Application of impedance inversion based on BP neural network[J].Progress in Geophysics,2005,20(1):34-37

[16] ZHANG Y,DUAN L,ROBERTS G.True amplitude reverse time migration:from eflectivity to velocity and impedance perturbations[J/OL].[2016-10-01].http:∥www.cgg.com/technicalDocuments/cggv_0000014832.pdf

[17] HEMON C.Wave equations and models[J].Geophysical Prospecting,1978,26(4):790-821

[18] BAYSAL E,KOSLOFF D,SHERWOOD J.Reverse time migration[J].Geophysics,1983,48(11):1514-1524

[19] MCMECHAN G.Migration by extrapolation of time-dependent boundary values[J].Geophysical Prospecting,1983,31(3):413-420

[20] WHITMORE D.Iterative depth migration by backward time propagation[J].Expanded Abstracts of 53rdAnnual Internat SEG Mtg,1983:382-385

[21] ETGEN J,GRAY S,ZHANG Y.An overview of depth imaging in exploration geophysics[J].Geophysics,2009,74(6):WCA5-WCA17

[22] LEVEILLE J,JONES I,ZHOU Z,et al.Subsalt imaging for exploration,production,and development:a review[J].Geophysics,2011,76(5):WB3-WB20

[23] CLAERBOUT J.Toward a unified theory of reflector mapping[J].Geophysics,1971,36(3):467-481

[24] ROBEIN E.Seismic imaging:a review of the techniques,their principles,merits and limitations[M].Houten:EAGE Publications,2010:1-244

[25] 陳生昌,周華敏.再論地震數(shù)據(jù)偏移成像[J].地球物理學(xué)報,2016,59(2):643-654 CHEN S C,ZHOU H M.Re-exploration to migration of seismic data[J].Chinese Journal of Geophysics,2016,59(2):643-654

[26] 陳生昌,陳國新.時間二階積分波場的全波形反演[J].地球物理學(xué)報,2016,59(10):3765-3776 CHEN S C,CHEN G X.Full waveform inversion of the second-order time integral wavefield[J].Chinese Journal of Geophysics,2016,59(10):3765-3776

(編輯：陳 杰)

Arelativeimpedancevariationimagingofpre-stackseismicreflectiondatabasedonreflectionwaveequation

CHEN Shengchang1,ZHOU Huamin2

(1.SchoolofEarthSciences,ZhejiangUniversity,Hangzhou310027,China;2.ChangjiangRiverScientificResearchInstituteofChangjiangWaterResourcesCommission,Wuhan430010,China)

Impedance of subsurface media is an important parameter for lithologic processing and interpretation of seismic data.In this paper,a relative impedance variation imaging in depth domain based on the reflection wave equation is proposed.Starting from the variable density acoustic wave equation,the proposed method derives a linear propagation equation of the primary reflection waves based on the relative impedance variation through high frequency approximation.Then,on the basis of the linear propagation equation,it derives the relative impedance variation imaging method and its formula by using the linear inversion theory.The relative impedance variation imaging results obtained by the proposed method are deficient in fidelity and resolution,owing to the limitation of the frequency band of the reflected seismic data and to the inherent limitations of the imaging methods,which are approximate solutions to the waveform linear inversion problem.However,the computational complexity is equivalent to that of reverse time migration.The result of the model tests on synthetic data showed the effectiveness of the proposed method.

pre-stack seismic reflection data,high frequency approximation,reflection wave equation,linear inversion,relative impedance variations,imaging

P631

：A

1000-1441(2017)05-0651-07DOI:10.3969/j.issn.1000-1441.2017.05.005

陳生昌,周華敏.基于反射波動方程的疊前地震反射數(shù)據(jù)波阻抗相對變化成像研究[J].石油物探,2017,56(5):657

CHEN Shengchang,ZHOU Huamin.A relative impedance variation imaging of prestack seismic reflection data based on reflection wave equation

[J].Geophysical Prospecting for Petroleum,2017,56(5):657

2016-10-11;改回日期：2017-04-06。

陳生昌(1965—),男,教授,博士生導(dǎo)師,主要從事勘探地球物理和計算地球物理研究。

國家自然科學(xué)基金項目(41074133,41374001)、國家科技重大專項(2016ZX05033002)共同資助。

This research is financially supported by the National Natural Science Foundation of China (Grant Nos. 41074133,41374001) and National Science and Technology Major Project of China (Grant No. 2016ZX05033002).