基于波動方程偏移的地震數據地層成像

陳生昌 劉亞楠 李代光

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

0 引言

隨著油氣勘探開發目標由構造型油氣藏轉化為巖性油氣藏、構造巖性復合型油氣藏和非常規油氣藏，對地震數據反演成像技術的要求越來越高，不僅對地震數據全波形反演技術實用化的呼聲越來越強烈，而且要求以地下構造成像為目標的波動方程偏移方法應保幅、保真，還希望地震數據反演結果能直接反映地層的物性變化。

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

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

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

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

在常規油氣地震勘探中，地表激發的地震波向地下傳播，當地震波遇到非均勻體或地層分界面時會產生散射或反射等次生波，再向上傳播到地表被檢波器接收。地震波在地下的傳播和散射(或反射)可用波動方程描述，如變密度聲波方程

=s(t)δ(x-xs)

(1)

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

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

u=F(ρ,v)

(2)

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

(3)

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

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

=s(t)δ(x-xs)

(4)

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

u=ub+gpub

(5)

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

(6)

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

up=u-ub=gpub

(7)

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

(8)

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

(9)

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

(10)

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

(11)

對于非均勻體產生的擾動波場，陳生昌等[20]根據非均勻體分布范圍Ω的尺寸與地震波長之間的相對大小關系，把相對于地震波長為小尺度的非均勻體產生的擾動波場稱為散射波場，相應的非均勻體稱為散射體； 而把相對于地震波長為大尺度的非均勻體產生的擾動波場稱為反射波場，相應的非均勻體稱為反射體。為了得到基于反射體波阻抗擾動的反射波場，根據文獻[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)，為背景介質波阻抗； ΔI(x)=Δ[v(x)ρ(x)]=ρ(x)Δv(x)+v(x)Δρ(x)，為波阻抗擾動。式(12)即為基于波阻抗相對擾動的地震反射波場線性正演表示，對應的波動方程為

(14)

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

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

2 波動方程地層成像

地震數據的正演是地震數據反演的基礎。利用上述基于波阻抗相對擾動的地震數據線性正演公式，結合線性反演理論，可構建以波阻抗擾動為目標的地震數據地層成像方法。

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

(15)

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

(16)

則式(15)可寫為

(17)

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

d=Lm

(18)

由式(18)中各個矩陣的數學物理性質可知，求解式(18)中的m相當于反演源(產生反射波場的虛源)[29]。根據線性最小二乘反演理論，式(18)的最小二乘解為

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

(19)

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

ma=L*d

(20)

式中ma為m的近似結果。

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

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

根據上述的推導和論述，關于IR(x,σ)成像的計算步驟歸納如下：

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

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

(21)

(22)

式中Tmax為地震數據的記錄長度。

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

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

(23)

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

(24)

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

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

(25)

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

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

(26)

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

(27)

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

(28)

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

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

基于上述認識，將本文提出的方法稱為基于波動方程偏移的地震數據地層成像方法。

3 數值試驗

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

3.1 標量波動方程的速度相對擾動成像

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

圖1 Sigsbee2A部分速度模型

圖2 Sigsbee2A背景速度模型

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

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

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

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

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

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

圖6 隨機層狀速度模型

圖7 隨機層狀密度模型

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

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

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

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

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

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

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

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

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

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

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

4 結論與討論

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

本文的地層成像具有與逆時偏移相同的計算步驟和相當的計算量，適用于逆時偏移的地震數據同樣適用于本文提出的地層成像方法，因此本文的地層成像方法對于實際地震資料具有與逆時偏移一樣的實用性。