基于波動方程偏移的地震數(shù)據(jù)地層成像

陳生昌 劉亞楠 李代光

(浙江大學(xué)地球科學(xué)學(xué)院，浙江杭州 310027)

0 引言

隨著油氣勘探開發(fā)目標(biāo)由構(gòu)造型油氣藏轉(zhuǎn)化為巖性油氣藏、構(gòu)造巖性復(fù)合型油氣藏和非常規(guī)油氣藏，對地震數(shù)據(jù)反演成像技術(shù)的要求越來越高，不僅對地震數(shù)據(jù)全波形反演技術(shù)實用化的呼聲越來越強烈，而且要求以地下構(gòu)造成像為目標(biāo)的波動方程偏移方法應(yīng)保幅、保真，還希望地震數(shù)據(jù)反演結(jié)果能直接反映地層的物性變化。

以獲取地下波阻抗空間變化的地層物性反演成像技術(shù)自20世紀(jì)70年代以來得到了很大的發(fā)展，如疊后地震波阻抗反演[1-6]、基于Zoeppritz方程及其近似的彈性阻抗反演[7-13]和全波形反演[14-16]。疊后波阻抗反演雖然具有計算效率高的優(yōu)勢，也是當(dāng)前生產(chǎn)中比較常用的技術(shù)，但由于地震數(shù)據(jù)水平疊加處理固有的不足，因此疊后波阻抗反演方法在復(fù)雜地區(qū)的應(yīng)用受到限制。基于Zoeppritz方程及其近似的AVO和AVA反演方法直接在疊前數(shù)據(jù)道集上進(jìn)行，可有效地利用地震振幅隨炮檢距或入射角變化的信息，能用于多種物性和巖性參數(shù)的反演。由于Zoeppritz方程是一種描述平面波入射到無窮大平界面的反射、透射和波型轉(zhuǎn)換的理論[17]，因此該方程對于實際地震勘探中常見的非平面波和非平界面可能會不適應(yīng)。利用基于物性參數(shù)(如波阻抗)的波動方程全波形反演，數(shù)理嚴(yán)謹(jǐn)，可完整地利用地震波場的時空變化信息，但計算量大，對地震數(shù)據(jù)要求高(如要求數(shù)據(jù)有盡可能低的低頻信息和大炮檢距數(shù)據(jù))[18]，是油氣地震勘探近30年的研究熱點和亟待突破的技術(shù)。

給定一個具有地震波運動學(xué)特征準(zhǔn)確的地下介質(zhì)光滑模型，得到以地下反射率為模型參數(shù)的地震數(shù)據(jù)線性正演公式，地震偏移可視為一種有關(guān)地下反射率的線性反演[19-20]。Zhang等[21-22]基于真振幅逆時偏移提出了一種速度和波阻抗反演方法，并認(rèn)為角度域共成像點道集的小角度成像結(jié)果主要反映波阻抗的變化，大角度成像結(jié)果主要反映速度的變化，為利用波動方程偏移實現(xiàn)地下波阻抗反演提供了思路。但是該方法所用的真振幅逆時偏移中的成像條件是源于Bleistein等[19，23]利用Kirchhoff近似導(dǎo)出的成像公式。而Kirchhoff近似中反射波與入射波間的關(guān)系僅適合平面波小于臨界角入射到無窮大平界面的特殊情況，而不適合實際地震勘探中常見的非平面波和非平界面。陳生昌等[24-26]曾利用推導(dǎo)出的反射波線性依賴于阻抗相對擾動的反射波方程，開展了有關(guān)波阻抗相對變化成像、波阻抗擾動的近似反演和最小二乘反演方法的初步研究，為利用反射地震數(shù)據(jù)進(jìn)行地層成像打下了基礎(chǔ)。

針對以上有關(guān)波阻抗或地層物性參數(shù)反演成像存在的問題和關(guān)于地震數(shù)據(jù)波動方程偏移成像的認(rèn)識，以及對不同模型參數(shù)的地震數(shù)據(jù)正演公式會產(chǎn)生不同的反演目標(biāo)和反演算法的認(rèn)知[27]，本文基于波動方程偏移的概念，利用用于地震數(shù)據(jù)偏移的光滑介質(zhì)模型，根據(jù)地震波傳播理論導(dǎo)出了基于地下反射體波阻抗相對擾動的地震數(shù)據(jù)線性正演公式； 然后在此基礎(chǔ)上結(jié)合線性反演理論，研究利用地震數(shù)據(jù)對地下反射體的波阻抗相對擾動的近似反演，構(gòu)建以波阻抗變化為主要目標(biāo)的地震數(shù)據(jù)地層成像方法。將提出的地震數(shù)據(jù)地層成像方法分別應(yīng)用于標(biāo)量波方程描述的地震數(shù)據(jù)和聲波方程描述的地震數(shù)據(jù)，均取得了理想的成像結(jié)果。

1 基于波阻抗相對擾動的線性正演

在常規(guī)油氣地震勘探中，地表激發(fā)的地震波向地下傳播，當(dāng)?shù)卣鸩ㄓ龅椒蔷鶆蝮w或地層分界面時會產(chǎn)生散射或反射等次生波，再向上傳播到地表被檢波器接收。地震波在地下的傳播和散射(或反射)可用波動方程描述，如變密度聲波方程

=s(t)δ(x-xs)

(1)

式中：u(x,t;xs)為壓力波場；ρ(x)為密度場；v(x)為速度場；s(t)為震源函數(shù)；x=(x,y,z),為空間坐標(biāo)；xs為震源位置。

式(1)中，波場u(x,t;xs)與右端震源項呈線性關(guān)系(震源項不包含波場u(x,t;xs))，與地下介質(zhì)密度ρ(x)和速度v(x)呈非線性關(guān)系，這種非線性關(guān)系可表示為

u=F(ρ,v)

(2)

式中F是關(guān)于ρ和v的非線性函數(shù)。將隨空間變化的密度ρ(x)和速度v(x)用背景模型和擾動模型表示，即

(3)

式中：ρb(x)和vb(x)分別為介質(zhì)的背景密度和背景速度； Δρ(x)和Δv(x)分別為相對于背景密度和背景速度的擾動。本文不僅要求ρb(x)和vb(x)光滑變化，還要求在背景密度和背景速度下地震波的運動學(xué)特征與真實介質(zhì)中地震波一致。

如果取式(1)中的密度和速度為背景密度和背景速度，有

=s(t)δ(x-xs)

(4)

式中ub(x,t;xs)為背景密度和背景速度下的壓力波場，也稱為背景波場。如果式(4)右端的震源函數(shù)s(t)為脈沖函數(shù)，則該方程對應(yīng)的解即為背景密度和背景速度下的Green函數(shù)g(x,t;xs)。利用背景波場ub(x,t;xs)和Green函數(shù)g(x,t;xs)，并根據(jù)Lippmann-Schwinger方程的級數(shù)展開式和一階Born近似，式(2)可近似為

u=ub+gpub

(5)

式中p為式(1)中的波動算子與式(4)中的波動算子之差，即

(6)

式(5)右端第二項稱為由速度和密度擾動產(chǎn)生的一階Born近似擾動波場，有

up=u-ub=gpub

(7)

與式(5)所對應(yīng)的波動方程為

(8)

式中av和aρ分別為速度和密度的相對擾動，即

(9)

由于ρb(x)和vb(x)的光滑性，ub(x,t;xs)不包含散射波和反射波，因此在地表觀測情況下，有ub(xr,t;xs)=0，其中xr為檢波點坐標(biāo)。對于式(7)中一階Born近似的擾動波場，有

(10)

式中：Ω表示av和aρ的空間分布范圍； “*t”表示時間方向的褶積。利用分部積分法，式(10)可進(jìn)一步化簡為[11]

(11)

對于非均勻體產(chǎn)生的擾動波場，陳生昌等[20]根據(jù)非均勻體分布范圍Ω的尺寸與地震波長之間的相對大小關(guān)系，把相對于地震波長為小尺度的非均勻體產(chǎn)生的擾動波場稱為散射波場，相應(yīng)的非均勻體稱為散射體； 而把相對于地震波長為大尺度的非均勻體產(chǎn)生的擾動波場稱為反射波場，相應(yīng)的非均勻體稱為反射體。為了得到基于反射體波阻抗擾動的反射波場，根據(jù)文獻(xiàn)[26，28]，式(11)可表示為

(12)

式中：uR(xr,t;xs)為采集到的反射波場；σ為反射波傳播方向與入射波傳播方向間的反射開角；IR(x,σ)稱為波阻抗相對擾動，有IR(x,σ)=av+aρ×(1+cosσ)。對于沿反射面法向的入射及其反射，有σ=0，則

IR(x,σ=0)=av+2aρ

(13)

式中：Ib(x)=vb(x)ρb(x)，為背景介質(zhì)波阻抗； ΔI(x)=Δ[v(x)ρ(x)]=ρ(x)Δv(x)+v(x)Δρ(x)，為波阻抗擾動。式(12)即為基于波阻抗相對擾動的地震反射波場線性正演表示，對應(yīng)的波動方程為

(14)

在式(14)中，波場uR(x,t;xs)與右端源項成線性關(guān)系。

上述線性關(guān)系依賴于反射體波阻抗相對擾動的地震反射波場表達(dá)式(式(12))和式(14)，是本文利用反射地震數(shù)據(jù)獲取地下波阻抗空間變化的地層成像方法的理論基礎(chǔ)。

2 波動方程地層成像

地震數(shù)據(jù)的正演是地震數(shù)據(jù)反演的基礎(chǔ)。利用上述基于波阻抗相對擾動的地震數(shù)據(jù)線性正演公式，結(jié)合線性反演理論，可構(gòu)建以波阻抗擾動為目標(biāo)的地震數(shù)據(jù)地層成像方法。

式(12)兩邊對時間做Fourier變換，有

(15)

式中：上劃“～”表示在頻率域；ω表示圓頻率。令

(16)

則式(15)可寫為

(17)

將式(17)寫為矩陣形式，有

d=Lm

(18)

由式(18)中各個矩陣的數(shù)學(xué)物理性質(zhì)可知，求解式(18)中的m相當(dāng)于反演源(產(chǎn)生反射波場的虛源)[29]。根據(jù)線性最小二乘反演理論，式(18)的最小二乘解為

m=L-1d=(L*L)-1L*d

(19)

式中L-1、L*分別為傳播矩陣L的最小二乘逆和伴隨矩陣。在地震數(shù)據(jù)的反演中，特別是三維地震數(shù)據(jù)的反演中，矩陣L的規(guī)模通常非常大，(L*L)-1的計算在當(dāng)前的計算機條件下難以直接實現(xiàn)。因此，在地震數(shù)據(jù)的反演中通常把式(19)近似(即用L*代替L-1)為

ma=L*d

(20)

式中ma為m的近似結(jié)果。

式(20)相比于式(19)，不僅計算量大為減少，而且計算更穩(wěn)定。式(19)通常被稱為對m的最小二乘反演，式(20)為對m的近似反演或成像。地震數(shù)據(jù)的逆時偏移就是采用類似式(20)的運算實現(xiàn)對地下反射率的成像。式(20)中L*對波場d的運算表示波場d的伴隨傳播。這種波場的伴隨傳播，在地震數(shù)據(jù)的反演成像中，是通過波場的逆時傳播實現(xiàn)的[16]。

由于矩陣m中各個元素的組成不僅包含有x處待成像的波阻抗相對擾動，還含有頻率域的二階導(dǎo)數(shù)和x處已知的背景波場以及背景密度與背景速度，因此需要從ma中消除背景波場、背景密度和背景速度，并處理好頻率域的二階導(dǎo)數(shù)運算，才能得到x處波阻抗相對擾動的成像結(jié)果。

根據(jù)上述的推導(dǎo)和論述，關(guān)于IR(x,σ)成像的計算步驟歸納如下：

(1)應(yīng)用式(3)計算背景波場ub(x,t;xs)；

(2)利用伴隨傳播對式(14)右端虛源項成像，即地表波場uR(xr,t;xs)的逆時傳播

(21)

(22)

式中Tmax為地震數(shù)據(jù)的記錄長度。

上述波阻抗相對擾動成像的計算步驟與地震數(shù)據(jù)逆時偏移基本相同，僅在第(3)步的成像公式存在差別。本文的成像公式(式(22))是根據(jù)本文的地震反射波場線性正演公式(式(12))和反射波動方程(式(14))得到的。

為了得到反射開角σ域的成像結(jié)果IR(x,σ)，在應(yīng)用成像公式(式(22))之前，需要先對x處的波場uR(x,t;xs)和ub(x,t;xs)進(jìn)行角度分解[30-31]，得到角度域的波場uR(x,t,θR;xs)和ub(x,t,θI;xs)，其中θR為反射角(即z軸方向沿逆時針與反射波傳播方向間的夾角)，θI為入射角(即z軸方向沿逆時針與入射波傳播方向間的夾角)，有σ=π+θI-θR。IR(x,σ)可進(jìn)一步表示為

(23)

因此式(22)的成像公式在角度域可寫為

(24)

對于式(24)，如果σ比較小，則成像結(jié)果主要反映波阻抗的相對擾動ΔI(x)/Ib(x)； 如果σ比較大，則成像結(jié)果主要反映速度的相對擾動Δv(x)/vb(x)。利用角度域成像公式(式(24))得到的角度域共成像道集是進(jìn)一步物性和巖性分析、反演的基礎(chǔ)。

如果在成像過程中不考慮密度變化，僅考慮地下介質(zhì)的速度變化，則可用標(biāo)量波方程代替聲波方程描述地震波，即

(25)

對應(yīng)的地層性成像結(jié)果是地下介質(zhì)的速度相對擾動av(x)，相應(yīng)的計算步驟為：

(1)計算背景波場ub(x,t;xs)

(26)

(2)地表反射波場uR(xr,t;xs)的逆時傳播

(27)

(3)速度相對擾動av(x)的成像

(28)

由于標(biāo)量波方程中的波場是標(biāo)量，不涉及方向，因此利用標(biāo)量波波動方程的物性成像方法只能得到角度域平均的速度相對擾動，而不能得到類似于式(24)的角度域成像結(jié)果。

式(3)、式(21)、式(22)和式(24)表示以波阻抗相對擾動為目標(biāo)的地震數(shù)據(jù)地層成像方法，式(26)～式(28)為以速度相對擾動為目標(biāo)的地震數(shù)據(jù)地層成像方法。本文的地震數(shù)據(jù)地層成像方法雖然是以地震數(shù)據(jù)線性表示為基礎(chǔ)，但實現(xiàn)過程與波動方程逆時偏移一致，都需波場外推和成像。如果不進(jìn)行波場角度分解的角度域成像，則地層成像方法的計算量與逆時偏移基本相當(dāng)。此外，本文的地層成像方法對背景模型的要求與逆時偏移也是相同的，即要求有一個空間光滑變化的背景模型，背景模型的地震波運動學(xué)特征與真實模型一致。

基于上述認(rèn)識，將本文提出的方法稱為基于波動方程偏移的地震數(shù)據(jù)地層成像方法。

3 數(shù)值試驗

為了驗證本文的地層成像方法的正確性和有效性，對兩個模型的合成數(shù)據(jù)進(jìn)行試驗。由于方法的計算步驟和計算公式以及對模型的要求與逆時偏移相同，因此主要驗證方法的正確性。

3.1 標(biāo)量波動方程的速度相對擾動成像

假定地下介質(zhì)的密度為常數(shù)，僅考慮地下介質(zhì)的速度變化，使用標(biāo)量波動方程描述地震波的傳播。圖1為用于合成地震數(shù)據(jù)的速度模型，是部分Sig-sbee2A模型。模型尺寸為3km×3km，網(wǎng)格間距為15m。應(yīng)用波動方程有限差分方法合成地震數(shù)據(jù)，震源為主頻20Hz的Ricker子波，炮間距為30m，共101炮； 道間距為15m，每炮有201道。圖2為圖1經(jīng)過平滑、用于速度相對擾動成像的背景速度模型，圖3為由圖1和圖2得到的速度相對擾動的理論模型。

圖1 Sigsbee2A部分速度模型

圖2 Sigsbee2A背景速度模型

利用本文提出的標(biāo)量波地震數(shù)據(jù)地層性成像方法，對合成地震數(shù)據(jù)進(jìn)行速度相對擾動的地層成像(圖4)。對比圖4與圖3可見，圖4的地層成像結(jié)果準(zhǔn)確地反映了圖3的速度相對擾動變化。由于模型邊部的地震波照明能量不足，致使模型邊部的成效果較差。為了與模型理論值進(jìn)行仔細(xì)對比，提取圖3和圖4中x=1500m處的單道曲線，如圖5所示，可見，成像結(jié)果中的速度相對擾動與模型的理論速度相對擾動具有較高的一致性。由此可知，本文提出的速度相對擾動成像方法可以有效地實現(xiàn)對地層速度擾動的成像。

圖3 Sigsbee2A模型理論速度相對擾動

圖4 Sigsbee2A模型的速度相對擾動成像結(jié)果

圖5 Sigsbee2A模型速度擾動成像結(jié)果(紅)與理論值(藍(lán))的單道對比

3.2 聲波方程的聲阻抗相對擾動成像

同時考慮地下介質(zhì)的密度和速度變化，使用聲波方程描述地震波的傳播。圖6和圖7分別為用于合成地震數(shù)據(jù)的隨機層狀速度和密度模型，尺寸是3000m×1500m，兩個方向的網(wǎng)格間距均為10m。應(yīng)用波動方程有限差分方法合成地震數(shù)據(jù)，震源是主頻為20Hz的Ricker子波，炮間距為20m，共151炮； 道間距為10m，每炮有301道。圖8和圖9分別為用于地層成像的光滑背景速度模型和光滑背景密度模型。

圖6 隨機層狀速度模型

圖7 隨機層狀密度模型

圖8 隨機層狀模型的背景速度

圖9 隨機層狀模型的背景密度

根據(jù)給定的速度、密度模型和背景模型可計算得到波阻抗相對擾動理論模型(圖10)和速度相對擾動理論模型(圖11)。

圖10 隨機層狀模型的理論波阻抗相對擾動

圖11 隨機層狀模型的理論速度相對擾動

圖12 隨機層狀模型角度域平均波阻抗相對擾動成像結(jié)果

圖13 隨機層狀模型角度域平均的波阻抗相對擾動成像結(jié)果(紅色)與理論波阻抗相對擾動(藍(lán)色)的單道對比

圖14 隨機層狀模型x=1500m處的角度域共成像點道集

圖15 隨機層狀模型x=1500m處σ/2=0°對應(yīng)的波阻抗相對擾動成像結(jié)果(紅色)與理論波阻抗相對擾動(藍(lán)色)的單道對比

圖16 隨機層狀模型x=1500m處σ/2=50°對應(yīng)的波阻抗相對擾動成像結(jié)果(紅色)與理論速度相對擾動(藍(lán)色)的單道對比

由上述數(shù)值實驗結(jié)果可知，本文提出的波阻抗相對擾動成像方法可以有效地實現(xiàn)對地層波阻抗擾動的成像，得到的角度域平均成像結(jié)果和角度域共成像點道集不同角度的成像結(jié)果可作為地層波阻抗相對擾動或速度相對擾動的估計。本文的地層成像方法只需地震波運動學(xué)特征準(zhǔn)確的地下介質(zhì)光滑模型，就能具有與逆時偏移方法一樣的實用性。

4 結(jié)論與討論

通過引入具有地震波運動學(xué)特征準(zhǔn)確的地下介質(zhì)光滑模型，波動方程偏移是一種基于地震波傳播理論的地震數(shù)據(jù)線性反演，其目標(biāo)是對地層分界面反射率的成像。本文根據(jù)地震數(shù)據(jù)波動方程偏移成像的思想，利用用于偏移的地下介質(zhì)光滑模型，推導(dǎo)了基于地下地層波阻抗相對擾動的地震數(shù)據(jù)線性正演公式，結(jié)合線性反演理論，構(gòu)建了面向地層波阻抗相對擾動的地震數(shù)據(jù)地層成像方法。如果在成像中不對波場進(jìn)行角度分解，地層成像方法對于常規(guī)小炮檢距觀測系統(tǒng)的反射地震數(shù)據(jù)可以快捷地獲取角度域平均的地層成像結(jié)果。如果在成像中對波場進(jìn)行角度分解，地層成像方法雖然可以得到角度域共成像點道集和不同角度的地層成像結(jié)果，但計算量大增，尤其是三維地震數(shù)據(jù)的地層成像。波阻抗相對擾動的角度域共成像點道集可為發(fā)展精細(xì)的物性巖性分析、反演奠定基礎(chǔ)。地層成像方法的數(shù)值試驗結(jié)果表明，標(biāo)量波方程對應(yīng)的速度相對擾動成像和聲波方程對應(yīng)的波阻抗相對擾動成像均取得了理想的效果。

本文的地層成像具有與逆時偏移相同的計算步驟和相當(dāng)?shù)挠嬎懔?，適用于逆時偏移的地震數(shù)據(jù)同樣適用于本文提出的地層成像方法，因此本文的地層成像方法對于實際地震資料具有與逆時偏移一樣的實用性。