地震散射偏移方法

胡自多 雍學善 劉 威 陳生昌 王艷香

(①中國石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020； ②浙江大學地球科學學院，浙江杭州 310027)

0 引言

地震數據偏移成像是當前油氣勘探開發(fā)中獲取地下三維地質結構的主要方法[1-2]。自20世紀70年代初期，Claerbout[3-4]提出基于波動方程的現(xiàn)代偏移成像概念以來，有關偏移成像的研究與應用一直是地震數據處理技術的熱點之一。在Claerbout波動方程偏移成像概念的基礎上，Berkhout[5]提出適用于地震數據偏移成像的“WRW”地震反射波傳播模型。經過40余年的發(fā)展，地震數據偏移已從疊后走向了疊前，從二維走向了三維，從時間域走向了深度域，從各向同性介質走向了各向異性介質，從主要關注相位信息的構造成像走向了相位和振幅并舉的保幅成像[6-21]。

當前的地震數據偏移不區(qū)分散射和反射[22]，在偏移中主要用反射的概念，而在最小二乘偏移中主要應用散射的概念[23]。偏移的目標是對地下非均勻體邊界位置進行成像，即對產生散射的散射體位置和產生反射的反射面位置進行成像[24-31]。根據Berkhout的“WRW”概念模型，共炮道集地震數據偏移分兩步實現(xiàn)，一是基于波動方程的震源波場正向傳播和記錄波場反向傳播；二是利用正向傳播得到的入射波場和反向傳播得到的散射波場或反射波場實現(xiàn)對非均勻體空間位置的成像，即實現(xiàn)對非均勻體的構造成像。常規(guī)波動方程偏移成像中，沒有考慮散射波與反射波在傳播上的不同，散射波或反射波成像結果存在相位(極性)偏差。

本文首先討論了常規(guī)波動方程偏移方法存在的不足，然后從波動方程的擾動形式出發(fā)，根據地下非均勻體大小與波長之間的關系定義產生散射波的散射體和產生反射波的反射體，在具有地震波運動學特征準確的偏移速度模型下，建立一次散射波數據的線性正演方程，再利用線性反演理論推導出對散射體空間位置進行成像的地震散射數據偏移計算公式，最后應用國際上常用的Sigbee2A模型數值模擬數據驗證本文所提出的散射數據偏移方法。

1 常規(guī)波動方程偏移的主要不足

共炮集地震數據的逆時偏移，具體計算步驟概括如下：

(1)利用偏移速度模型計算入射波場

(1)

(2)利用逆時外推重建地下反射波場(或散射波場)

(2)

(3)利用成像條件(反射系數成像條件)進行波場成像

(3)

在地震波運動學特征(旅行時)準確的偏移速度模型下，利用式(1)和式(2)可得到入射波波場和重建的反射波場(或散射波場)，但反射系數成像公式(式(3))僅適用于無窮大平邊界和臨界角下平面波的反射，要求ur(x,t;xs)與ui(x,t;xs)具有相同的波形。在彎曲邊界或非平面波的情況下，反射系數與頻率有關[32]，即此時ur(x,t;xs)與ui(x,t;xs)的波形不同，會存在相位差異。由于產生散射波的非均勻體和無特定傳播方向的散射波場不滿足成像公式(式(3))所要求的“無窮大平反射邊界和臨界角下的平面波”條件，因此式(1)～式(3)構成的逆時偏移方法應用于散射數據是不合適的，得到的偏移成像結果會存在相位誤差。

為了構建適合散射數據的偏移成像計算公式，首先必須建立適合偏移成像的散射數據正演方程。

2 散射波正演方程

為了推導的簡便和不失一般性，本文僅考慮常密度聲波方程，即標量波動方程

(4)

式中：u為地震波場；v為速度。根據擾動理論，將速度表示為背景速度與速度擾動之和，即

v(x)=vb(x)+Δv(x)

(5)

波場表示為背景波場與擾動波場之和，即

u(x,t)=ub(x,t)+Δu(x,t)

(6)

式(4)可寫為擾動形式的波動方程形式

(7)

對于背景波場ub(x,t)，有

(8)

如果把式(7)的右端項視為產生擾動場的虛源項，是由波場與速度擾動相互作用產生的，因此與速度擾動體的尺寸和速度擾動量有關。為了合理地描述速度擾動體尺寸對擾動波場的作用，以地震波場主頻對應的波長度量速度擾動體尺度。

如果速度擾動Δv(x)的尺寸或者其空間變化的特征尺度a相對于波長λ較小，即a/λ≤1，則本文把該情況下的速度擾動體視為局部散射體，相應地式(7)的右端項可近似為點源項，產生的擾動波場可視為散射波場。如果速度擾動Δv(x)的尺寸或者其空間變化的特征尺度a相對于波長λ較大，即a/λ>1，則速度擾動體是一種在空間具有一定延續(xù)度的非均勻體，本文稱之為反射體，相應地式(7)的右端項可視為一種集成源項，產生具有特定方向性的反射波場，本文不考慮這種情況。上述的“速度擾動Δv(x)的尺寸或者其空間變化的特征尺度相對于波場的波長比較小”，在物理實質上是一種局部近似條件，本文稱為局部散射近似條件。

在局部散射近似條件下，擾動波場Δu(x，t)轉化為散射波場us(x，t)，即Δu(x，t)→us(x，t)，相應地把波場擾動方程(式(7))改寫為散射波正演方程

(9)

再利用散射波場相對于背景波場的小擾動近似(一階Born近似)，式(9)退化為描述散射波場與速度擾動Δv(x)成線性關系的一次散射波正演方程，即

(10)

由于現(xiàn)有的偏移方法不區(qū)分散射和反射，因此在波動方程的最小二乘偏移和常規(guī)偏移方法的正演表達推導中會出現(xiàn)與上述散射波正演方程類似的推導[23]，并把速度相對擾動Δv(x)/vb(x)定義為反射系數。由于反射系數在數學上是一個Dirac函數，而速度相對擾動在數學上是一個階躍函數。因此本文認為在偏移中不區(qū)分散射與反射，并直接用Born近似下的散射波正演方程作為波動方程的最小二乘偏移和常規(guī)偏移方法的正演表達，在數學和物理意義上欠嚴謹。

3 散射數據偏移

在局部散射近似和一階Born近似條件下推導出的一次散射波正演方程式(10)是進行一次散射波模擬和一次散射波反演成像的基礎方程，也是本文散射數據偏移的基礎方程。

為便于推導，利用式(8)所對應的Green函數將式(10)轉化為積分形式，即

(11)

(12)

(13)

線性積分方程(式(11))可寫為

(14)

因此求解式(11)對應為一個線性反演源問題(Inverse Source Problem)[33]。把式(14)寫為矩陣方程形式，有

d=Gm

(15)

式中：d為一次散射數據組成的矩陣；G為由Green函數構成的矩陣；m為虛源項形成的矩陣。利用最小二乘反演，可得到m的解

m=(GTG)-1GTd

(16)

由于(GTG)-1的數值計算難度問題，把式(16)近似為

m=GTd

(17)

利用式(17)得到的近似m，是對式(13)的虛源s(y,t;xs)的近似求解。再利用由式(8)得到的背景波場ub(y,t;xs)，就可得到對速度擾動Δv(y)的近似估計

(18)

考慮到式(18)除法運算的穩(wěn)定性問題，根據常規(guī)偏移中把反褶積成像條件轉化為互相關成像條件的處理辦法，把式(18)中的分母項轉化到分子項，得到散射波互相關成像條件公式，即

(19)

由于式(19)與式(17)近似，利用上述公式得到的速度擾動反演結果雖然只是真實速度擾動的近似，但由于速度擾動體的局部假定，該反演結果還是能清楚反映速度擾動體(散射體)的空間位置，即實現(xiàn)對散射體空間位置的成像。

式(17)是利用Green函數矩陣的轉置，把散射波場反向傳播到產生散射波的虛源處，在數值計算中可利用地震波的逆時外推方法，把記錄的散射波場逆時外推到入射波場(背景波場)作用于散射體產生散射虛源時刻的波場。式(19)中的s(y,t;xs)與背景波場ub(y,t;xs)具有相同的時間(即入射波作用于散射體時的時間)。

綜上所述，可以得到下述對散射體空間位置進行成像的計算公式，即散射數據的偏移計算公式。

(1)根據式(17)，得到散射數據的逆時傳播近似反演散射虛源的計算公式

(20)

(2)根據式(8)可得地下背景波場計算公式

(21)

(3)根據式(19)，得到對散射體空間位置進行成像的計算公式

(22)

式(20)是散射數據從Tmax到零時刻的反向傳播，式(21)是地震波場從零時刻到Tmax的正向傳播。

對比本文推導的散射數據偏移方法與常規(guī)反射波偏移方法，可見兩者的主要差異是成像條件不同(式(22)與式(3))。式(22)是基于一次散射波方程式(10)推導得到的，而式(3)是基于無窮大平邊界和平面波的反射波傳播理論推導得到的。散射數據偏移方法是對速度擾動的近似成像，由于速度擾動體在地震波長尺度下的空間局部性，速度擾動Δv(x)在空間上可近似為一個Dirac函數，因此通過對速度擾動的近似成像實現(xiàn)對速度擾動體的位置成像[29]，而常規(guī)逆時偏移方法是對反射面的反射系數成像。其次是散射數據的逆時傳播方程(式(20)與式(2))不同，式(20)是非齊次波動方程(有源項)的波場逆時傳播，式(2)是齊次波動方程(無源項)邊值條件的波場逆時傳播。本文利用局部散射近似條件推導出的散射數據逆時偏移方法，與當前不區(qū)分散射和反射而直接把常規(guī)逆時偏移方法作為散射數據偏移方法相比，更嚴謹合理。

由于式(22)較式(3)增加了對入射波場的時間二階導數算子(該算子不僅可以提高成像的分辨率，相位也會發(fā)生180°旋轉)，因此本文的散射數據偏移方法較常規(guī)散射數據偏移方法在理論上會有更高的分辨率和準確的相位。對比本文提出的散射數據偏移方法的計算公式(式(20)～式(22))與常規(guī)逆時偏移方法的計算公式(式(1)～式(3))可知，本文算法的運算成本與常規(guī)方法相當。

4 數值試驗

為了驗證本文基于散射波方程推導出的散射偏移方法的正確性與有效性，應用Sigbee2A模型進行偏移成像試驗。該模型(圖1)在沉積地層中包含高速鹽體和許多高速散射體(圖1深層中分布的黑色孤立體)，非常適合對本文散射偏移方法的驗證。圖2為用于偏移成像試驗的背景速度模型。模型大小為1067×1201個網格點，水平和垂直方向的網格間距分別為22.86m和7.62m，觀測方式為海上拖纜觀測系統(tǒng)。共模擬500炮，炮間距為45.72m，拖纜長度為7932.42m，檢波點間距為22.86m，每炮道數不等，最少為57道，最多為348道，時間采樣間隔為8ms，時間采樣點數為1500。

分別應用常規(guī)逆時偏移成像方法和本文提出的散射地震偏移成像方法對Sigbee2A模型數據進行偏移成像(圖3、圖4)。

圖1 Sigbee2A速度模型縱、橫坐標為網格點序號，圖2～圖7同

圖2 Sigbee2A偏移速度模型

對比圖1、圖3和圖4可以看出，不論常規(guī)逆時偏移方法還是本文的散射偏移方法都能對Sigbee2A模型中的高速鹽體、反射層和高速散射體進行清晰成像，相對于常規(guī)逆時偏移方法的成像結果(圖3)，散射波偏移方法的成像結果(圖4)中同相軸較細、分辨率較高。由圖3、圖4中的鹽體邊界成像結果可以看出，兩者是反相的，存在180°的相位差異，這與本文的理論是一致的，但兩者的相位都是錯的[34]，都存在90°的相位誤差。對于反射體邊界的偏移成像需要使用針對反射數據的偏移成像方法[18，34]，才能獲得正確相位的反射面偏移成像結果(圖3相位旋轉90°，圖4相位旋轉-90°)。對于散射體的成像，散射偏移成像方法得到的散射體的偏移成像結果要好于常規(guī)偏移成像方法。

為了仔細對比散射偏移方法和常規(guī)逆時偏移方法的成像效果，截取圖1速度模型和圖3、圖4兩種偏移成像結果的左邊部分進行局部放大，如圖5～圖7所示。常規(guī)逆時偏移對散射體的成像出現(xiàn)了反極性，高速散射體成像為低速散射體(圖6)；在圖7的散射偏移結果中，散射體極性正確，高速散射體成像為高速散射體。與圖6所示的常規(guī)逆時偏移成像結果相比，散射偏移結果中高速散射體不僅分辨率高，而且位置準確、相位(極性)正確。這一結果證明了本文提出的散射偏移方法對散射體成像的有效性。

圖3 Sigbee2A模型的常規(guī)逆時偏移結果

圖4 Sigbee2A模型的散射偏移結果

圖5 速度模型的局部放大

圖6 常規(guī)逆時偏移成像結果的局部放大

圖7 散射偏移成像結果的局部放大

5 結論

地震散射波是地下相對于地震波長較小的局部非均勻體產生的，而地震反射波是地下相對于地震波長在空間上具有一定延續(xù)度的非均勻體產生的。不考慮散射波場與反射波場之間的特征差異，直接把基于反射波概念的逆時偏移方法應用于散射數據的偏移成像，缺乏理論上的嚴謹性。本文利用波動方程的擾動形式，通過地下非均勻體與地震波長之間的相對關系，得到了描述一次散射波的非齊次波動方程，并把該方程作為散射數據偏移成像的線性正演方程；然后再利用反源問題的求解方法反演散射波的地下散射虛源，推導出散射數據的逆時偏移公式。對于散射體產生的散射數據，本文提出的散射偏移方法相對于當前的常規(guī)逆時偏移方法，可以得到更高分辨率的散射體成像結果和準確的相位(極性)，并在Sigbee2A模型高速散射體的成像中取得理想效果。