基于Marchenko理論壓制自由表面多次波方法研究

王小衛(wèi)， 王孝， 謝俊法， 曾華會(huì)， 趙玉合

中國石油勘探開發(fā)研究院西北分院， 蘭州 730020

0 引言

對于深海地震數(shù)據(jù)，自由表面多次波衰減是地震資料處理過程中的重點(diǎn)和難點(diǎn).多次波的存在，特別是與水層相關(guān)多次波的存在常常會(huì)對接收到的有效波產(chǎn)生破壞，對地下地質(zhì)構(gòu)造成像、速度分析的準(zhǔn)確性等環(huán)節(jié)造成影響(牛濱華等, 2002；路鵬飛等，2020；王常波等，2020).

預(yù)測反褶積方法是早期地震資料處理中壓制短周期多次波和海上鳴震的有效處理方法(Robinson, 1957)，但是對于地下介質(zhì)構(gòu)造較為復(fù)雜的地區(qū)，預(yù)測反褶積對海上鳴震等短周期多次波壓制效果較差.Kennet(1979)提出了反饋迭代法，Berkhout(1982)進(jìn)一步對反饋迭代法進(jìn)行了方法改進(jìn)，使其適用于任何表面多次波模型.表面相關(guān)多次波去除方法(SRME)方法基于反饋迭代模型，在無需已知地下介質(zhì)信息前提下利用地震數(shù)據(jù)本身預(yù)測多次波，是當(dāng)前工業(yè)上成熟應(yīng)用的多次波壓制技術(shù)(Verschuur, 1991; Verschuur et al., 1992; Berkhout and Verschuur, 1997; 王維紅等, 2007; 李學(xué)聰?shù)? 2010; 石穎和王維紅, 2012; 石穎等, 2013; 徐嘉亮等，2018).對于近炮檢距缺失的地震資料，使用SRME方法需要對缺失數(shù)據(jù)進(jìn)行波場外推來獲得，數(shù)據(jù)重構(gòu)過程中將會(huì)產(chǎn)生一定誤差，影響表面多次波的壓制結(jié)果.Moore和Bisley(2006)、Martin等(2011)利用模型驅(qū)動(dòng)方法壓制自由表面多次波，基于模型的自由表面多次波壓制方法在部分短周期多次波壓制方面取得了較好的效果，但是對地震數(shù)據(jù)之外的信息依賴程度過高(趙昌壘等, 2013).基于稀疏反演的一次波估計(jì)方法EPSI(Estimation of Primaries by Sparse Inversion)是建立在與SRME相同的理論基礎(chǔ)上，通過直接估計(jì)一次波脈沖響應(yīng)和地震子波，進(jìn)而重構(gòu)出一次波和自由表面多次波，避免了自適應(yīng)相減步驟，更好地保護(hù)了有效信號(hào)(Van Groenestijn and Verschuur, 2009a,b).但是該方法是一種迭代反演方法，在迭代反演過程中存在計(jì)算量較大的問題，同時(shí)該方法計(jì)算時(shí)忽略震源的方向特性，假設(shè)所有炮具有相同的震源子波，這在實(shí)際數(shù)據(jù)中并不一定時(shí)時(shí)滿足，進(jìn)而影響預(yù)測效果.

基于Marchenko理論壓制多次波方法，主要利用地下介質(zhì)中虛源點(diǎn)的格林函數(shù)和地表得到的反射響應(yīng)進(jìn)行多次波預(yù)測.Broggini和Snieder(2012)將該方法引入地球物理領(lǐng)域.其本質(zhì)是通過虛源點(diǎn)與地表之間的直達(dá)波記錄和地表反射響應(yīng)反演得到上下行格林函數(shù).Broggini等(2012)、Wapenaar等(2013)將該方法擴(kuò)展到二維和三維.Marchenko方法輸入數(shù)據(jù)不包含自由表面多次波，因此計(jì)算得到的格林函數(shù)也不包含與自由表面相關(guān)的多次波.本文對Marchenko方程進(jìn)行了改進(jìn)，在自由表面多次波存在的情況下得到格林函數(shù)場，得到包含一次波、層間多次波和自由表面多次波的格林函數(shù).該方法能夠預(yù)測關(guān)于海水面的自由表面多次波，進(jìn)而利用自適應(yīng)匹配濾波方法使自由表面多次波得到有效壓制.

1 方法原理

Marchenko方法實(shí)質(zhì)上是通過虛源點(diǎn)與地表之間的直達(dá)波記錄和地表反射響應(yīng)反演得到格林函數(shù)的上下行波場.本文引入一個(gè)假想介質(zhì)，假想介質(zhì)中的波場如圖1a所示，實(shí)際介質(zhì)中波場如圖1b所示.圖1a中，在Di深度以下為均勻半空間的假想介質(zhì)，引入聚焦函數(shù)f1，f1在Di深度以下無反射.在實(shí)際介質(zhì)和假想介質(zhì)中，表面邊界D0以上無反射.Wapenaar等(2014)將實(shí)際介質(zhì)的上行和下行格林函數(shù)與地表的反射響應(yīng)聯(lián)系起來得到：

×R0(x″0,x0,ω)dx，

(1)

×R0(x″0,x0,ω)dx,

(2)

圖1 假想介質(zhì)中波場(a)和實(shí)際介質(zhì)中波場(b)傳播情況Fig.1 Wave field propagation in imaginary medium (a) and real medium (b)

其中，G0表示實(shí)際介質(zhì)中的格林函數(shù)，R0表示地表反射響應(yīng)，反射響應(yīng)R0不包含與自由表面相關(guān)的多次波.在本文中上標(biāo)“+”表示下行波場，上標(biāo)“-”表示上行波場.對方程(1)、(2)求和，并利用炮點(diǎn)-檢波點(diǎn)互易關(guān)系對格林函數(shù)求解：

G0(x″0,x′i,ω)=f2(x′i,x″0,ω)*

(3)

聚焦函數(shù)f2定義為：

(4)

在地震波直接到達(dá)聚焦點(diǎn)之前沒有能量到達(dá)，因此格林函數(shù)上下行波場為0，將其作為迭代方法的初始條件，將聚焦點(diǎn)作為虛源點(diǎn)，對虛源點(diǎn)直接到達(dá)檢波點(diǎn)的記錄進(jìn)行時(shí)間反轉(zhuǎn)作為聚焦函數(shù)下行波場的初始值，構(gòu)成迭代方案的基礎(chǔ)完成迭代更新.

在介紹由包含自由表面多次波情況下的反射響應(yīng)得到格林函數(shù)之前，本文先引入R、R0、G、G0的概念，R為自由表面存在時(shí)記錄的反射波，R0為無自由表面介質(zhì)的反射響應(yīng).G為存在自由表面時(shí)位于介質(zhì)中某一點(diǎn)的虛源點(diǎn)在表面上的格林函數(shù)，G0是不存在自由表面時(shí)的格林函數(shù).Marchenko方法為不存在自由表面的情況下聚焦波場的理論.也就是說，從R0中求得格林函數(shù)G0(Rose, 2002; Broggini et al., 2012; Wapenaar et al., 2013).

Wapenaar等(2004)提出將存在和不存在自由表面的介質(zhì)傳播算子聯(lián)系起來，從而通過G0計(jì)算出G(自由曲面存在時(shí)的格林函數(shù))，其頻域表達(dá)式為：

×R(x0,x″0,ω)dx，

(5)

其中，D0為采集面，x0和x′i為沿D0和Di(D0以下的任意深度平面)的空間位置，R為D0界面處下行入射波場的反射響應(yīng).在該方法中，本文將傳輸響應(yīng)替換為相應(yīng)的格林函數(shù)G或G0，因?yàn)楦窳趾瘮?shù)是從聚焦點(diǎn)到表面的總傳輸波場.需要注意的是，這種方法遵循由R進(jìn)行SRME得到R0，進(jìn)而獲得G0，最后得到G的環(huán)節(jié).然而，本文可以直接從測量的反射數(shù)據(jù)R通過計(jì)算得到自由表面存在時(shí)的格林函數(shù)G，即：直接由反射數(shù)據(jù)R求取格林函數(shù)，避免了SRME和G0的求取.本文將Marchenko方法由R0求取G0的公式推廣到包含自由表面多次波(由R求G).來自自由表面的反射包含在聚焦方法中，類似于Wapenaar等(2004)對自由表面多次波的處理.

本文定義水平和深度空間坐標(biāo)，如：x0=(xH,x3,0)、xH表示在深度x3,0處水平坐標(biāo)(x1,x2).本文定義了波動(dòng)方程在介質(zhì)中聚焦于一點(diǎn)的解，本文稱之為聚焦函數(shù)f1和f2.f1函數(shù)所涉及的波在x′i處聚焦在一個(gè)定義的深度水平(?Di)上，對于在采集表面(?D0)上x0的入射波和傳出波，f2函數(shù)是一個(gè)對于入射波和傳出波在?Di上且在D0上聚焦的波的解，如圖1a.

在?D0以上帶有自由表面的實(shí)際非均勻模型，如圖2所示.與Wapenaar等(2013)中沒有自由表面的模型不同，本文認(rèn)為來自自由表面的反射是下行源，類似于Wapenaar等(2004)的研究.在圖2中，本文描述了波場中的上行和下行分量.在?D0上的波場的下行分量為G+(x0,x″0,ω)=δ(xH-x″H)+rR(x0,x″0,ω)，等號(hào)右邊包含了下行的脈沖源和自由表面的反射.下行脈沖源δ(xH-x″H)是一個(gè)2D狄拉克δ函數(shù)，x″H為聚焦函數(shù)f2的聚焦點(diǎn)橫坐標(biāo).

圖2 實(shí)際非均勻介質(zhì)中自由表面存在時(shí)的波場Fig.2 Wave field in real inhomogeneous medium with free surface

分別將聚焦函數(shù)(f1和f2)之間的波場和聚焦函數(shù)與實(shí)際介質(zhì)之間的波場聯(lián)系起來.那么，自由表面存在時(shí)對應(yīng)的格林函數(shù)與聚焦函數(shù)f2具有關(guān)系式為：

G(x′i,x″0,ω)=f2(x′i,x″0,ω)*

(6)

對比方程(6)和方程(3)，方程(6)右邊的最后一項(xiàng)表示從自由表面反射的波場信息.方程(6)中的反射響應(yīng)R包含自由表面的反射波，而表達(dá)式(3)包含的是沒有自由表面介質(zhì)的反射響應(yīng)R0.其中，上行格林函數(shù)和下行格林函數(shù)表達(dá)式分別如式(7)和式(8)所示：

(7)

(8)

對式(7)、(8)計(jì)算得到下行格林函數(shù)進(jìn)行去直達(dá)波處理，將地下虛源點(diǎn)與地表接收點(diǎn)之間去直達(dá)波后的格林函數(shù)下行波場和完整的上行波場做褶積，預(yù)測出與自由表面相關(guān)的多次波：

(9)

2 模型測試

2.1 模型測試一

本文利用三層水平層狀速度模型如(圖3)所示驗(yàn)證本文提出方法的正確性.速度模型大小為2000 m×2500 m，水平方向有200個(gè)網(wǎng)格點(diǎn)，垂直方向上有250個(gè)網(wǎng)格點(diǎn)，網(wǎng)格間距均為10 m.共放201炮，每炮有201個(gè)檢波點(diǎn)接收，炮檢距和檢波點(diǎn)間距為10 m.所用震源子波為主頻25 Hz的雷克子波，反射時(shí)間記錄長度為3.5 s，采樣間隔為1 ms.第一個(gè)和最后一個(gè)震源位置分別位于模型起始處.構(gòu)造上下行格林函數(shù)的虛源點(diǎn)均勻分布于速度模型第一層的400 m處和第二層的1200 m處，點(diǎn)間距為10 m，如圖3中的紅點(diǎn)所示.

圖3 速度模型Fig.3 Velocity model

圖4a為正演的第100炮地震記錄，圖4b為利用Marchenko方法預(yù)測的自由表面多次波，圖4c為地震記錄壓制自由表面多次波的數(shù)據(jù).從圖4中可知Marchenko方法預(yù)測的表面多次波在走時(shí)上與實(shí)際的表面多次波有很好的一致性，振幅和頻率略有差異，該差異可由多道匹配自適應(yīng)相減方法來消除.如圖4c所示，經(jīng)過自適應(yīng)匹配相減后，正演地震記錄中自由表面多次波得到有效壓制(紅色箭頭所指)，從而說明利用本文提出的基于Marchenko自由表面多次波壓制方法可有效壓制自由表面多次波.

圖4 第100炮地震記錄(a) 多次波壓制前； (b) 預(yù)測的多次波； (c) 多次波壓制后.Fig.4 Seismic record of the 100th shot(a) Before multiple suppression; (b) Predicted multiples; (c) After multiple suppression.

2.2 模型測試二

為了驗(yàn)證該方法的有效性，本文在測試水平模型后對SMARRT模型進(jìn)行測試.SMARRT模型正演參數(shù)為：震源和檢波點(diǎn)都位于海水面，采樣時(shí)間間隔為0.008 s，記錄長度為5 s，共1387炮，每炮361道，炮間距為75 m，模型左右下界面采用吸收邊界.圖5為SMARRT速度模型，紅點(diǎn)為虛源點(diǎn)位置.

圖5 SMARRT速度模型Fig.5 SMARRT velocity model

圖6為利用Marchenko自由表面多次波壓制方法壓制多次波前后單炮對比圖，圖6a為炮號(hào)4500的原始炮集，圖6b為預(yù)測的多次波，圖6c為壓制多次波后的地震記錄.從圖中可以看到基于Marchenko自由表面多次波壓制方法可以有效壓制表面多次波(紅色箭頭所指)而無損有效信號(hào)，提高了地震資料的信噪比.

圖6 Marchenko方法多次波壓制前后單炮記錄(a) 多次波壓制前； (b) 預(yù)測的多次波； (c) 多次波壓制后.Fig.6 Single shot record before and after multiple suppression by Marchenko method(a) Before multiple suppression; (b) Predicted multiples; (c) After multiple suppression.

圖7為利用Marchenko自由表面多次波壓制方法壓制多次波前后抽取的零偏移距剖面對比圖，圖7a為原始地震剖面，圖7b為預(yù)測的多次波數(shù)據(jù)剖面，圖7c為壓制多次波后的地震剖面.從圖中可以看到自由表面多次波得到有效壓制(如紅色箭頭和矩形框內(nèi)所示).由該實(shí)例說明基于Marchenko自由表面多次波壓制方法對于復(fù)雜模型具有較強(qiáng)的適應(yīng)性，通過進(jìn)一步研究有望應(yīng)用于實(shí)際數(shù)據(jù).

圖7 Marchenko方法多次波壓制前后零偏移距剖面(a) 多次波壓制前； (b) 預(yù)測的多次波； (c) 多次波壓制后.Fig.7 Zero offset profile before and after multiple suppression by Marchenko method(a) Before multiple suppression; (b)Predicted multiples; (c) After multiple suppression.

3 結(jié)論

本文提出一種改進(jìn)的自聚焦Marchenko多次波壓制方法，該方法不需要去除海水面相關(guān)的自由表面多次波，同時(shí)能夠?qū)ψ杂杀砻娑啻尾ㄟM(jìn)行壓制.傳統(tǒng)壓制自由表面多次波方法(如SRME)需要對近道缺失的地震數(shù)據(jù)進(jìn)行補(bǔ)充，同時(shí)不能壓制地下強(qiáng)阻抗界面的層間多次波；壓制層間多次波方法也都需要先驗(yàn)地下模型，計(jì)算效率有待提高.本文提出的方法能夠改善上述方法的缺陷，不需要對近道數(shù)據(jù)進(jìn)行補(bǔ)充，也不需要地下構(gòu)造的速度模型，僅僅需要背景速度場及迭代更新策略就可以使目標(biāo)函數(shù)得到收斂.由模型測試證明本文提出的方法是有效可行的，該方法能夠?qū)ψ杂杀砻娑啻尾ㄟM(jìn)行有效壓制.

本文方法在構(gòu)建的上下行格林函數(shù)場中，同時(shí)包括自由表面多次波和層間多次波，下一步利用復(fù)雜介質(zhì)模型測試本文方法對同時(shí)壓制自由表面多次波和層間多次波的有效性，該方法有望在實(shí)際海洋數(shù)據(jù)當(dāng)中得到廣泛應(yīng)用.