王維紅， 井洪亮

東北石油大學(xué) 地球科學(xué)學(xué)院， 大慶 163318

?

基于稀疏反演三維表面多次波壓制方法

王維紅， 井洪亮

東北石油大學(xué) 地球科學(xué)學(xué)院， 大慶 163318

三維表面多次波壓制是海洋地震資料預(yù)處理中的重要研究課題，基于波動(dòng)理論的三維表面多次波壓制方法(3D SRME)是數(shù)據(jù)驅(qū)動(dòng)的方法，理論上來(lái)說(shuō)，可有效壓制復(fù)雜構(gòu)造地震數(shù)據(jù)表面多次波.但該方法因?qū)υ嫉卣饠?shù)據(jù)采集要求高而很難在實(shí)際資料處理中廣泛應(yīng)用.本文基于貢獻(xiàn)道集的概念，將稀疏反演方法引入到表面多次波壓制中，應(yīng)用稀疏反演代替橫測(cè)線積分求和，無(wú)需對(duì)橫測(cè)線進(jìn)行大規(guī)模重建，進(jìn)而完成三維表面多次波預(yù)測(cè)，這樣可有效解決實(shí)際三維地震數(shù)據(jù)橫測(cè)線方向稀疏的問(wèn)題.基于縱測(cè)線多次波積分道集為拋物線的假設(shè)，為保證預(yù)測(cè)后三維表面多次波和全三維數(shù)據(jù)預(yù)測(cè)的多次波在運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征上基本一致，文中對(duì)預(yù)測(cè)數(shù)據(jù)實(shí)施基于穩(wěn)相原理的相位校正.理論模型和實(shí)際數(shù)據(jù)的測(cè)試結(jié)果表明，本文基于稀疏反演三維表面多次波壓制方法可在橫測(cè)線稀疏的情況下，有效壓制三維復(fù)雜介質(zhì)地震資料中的表面多次波，從而更好地提高海洋地震資料的信噪比，為高分辨率地震成像提供可靠的預(yù)處理數(shù)據(jù)保障.

稀疏反演; 表面多次波; 貢獻(xiàn)道集; 多次波預(yù)測(cè)

1 引言

多次反射波常出現(xiàn)于海洋地震勘探資料中，而且通常能量很強(qiáng)，從而對(duì)一次波的處理和成像造成強(qiáng)烈干擾，所以對(duì)當(dāng)前廣泛應(yīng)用的成像方法而言，多次波壓制預(yù)處理尤為重要(Jakubowicz，1998；李翔和胡天躍，2009；王保麗等，2014).在眾多的多次波壓制方法中，基于數(shù)據(jù)驅(qū)動(dòng)的3D SRME(Surface-related multiple elimination)方法被認(rèn)為是最有可能成功預(yù)測(cè)和壓制表面多次波的方法(Dragoset et al.，2010).基于Berkhout(1980)提出的數(shù)據(jù)矩陣的思想，該方法最初由Verschuur等提出(1992).依據(jù)惠更斯原理(Berkhout and Verschuur,1997; Verschuur and Berkhout, 1997; Sun, 1999)，SRME方法分為兩步來(lái)實(shí)現(xiàn)：一是多次波預(yù)測(cè)；二是多次波自適應(yīng)相減.而多次波預(yù)測(cè)過(guò)程又可分為兩步：一是計(jì)算多次波的貢獻(xiàn)道集(Kabir and Verschuur,1995；van Dedem and Verschuur, 2001；Hokstad and Sollie, 2003)；二是對(duì)貢獻(xiàn)道集進(jìn)行疊加，形成預(yù)測(cè)的多次波道集(van Dedem and Verschuur, 2005).SRME方法雖基于數(shù)據(jù)驅(qū)動(dòng)，可處理復(fù)雜構(gòu)造地震資料，但有效預(yù)測(cè)多次波的一個(gè)重要假設(shè)是輸入全地震數(shù)據(jù)(Verschuur，1991；Dragoset et al.，2006)，也就是說(shuō)每個(gè)檢波點(diǎn)位置都有炮點(diǎn)(炮點(diǎn)間隔等于檢波點(diǎn)間隔).常規(guī)觀測(cè)系統(tǒng)采集的海洋二維地震資料，在進(jìn)行簡(jiǎn)單的數(shù)據(jù)規(guī)則化后，上述假設(shè)條件就基本得到滿足，所以在二維地震資料多次波壓制處理中，SRME方法得到成功應(yīng)用.

從原理上來(lái)講，2D SRME方法可直接擴(kuò)展到3D，而且三維 SRME方法避免了二維算法的所有假設(shè)(Matson and Abma，2005；Dragoset et al.，2010).然而，當(dāng)前海洋三維地震資料采集的縱測(cè)線稀少，因此橫測(cè)線方向的采樣間隔較大，這就存在計(jì)算成本高、存儲(chǔ)以及空間采樣不足等問(wèn)題，致使全三維SRME技術(shù)一直未能得到有效應(yīng)用.

Verschuur(1991)和Ross 等(1999)等將三維地震資料按照二維方法進(jìn)行多次波壓制預(yù)處理，也就是說(shuō)對(duì)不同的檢波線依次進(jìn)行處理，在完成簡(jiǎn)單的數(shù)據(jù)規(guī)則化后，可完成依據(jù)波動(dòng)理論的多次波壓制處理.很顯然，該方法不能充分考慮三維介質(zhì)的傾角和方位角等信息，對(duì)于地層構(gòu)造簡(jiǎn)單的地區(qū)可以應(yīng)用，但是對(duì)于稍微復(fù)雜的地層，該方法在多次波預(yù)測(cè)的走時(shí)、相位和振幅等諸多方面都存在很大誤差，難以取得理想的多次波預(yù)處理結(jié)果(van Dedem and Verschuur，2005).因此，進(jìn)行全三維SRME實(shí)用方法研究是地震勘探的迫切需求.van Dedem和 Verschuur(2001)、Hokstad和Sollie(2003)提出貢獻(xiàn)道集概念，為3D SRME的有效實(shí)現(xiàn)探索新的思路.本文借鑒貢獻(xiàn)道集的思想，利用基于稀疏反演方法實(shí)現(xiàn)3D SRME表面多次波壓制.

依據(jù)實(shí)際海洋地震資料采集的特點(diǎn)計(jì)算三維地震資料表面多次波的貢獻(xiàn)道集.先對(duì)縱測(cè)線進(jìn)行道集記錄的簡(jiǎn)單外推和規(guī)則化后，形成可按照二維方法求取貢獻(xiàn)道集的地震數(shù)據(jù)，橫測(cè)線方向數(shù)據(jù)很稀疏，但是橫測(cè)線方向地震同相軸可視為雙曲線或拋物線，因此可將該同相軸沿給定的曲線路徑進(jìn)行積分，進(jìn)而可進(jìn)行反演計(jì)算，以代替橫測(cè)線方向的求和過(guò)程，降低3D SRME對(duì)運(yùn)算數(shù)據(jù)的要求，避免數(shù)據(jù)重建的處理，減少貢獻(xiàn)道集菲涅爾帶求和過(guò)程中的空間假頻.由于拋物Radon變換求解矩陣小，頻率域解耦和計(jì)算效率高等優(yōu)點(diǎn)，本文針對(duì)三維地震數(shù)據(jù)橫測(cè)線稀疏的情況，并依據(jù)其拋物線型同相軸的假設(shè)，進(jìn)行稀疏反演三維表面多次波壓制方法研究，并對(duì)理論和實(shí)際模型進(jìn)行了試算，顯示了多次波壓制方法的效果和實(shí)用性.

2 三維SRME表面多次波壓制基本原理

Verschuur(1992)等提出了壓制表面多次波的2D SRME級(jí)數(shù)展開法，即：

P(ω)=D(ω)-A(ω)D2(ω)+A2(ω)D3(ω)-…，

(1)

其中，P(ω)為只含有效波的地震數(shù)據(jù)，D(ω)為含多次波的地震數(shù)據(jù)，A(ω)為與表面反射系數(shù)r0和震源子波S(ω)相關(guān)的表面算子，公式為

(2)

其中，r0為表面反射系數(shù)，S-1(ω)為逆源子波.

從式(1)和(2)可知，2DSRME方法依賴于原始地震數(shù)據(jù)、表面反射系數(shù)以及震源子波信息.但是實(shí)際上震源子波未知，Verschuur等(1992)假設(shè)多次波壓制后地震數(shù)據(jù)能量最小，利用最小二乘方法實(shí)現(xiàn)了多次波的自適應(yīng)衰減.

理論上講，級(jí)數(shù)展開法要求預(yù)測(cè)所有階表面多次波，而且要求數(shù)據(jù)矩陣的多次相乘計(jì)算.基于多次波形成的物理原理，Berkhout和Verschuur(1997)等給出了多次波壓制的SRME迭代計(jì)算法，迭代形式如下：P(i+1)(ω)=D(ω)-A(ω)P(i)(ω)D(ω)，

(3)

其中，i代表迭代次數(shù).一般來(lái)說(shuō)，首次迭代時(shí)，取P(1)(ω)=D(ω)，類似于級(jí)數(shù)展開法，通過(guò)預(yù)測(cè)多次波與總波場(chǎng)的最小二乘匹配可估計(jì)出震源子波，進(jìn)而實(shí)現(xiàn)多次波壓制處理.

在公式(3)中，定義IR(ω)代表第i次迭代的多次波壓制結(jié)果，即IR(ω)=P(i)(ω)，定義M(ω)為預(yù)測(cè)多次波.則多次波預(yù)測(cè)公式可表述為

(4)

對(duì)任一頻率，依據(jù)式(4)預(yù)測(cè)多次波需要計(jì)算矩陣乘法，也就是對(duì)矩陣IR(ω)第i行和矩陣D(ω)的第j列進(jìn)行內(nèi)積計(jì)算.

由公式(4)，單道二維地震數(shù)據(jù)的多次波預(yù)測(cè)可進(jìn)一步表示為

(5)

將二維多次波預(yù)測(cè)方法直接擴(kuò)展到三維，得到3DSRME多次波預(yù)測(cè)公式為

×D(xk,yk,ω;xs,ys)，

(6)

其中，xr、yr、xs、ys分別表示檢波點(diǎn)與炮點(diǎn)的空間坐標(biāo)，xk、yk分別表示求和空間范圍內(nèi)震源和檢波點(diǎn)位置的坐標(biāo)變量.

依據(jù)貢獻(xiàn)道集的概念，式(6)中的求和可進(jìn)一步用貢獻(xiàn)道集作以說(shuō)明，若令Mxy為多次波貢獻(xiàn)道集，則有

Mxy(xr,yr,xs,ys,ω;xk,yk)=IR(xr,yr,ω;xk,yk)

×D(xk,yk,ω;xs,ys)，

(7)

上式也可寫成褶積表達(dá)式為

mxy(xr,yr,xs,ys,t;xk,yk)=ir(xr,yr,t;xk,yk)

*d(xk,yk,t;xs,ys)，

(8)

因此，應(yīng)用貢獻(xiàn)道集的思想，地震數(shù)據(jù)多次波預(yù)測(cè)的數(shù)學(xué)描述為

(9)

其中，mxy為多次波貢獻(xiàn)道集的時(shí)空域表示.不難發(fā)現(xiàn)，SRME方法壓制多次波可分為三個(gè)步驟：即貢獻(xiàn)道集計(jì)算、貢獻(xiàn)道集疊加和多次波自適應(yīng)相減.

設(shè)計(jì)如圖1a所示的三維地質(zhì)模型，測(cè)試分析二維和三維SRME方法的多次波預(yù)測(cè)和壓制效果.模型最上層為海水，下面為含斷層、陡傾角和尖滅等復(fù)雜構(gòu)造的地層.深度域速度模型三個(gè)空間方向的網(wǎng)格數(shù)均為400，網(wǎng)格間距為5m.設(shè)計(jì)Inline方向和Crossline方向均100個(gè)檢波器接收的觀測(cè)系統(tǒng)，也就是道間距和炮間距均為20m，震源子波為主頻25Hz的零相位Ricker子波，時(shí)間方向采樣間隔為4ms，記錄道長(zhǎng)為2s.按傳統(tǒng)三維SRME多次波預(yù)測(cè)對(duì)地震數(shù)據(jù)的要求，在每個(gè)檢波點(diǎn)都布設(shè)了炮點(diǎn)，共模擬10000炮地震記錄.圖1b為正演的三維某單炮地震記錄，可見(jiàn)能量較強(qiáng)的表面多次波.

文中應(yīng)用公式(7)對(duì)圖1所示的復(fù)雜模型數(shù)據(jù)進(jìn)行三維多次波的預(yù)測(cè)與壓制試算，同時(shí)也應(yīng)用公式(5)對(duì)抽取三維數(shù)據(jù)中一條二維測(cè)線(Inline50)進(jìn)行二維SRME的計(jì)算.2D和3DSRME算法預(yù)測(cè)和壓制多次波效果如圖2所示.圖2a為正演模擬的含多次波Inline50線共零偏移距剖面，從圖中可以看出含有能量較強(qiáng)的多次波.圖2b和2c分別為2DSRME和3DSRME算法的多次波預(yù)測(cè)結(jié)果，對(duì)比分析可知，3DSRME算法預(yù)測(cè)結(jié)果的同相軸更加連續(xù)，特別需要說(shuō)明的是圖中箭頭所示位置處，適應(yīng)三維復(fù)雜介質(zhì)的全三維SRME算法能夠清晰預(yù)測(cè)出多次波.如前所述，當(dāng)?shù)貙觾A角較大或者較為復(fù)雜時(shí)，2DSRME對(duì)多次波預(yù)測(cè)的各種假設(shè)已經(jīng)不符合實(shí)際地質(zhì)條件，也就是說(shuō)對(duì)于上述地震數(shù)據(jù)體，2DSRME算法預(yù)測(cè)的多次波出現(xiàn)較為明顯的走時(shí)誤差.3DSRME應(yīng)用全三維的預(yù)測(cè)算法，其預(yù)測(cè)的多次波位置符合實(shí)際數(shù)據(jù)情況，因此預(yù)測(cè)復(fù)雜構(gòu)造數(shù)據(jù)多次波時(shí)，3DSRME算法精度要明顯高于2DSRME算法.圖2d和2e分別為采用相同自適應(yīng)相減方法(石穎等，2013)得到的2DSRME和3DSRME多次波壓制結(jié)果，由上兩圖中箭頭所示位置可知，利用2DSRME方法多次波壓制不徹底，多次波能量有較大殘余，相比而言，3DSRME方法壓制多次波的效果較為理想.

圖1 三維復(fù)雜模型與正演模擬炮記錄 (a) 三維速度模型； (b) 正演模擬三維單炮記錄.Fig.1 Three dimensional complex geological model and modeling shot record (a) 3D velocity model; (b) The 3D shot record modeled by wave-equation method.

3 基于稀疏反演三維表面多次波壓制

對(duì)圖1所示的理論模型，設(shè)計(jì)觀測(cè)系統(tǒng)的炮間隔與檢波點(diǎn)間隔相等，因此模擬數(shù)據(jù)為全波場(chǎng)數(shù)據(jù)，符合SRME方法預(yù)測(cè)表面多次波的數(shù)據(jù)假設(shè)條件.對(duì)實(shí)際地震數(shù)據(jù)來(lái)說(shuō)，其Inline方向采樣較密，稍做預(yù)處理后，可滿足多次波預(yù)測(cè)計(jì)算的要求，但是由于觀測(cè)系統(tǒng)的限制，其Crossline方向通常較為稀疏，如果直接按照公式(9)所述方法預(yù)測(cè)多次波，將會(huì)產(chǎn)生嚴(yán)重的假頻，通常采用數(shù)據(jù)插值和規(guī)則化的方法，獲得全波場(chǎng)數(shù)據(jù)，再利用SRME方法計(jì)算，但是數(shù)據(jù)插值與規(guī)則化后，將產(chǎn)生更大的地震數(shù)據(jù)體，因此計(jì)算效率和存儲(chǔ)成本是制約SRME方法有效應(yīng)用的兩個(gè)重要問(wèn)題.不論是二維還是三維SRME方法，多次波預(yù)測(cè)的計(jì)算效率可通過(guò)應(yīng)用CPU/GPU異構(gòu)平臺(tái)加速計(jì)算來(lái)解決(石穎等，2013)，對(duì)海量地震數(shù)據(jù)的規(guī)則化預(yù)處理而言，波場(chǎng)存儲(chǔ)則具有很大的挑戰(zhàn)性.

由公式(9)可知，多次波預(yù)測(cè)由貢獻(xiàn)道集計(jì)算與疊加兩步組成.對(duì)于實(shí)際的三維地震資料而言，預(yù)測(cè)單道多次波的貢獻(xiàn)道集為三維數(shù)據(jù)體，可將貢獻(xiàn)道集計(jì)算分為兩個(gè)步驟.為此，定義如下數(shù)據(jù)向量：

(10)

公式(10)表示沿主線方向進(jìn)行疊加求和后得到的數(shù)據(jù)向量.理論上，再按照Crossline方向進(jìn)行疊加求和就得到最終預(yù)測(cè)的多次波.實(shí)際上，Crossline方向數(shù)據(jù)通常較為稀疏，無(wú)法進(jìn)行求和運(yùn)算.vanDedem和Verschuur(2001)證明公式(10)表示的數(shù)據(jù)的時(shí)距曲線可近似表示為如下雙曲線形式：

(11)

其中，τ為截距時(shí)間，ya為雙曲線頂點(diǎn)位置坐標(biāo)，v為速度，滿足v=vnmo/2，vnmo為正常時(shí)差校正速度.

將式(11)進(jìn)行Taylor展開，可得到拋物線方程為

(12)

其中，q=1/2τv2為曲率參數(shù)，公式(12)的假設(shè)條件是勘探目標(biāo)層的埋藏深度遠(yuǎn)大于Crossline方向的寬度，對(duì)多數(shù)實(shí)際三維地震數(shù)據(jù)，上述假設(shè)是可以滿足的.

圖2 2D和3D SRME預(yù)測(cè)及壓制多次波對(duì)比 (a) 含多次波的原始共偏移距剖面； (b) 2D SRME方法預(yù)測(cè)的多次波； (c) 3D SRME方法預(yù)測(cè)的多次波； (d) 基于2D SRME方法預(yù)測(cè)的多次波壓制結(jié)果； (e) 基于3D SRME方法預(yù)測(cè)的多次波壓制結(jié)果.Fig.2 Comparison of multiple prediction and suppression using 2D and 3D SRME (a) The modeled common offset section; (b) Predicted surface-related multiple by 2D SRME; (c) Predicted surface-related multiple by 3D SRME; (d) Multiple suppression based on 2D SRME multiple prediction; (e) Multiple suppression based on 3D SRME multiple prediction.

對(duì)于頻率域的數(shù)據(jù)，公式(12)的拋物Radon變換形式可表示為

(13)

其中，qi(i=1,2,…,Nq)是曲率參數(shù)，yaj(j=1,2,…,Nya)為拋物線頂點(diǎn)位置，上述方程用矩陣形式表示為

M=LHD，

(14)

其中，L為算子矩陣，上角標(biāo)H代表復(fù)共軛轉(zhuǎn)置，D為數(shù)據(jù)矩陣.其反變換可表示為

(15)

公式(14)模型空間的估計(jì)可通過(guò)最小二乘反演問(wèn)題的計(jì)算來(lái)實(shí)現(xiàn).假設(shè)模型空間是稀疏的，在模型求解時(shí)應(yīng)用Cauchy概率密度函數(shù)和Gauss-Cauchy準(zhǔn)則，可得到稀疏的高分辨率的解(SacchiandUlrych,1995;Tradetal.,2003;HokstadandSollie, 2003).模型空間的最大后驗(yàn)估計(jì)，也就是其解的形式為

(16)

公式(16)為用于求解超定方程的加權(quán)最小二乘解(Wang, 2003).但本文求解的是數(shù)據(jù)空間點(diǎn)數(shù)小于模型空間的情況，也就是所謂的欠定問(wèn)題，根據(jù)Sacchi和Ulrych(1995)以及Hokstad和Sollie(2003)等給出的矩陣形式，可將公式(16)寫為

(17)

(18)

式中，mj為模型空間矩陣的第j個(gè)分量.由式(18)可知，式(17)所表示的為模型空間的非線性稀疏解，該方程可應(yīng)用迭代法進(jìn)行求解(SacchiandUlrych,1995;HokstadandSollie,2003).

應(yīng)用非線性拋物稀疏反演替代多次波貢獻(xiàn)道集稀疏的Crossline方向的疊加，可大幅減少三維數(shù)據(jù)的疊前規(guī)則化預(yù)處理.同時(shí)這里需要指出的是在得到預(yù)測(cè)的多次波數(shù)據(jù)時(shí)，不需進(jìn)行模型空間數(shù)據(jù)的反變換和求和運(yùn)算，而直接對(duì)稀疏的模型空間數(shù)據(jù)疊加即得到預(yù)測(cè)的多次波，但是這種三維多次波預(yù)測(cè)結(jié)果存在相位誤差.依據(jù)穩(wěn)相近似理論(Bleistein,1984;Wapenaar, 1992)，可應(yīng)用如下公式進(jìn)行相位校正：

(19)

模型空間數(shù)據(jù)求和后，應(yīng)用公式(19)進(jìn)行校正，就得到最終三維表面多次波預(yù)測(cè)數(shù)據(jù)為

(20)

上述基于稀疏反演的三維表面多次波預(yù)測(cè)主要包括四個(gè)步驟，即：沿Inline方向貢獻(xiàn)道集求和，Crossline方向貢獻(xiàn)道集數(shù)據(jù)的稀疏反演，反演模型空間數(shù)據(jù)相位校正以及校正后疊加，得到如式(20)所示的三維表面多次波預(yù)測(cè)結(jié)果.為了對(duì)比分析預(yù)測(cè)效果，文中也進(jìn)行了多次波的自適應(yīng)相減計(jì)算(BerkhoutandVerschuur, 1997;Abmaetal., 2005; 石穎等, 2013).不同于常規(guī)的SRME方法，基于稀疏反演的表面多次波預(yù)測(cè)的方法，可避免Crossline方向的數(shù)據(jù)規(guī)則化處理，有效緩解數(shù)據(jù)存儲(chǔ)的壓力，降低三維SRME方法對(duì)地震采集數(shù)據(jù)的要求.

4 數(shù)據(jù)測(cè)試分析

4.1 理論數(shù)據(jù)試算

對(duì)所形成的基于稀疏反演三維表面多次波預(yù)測(cè)方法，文中也應(yīng)用CPU/GPU并行加速技術(shù)提高三維數(shù)據(jù)體多次波預(yù)測(cè)和壓制的計(jì)算效率.為驗(yàn)證稀疏反演3D SRME方法，設(shè)計(jì)了含四層水平層狀介質(zhì)的理論模型，正演模擬時(shí)采用全3D地震數(shù)據(jù)采集形式，觀測(cè)系統(tǒng)包括56條炮線，在x方向和y方向分別放置56個(gè)檢波器，每條炮線等間隔放置56個(gè)震源點(diǎn)，每炮3136道接收，道間距和炮間距均為15 m，時(shí)間采樣間隔為4 ms.抽取地震數(shù)據(jù)中的5炮進(jìn)行顯示和分析(圖3)，圖3a為原始含多次波模型數(shù)據(jù)，圖3b為采用傳統(tǒng)2D SRME方法預(yù)測(cè)的多次波，圖3c為利用傳統(tǒng)3D SRME方法對(duì)全波場(chǎng)數(shù)據(jù)預(yù)測(cè)的多次波，圖3d為抽取全波場(chǎng)數(shù)據(jù)中四條測(cè)線(線間距為75 m)，利用稀疏反演方法預(yù)測(cè)的多次波.分析可知，相比于2D SRME方法，傳統(tǒng)3D SRME和基于稀疏反演的3D SRME方法改善了預(yù)測(cè)效果.此外，拋物稀疏反演的3D SRME方法避免了Crossline方向的大規(guī)模數(shù)據(jù)規(guī)則化，極大地降低了計(jì)算成本，節(jié)約了大部分存儲(chǔ)空間.

為更好理解本文給出的稀疏反演表面多次波預(yù)測(cè)原理，圖4給出了三維地震數(shù)據(jù)中單道多次波預(yù)測(cè)的對(duì)比分析.該方法可有效避免數(shù)據(jù)規(guī)則化，特別是稀疏Crossline方向數(shù)據(jù)的大規(guī)模重建處理.基于圖3抽出的四條橫測(cè)線數(shù)據(jù)(線間距為75 m)，在完成沿主線方向的貢獻(xiàn)道集計(jì)算與疊加后，圖4a給出了某一接收點(diǎn)處 Crossline貢獻(xiàn)道集剖面，在稀疏反演三維表面多次波預(yù)測(cè)計(jì)算中，該道集就是稀疏反演的輸入數(shù)據(jù)，根據(jù)公式(17)計(jì)算，稀疏反演的模型空間數(shù)據(jù)如圖4b所示.圖4c為單道多次波預(yù)測(cè)數(shù)據(jù)的對(duì)比，重點(diǎn)對(duì)比密集Crossline貢獻(xiàn)道集疊加(傳統(tǒng)3D SRME方法)和模型道集疊加(稀疏反演方法)多次波預(yù)測(cè)結(jié)果.其中第1道為正演模擬的含多次波的原始道數(shù)據(jù)；基于全波場(chǎng)數(shù)據(jù)，利用傳統(tǒng)3D SRME方法預(yù)測(cè)的多次波如圖4c中第2道所示；第3道為輸入的數(shù)據(jù)Inline網(wǎng)格不變，但只有4條主線(橫測(cè)線方向稀疏)，采用本文稀疏反演三維SRME方法預(yù)測(cè)的多次波.與第2道對(duì)比可知，稀疏反演方法預(yù)測(cè)多次波在能量級(jí)別與走時(shí)上和傳統(tǒng)3D SRME方法基本一致，但是相位出現(xiàn)了反轉(zhuǎn)現(xiàn)象，需要應(yīng)用公式(20)進(jìn)行穩(wěn)態(tài)相位校正，圖4c第4道為相位校正后的結(jié)果，和第2道對(duì)比可知，稀疏反演多次波預(yù)測(cè)結(jié)果經(jīng)相位校正后和傳統(tǒng)全數(shù)據(jù)SRME方法預(yù)測(cè)得到的多次波無(wú)明顯差別，相位校正后的數(shù)據(jù)體就是多次波自適應(yīng)減的輸入數(shù)據(jù).

圖3 表面多次波預(yù)測(cè)與對(duì)比 (a) 原始炮數(shù)據(jù)； (b) 2D SRME方法預(yù)測(cè)； (c) 傳統(tǒng)3D SRME方法預(yù)測(cè)； (d) 拋物稀疏反演3D SRME方法預(yù)測(cè).Fig.3 Surface-related multiple prediction comparison for different prediction approaches (a) Shot records with multiple； (b) Predicted multiple shot gathers using 2D SRME； (c) Predicted multiple shot gathers using traditional 3D SRME； (d) Predicted multiple using the proposed sparse inversion 3D SRME prediction algorithm.

4.2 實(shí)際數(shù)據(jù)計(jì)算

在理論模型驗(yàn)證的基礎(chǔ)上，本文也將研究形成的稀疏反演三維表面多次波預(yù)測(cè)與壓制算法應(yīng)用于中國(guó)某海洋探區(qū)的實(shí)際地震資料.該數(shù)據(jù)有6條Inline線(拖纜)，每條拖纜有144個(gè)檢波器，線間距為100 m，道間距為12.5 m，炮間距為50 m，炮線和檢波線平行，時(shí)間方向采樣間隔為4 ms.根據(jù)稀疏反演多次波預(yù)測(cè)對(duì)地震數(shù)據(jù)的要求，對(duì)Inline線方向數(shù)據(jù)進(jìn)行了規(guī)則化處理，規(guī)則化后炮間距也為12.5 m，計(jì)算所用炮數(shù)為204炮.規(guī)則化后數(shù)據(jù)Inline測(cè)線方向道間距為12.5 m，而Crossline方向的道間距為100 m，由于橫測(cè)線方向間距過(guò)大，傳統(tǒng)的3D SRME方法是無(wú)法預(yù)測(cè)準(zhǔn)確的多次波數(shù)據(jù)，而且會(huì)帶來(lái)嚴(yán)重的假頻干擾.

應(yīng)用稀疏反演三維表面多次波壓制方法對(duì)主線規(guī)則化后的海洋地震數(shù)據(jù)進(jìn)行了多次波預(yù)測(cè)和壓制

圖4 稀疏反演三維SRME方法預(yù)測(cè)多次波單道數(shù)據(jù)分析 (a) Crossline多次波貢獻(xiàn)道集(作為稀疏反演的輸入數(shù)據(jù))； (b) 稀疏反演的模型空間域剖面； (c) 原始單道數(shù)據(jù)(1)、密集道集疊加(2)、模型道疊加(3)與相位校正道(4)對(duì)比.Fig.4 Single trace analysis of multiple prediction by three-dimensional SRME using sparse inversion (a) The Crossline multiple contributions used as input for the sparse inverse processing; (b) Sparse inversion data in model space domain; (c) The single trace of original data (trace1), stack of dense gathers result (trace2), stack of model gathers result(trace3), and the phase correction trace (trace4).

圖5 三維實(shí)際海洋數(shù)據(jù)表面多次波壓制 (a) 原始地震數(shù)據(jù)共偏移距剖面； (b) 2D SRME方法預(yù)測(cè)多次波； (c) 稀疏反演3D SRME方法預(yù)測(cè)多次波； (d) 2D SRME方法多次波壓制； (e) 稀疏反演3D SRME方法多次波壓制.Fig.5 3D marine dataset for multiple suppression (a) The common offset section of original seismic data; (b) The multiple predicted by 2D SRME; (c) The multiple predicted by sparse inversion 3D SRME; (d) The demultiple section by 2D SRME; (e) The demultiple section by sparse inversion 3D SRME.

處理，為對(duì)比分析計(jì)算效果，文中也應(yīng)用二維SRME方法對(duì)該實(shí)際數(shù)據(jù)進(jìn)行了處理.圖5a為三維海洋數(shù)據(jù)體中的某條Inline線的含多次波的原始共偏移距剖面，可以看出該數(shù)據(jù)含有較強(qiáng)的海底一階、二階甚至高階的表面多次波.圖5b為2D SRME方法預(yù)測(cè)的表面多次波，圖5c為基于稀疏反演三維方法預(yù)測(cè)的表面多次波.對(duì)比可知，三維稀疏反演方法預(yù)測(cè)的表面多次波的同相軸連續(xù)性更好，與原始多次波數(shù)據(jù)在能量對(duì)比關(guān)系上具有更好的一致性，也就是說(shuō)可更好的與原始剖面的多次波同相軸進(jìn)行擬合，為多次波自適應(yīng)相減的有效處理提供了數(shù)據(jù)基礎(chǔ).

針對(duì)兩種方法的預(yù)測(cè)結(jié)果，文中采用前述理論模型測(cè)試的多次波壓制算法進(jìn)行處理，2D SRME及稀疏反演3D SRME多次波壓制結(jié)果分別如圖5d和5e所示.對(duì)圖中所示的雙程走時(shí)約為0.55 s位置處箭頭所示同相軸分析表明，2D SRME方法壓制多次波后，其海底一階表面多次波同相軸的能量仍很強(qiáng)，而3D稀疏反演的SRME方法中多次波能量得到很好的壓制，同樣，在雙程走時(shí)約為1.5 s和1.8 s附近的黑色箭頭也顯示3D稀疏反演SRME方法多次波壓制效果明顯優(yōu)于2D算法.雙程走時(shí)約為0.8 s處箭頭所示同相軸為一有效波同相軸，這一點(diǎn)從原始共偏移距剖面和預(yù)測(cè)的多次波剖面對(duì)比可以判斷出，由于2D SRME預(yù)測(cè)的多次波存在振幅和相位的誤差，使得自適應(yīng)多次波匹配相減無(wú)法有效實(shí)現(xiàn)，有效波的振幅受到一定程度的損傷，與2D SRME方法相比，基于全三維假設(shè)的3D稀疏反演的SRME方法可較為精確的預(yù)測(cè)出多次波走時(shí)，進(jìn)而可與原始數(shù)據(jù)中表面多次波有效擬合，可實(shí)現(xiàn)保幅的表面多次波壓制處理.

5 結(jié)論

相比于二維方法壓制三維地震數(shù)據(jù)的多次波，三維SRME方法考慮了地震波在地下介質(zhì)中傳播機(jī)理和實(shí)際情況，更為適合復(fù)雜構(gòu)造海洋地震資料的多次波壓制處理，但是該算法對(duì)數(shù)據(jù)的采集要求很高.實(shí)際采集地震數(shù)據(jù)的橫測(cè)線方向通常較為稀疏，為提高計(jì)算效率、避免大數(shù)據(jù)體的規(guī)則化處理以及規(guī)則化數(shù)據(jù)的存儲(chǔ)，本文基于貢獻(xiàn)道集的概念，將三維表面多次波貢獻(xiàn)道集求取分兩步實(shí)現(xiàn)，也就是在完成縱測(cè)線方向的多次波貢獻(xiàn)積分后，應(yīng)用高分辨率拋物Radon變換稀疏反演代替橫測(cè)線的多次波道集求和，進(jìn)而完成三維表面多次波的預(yù)測(cè).應(yīng)用CPU/GPU并行加速技術(shù)提高三維地震資料多次波預(yù)測(cè)和壓制的計(jì)算效率，可大幅縮短海洋地震資料預(yù)處理周期.對(duì)于稀疏反演的三維多次波預(yù)測(cè)數(shù)據(jù)和全三維全數(shù)據(jù)預(yù)測(cè)的多次波存在的相位差異，文中應(yīng)用穩(wěn)相原理近似完成預(yù)測(cè)多次波的相位校正.理論模擬和實(shí)際海洋勘探三維數(shù)據(jù)的試算表明，稀疏反演三維表面多次波壓制方法具有計(jì)算精度高、適用范圍廣和易于實(shí)現(xiàn)的特點(diǎn).

Abma R, Kabir N, Matson K H, et al. 2005. Comparisons of adaptive subtraction methods for multiple attenuation.TheLeadingEdge, 24(3): 277-280.

Berkhout A J. 1980. Seismic Migration: Imaging of Acoustic Energy by Wave Field Extrapolation. Part A: Theoretical Aspects. Elsevier.

Berkhout A J, Verschuur D J. 1997. Estimation of multiple scattering by iterative inversion, part I: Theoretical considerations.Geophysics, 62(5): 1586-1595.Bisley R, MacNeill M. 2008. Application of true-azimuth 3D SRME in the Northwest Shelf, Australia. ∥78th Annual International Meeting, SEG, Expanded Abstracts, 2436-2440.Bleistein N. 1984. Mathematical Methods for Wave Phenomena. Academic Press.

Dragoset B, Moore I, Kostov C. 2006. The impact of field-survey characteristics on surface-related multiple attenuation.GeophysicalProspecting, 54(6): 781-791. Dragoset B, Verschuur D J, Moore I, et al. 2010. A perspective on 3D surface-related multiple elimination.Geophysics, 75(5): 75A245-75A261.

Hokstad K, Sollie R. 2003. 3D surface-related multiple elimination using parabolic sparse inversion. ∥73rd Annual International Meeting, SEG, Expanded Abstracts, 1961-1964.

Jakubowicz H. 1998. Wave equation prediction and removal of interbed multiples. ∥68th Annual International Meeting, SEG, Expanded Abstracts, 1527-1530.

Kabir M N, Verschuur D J. 1995. Restoration of missing offsets by parabolic Radon transform.GeophysicalProspecting, 43(3): 347-368.

Kelamis P G, Verschuur D J. 2000. Surface-related multiple elimination on land seismic data—Strategies via case studies.Geophysics, 65(3): 719-734.

Li X, Hu T Y. 2009. Surface-related multiple removal with inverse scattering series method.ChineseJ.Geophy. (in Chinese), 52(6): 1633-1640, doi: 10.3969/j.issn.0001-5733.2009.06.026.Matson K H, Abma R. 2005. Fast 3D surface-related multiple elimination using azimuth moveout for multiples. ∥75th Annual International Meeting, SEG, Expanded Abstracts, 2064-2067.

Ross W S, Yu Y, Gasparotto F A. 1999. Travel time prediction and suppression of 3-D multiples.Geophysics, 64(1): 261-277.

Sacchi M D, Ulrych T J. 1995. High-resolution velocity gathers and offset space reconstruction.Geopyhsics, 60(4): 1169-1177.

Shi Y, Wang W H, Li Y, et al. 2013. 3D surface-related multiple prediction approach investigation based on wave equation.ChineseJ.Geophys. (in Chinese), 56(6): 2023-2032, doi: 10.6038/cjg20130623.

Sun Y. 1999. Anti-aliasing multiple prediction beyond 2-D. ∥69th Annual International Meeting SEG, Expanded Abstracts, 1338-1341

Trad D, Ulrych T, Sacchi M. 2003. Latest views of the sparse Radon transform.Geophysics, 68(1): 386-399.

van Dedem E J, Verschuur D J. 2001. 3D surface multiple prediction using sparse inversion. ∥71st Annual International Meeting, SEG, Expanded Abstracts, 1285-1288.

van Dedem E J, Verschuur D J. 2005. 3D surface-related multiple prediction: A sparse inversion approach.Geophysics, 70(3): V31-V43.

Verschuur D J. 1991. Surface-related multiple elimination, an inversion approach[Ph. D. thesis]. Netherlands: Delft University of Technology.Verschuur D J, Berkhout A J. 1997. Estimation of multiple scattering by iterative inversion, Part II: Practical aspects and examples.Geophysics, 62(5): 1596-1611.

Verschuur D J, Berkhout A J, Wapenaar C P A. 1992. Adaptive surface-related multiple elimination.Geophysics, 57(9): 1166-1177.

Wang B L, Sacchi M D, Yin X Y, et al. 2014. Multiple attenuation

based on amplitude preserving Radon transform.ChineseJ.Geophys. (in Chinese), 57(6): 1924-1933, doi: 10.6038/cjg20140623.

Wang Y H. 2003. Sparseness-constrained least-squares inversion: Application to seismic wave reconstruction.Geophysics, 68(5): 1633-1638.

Wapenaar C P A. 1992. The stationary phase method and its application in seismics and acoustics. Lecture notes: Akoestische Beeldvorming, Delft University of Technology.

附中文參考文獻(xiàn)

李翔, 胡天躍. 2009. 逆散射級(jí)數(shù)法去除自由表面多次波. 地球物理學(xué)報(bào), 52(6): 1633-1640, doi: 10.3969/j.issn.0001-5733.2009.06.026.

石穎, 王維紅, 李瑩等. 2013. 基于波動(dòng)方程三維表面多次波預(yù)測(cè)方法研究. 地球物理學(xué)報(bào), 56(6): 2023-2032, doi: 10.6038/cjg20130623.

王保麗, Sacchi M D, 印興耀等. 2014. 基于保幅拉東變換的多次波衰減. 地球物理學(xué)報(bào), 57(6): 1924-1933, doi: 10.6038/cjg20140623.

(本文編輯 張正峰)

3D surface-related multiple elimination based on sparse inversion

WANG Wei-Hong, JING Hong-Liang

CollegeofEarthSciences,NortheastPetroleumUniversity,Daqing163318,China

Three dimensional surface-related multiple elimination (SRME) is one of the important topics in the processing of seismic data from marine exploration,theoretically, the data-driven SRME based on wave-equation, can suppress all surface-related multiples from complex structure, both in 2D and 3D sense. But actually, because of the high requirement for seismic data acquisition, it is usually difficult to apply 3D SRME for field data demultiple processing.3D multiple suppression approach based on sparse inversion is analyzed.We classify multiple suppression method into two categories, filter and SRME method respectively. For the seismic data from complex geological structure, filter approach doesn′t work well. However, the data-driven SRME approach based on wave equation, can suppress multiple better, which has no requirement for velocity information.For SRME, the full wavefield data requirement is an important disadvantage, which cannot be meet for almost all marine field data. Therefore the data reconstruction is necessary for traditional multiple suppression using SRME.For current marine acquisition geometries, the data is densely sampled in the inline direction, but very sparsely in the crossline direction, we introduce contribution gather concept, and calculate sparsely sampled crossline multiple contribution by means of sparse inversion algorithm.Therefore, the data reconstruction is unnecessary before demultiple,comparing with traditional 3D SRME algorithm, which decrease storage cost greatly. Based on the assumption that the crossline time-distance curves are hyperbolic or parabolic, i.e., after integrated along inline direction after finished the first step processing of contribution gathers, we applied the phase correction algorithm which based on the principle of stationary phase approximation, predicted multiples by the proposed method and result by the full 3D SRME method are basically the same on kinematics and dynamics characteristics.Multiple prediction and adaptive subtraction are two crucial steps for 3D multiple suppression using SRME. We simulate complex model data to test the proposed 3D multiple prediction algorithm, and 10 thousands shot records are modelled, the shot interval and trace interval are both 25 m. The result comparison show the 3D multiple prediction approach can predict the multiple′s amplitude and phase correctly, and also, the subtraction result is superior than 2D algorithm. The horizontal four-layered media is also designed to test the sparse inversion 3D multiple suppression algorithm, there are 3136 shot records in total, the trace and shot intervals are both 25 m, the line interval is 75 m. The single trace and common-offset result show that the proposed approach can predict the multiple′s amplitude and traveltime correctly, and they are very close to the full data circumstance. The test on field data from some area in China show that the proposed sparse inversion method is applicable and effective, where the trace interval is 12.5 m, shot interval is 50 m, and the line interval is 100 m. After summation along inline direction, the partial integration data is transformed to Radon domain using apex-shifted Radon approach based on the assumption of hyperbolic or parabolic events, stacked the Radon imaging, and also applied the phase correction to the sparse inversion solution, the predicted multiple is acquired.After theoretical investigation and data tests, we have the following conclusions, (1) the proposed method is suitable for simple and complex 3D model data, (2) for real seismic data, the inline direction reconstruction is needed, (3) because only forward Radon transform, not inverse Radon transform is used,phase correction is demanded, and (4) the algorithm do not rely on the assumption of full data using sparse inversion approach, and also, widen the application extent of 3D multiple suppression.The result show that the proposed 3D multiple suppression algorithm can improve the S/N ratio during the course of preprocessing, and provide the high quality data for the subsequent high resolution imaging.

Sparse inversion; Surface-related multiple; Contribution gather; Multiple prediction

10.6038/cjg20150725.

國(guó)家自然科學(xué)基金項(xiàng)目(41474118),國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863)(2012AA061202)和大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(LP1509)聯(lián)合資助.

王維紅，男，1975年生，博士，副教授，主要從事地震資料數(shù)字處理方面的研究. E-mail: wwhsy@sina.com

10.6038/cjg20150725

P631

2015-04-02，2015-07-06收修定稿

王維紅， 井洪亮. 2015. 基于稀疏反演三維表面多次波壓制方法.地球物理學(xué)報(bào)，58(7):2496-2507,

Wang W H, Jing H L. 2015. 3D surface-related multiple elimination based on sparse inversion.ChineseJ.Geophys. (in Chinese),58(7):2496-2507,doi:10.6038/cjg20150725.