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

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

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

0 引言

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

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

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

1 方法原理

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

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

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

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

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

(5)

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

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

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

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

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

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

(6)

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

(7)

(8)

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

(9)

2 模型測試

2.1 模型測試一

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

圖3 速度模型Fig.3 Velocity model

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

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

2.2 模型測試二

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

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

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

圖6 Marchenko方法多次波壓制前后單炮記錄(a) 多次波壓制前； (b) 預測的多次波； (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為預測的多次波數(shù)據(jù)剖面，圖7c為壓制多次波后的地震剖面.從圖中可以看到自由表面多次波得到有效壓制(如紅色箭頭和矩形框內(nèi)所示).由該實例說明基于Marchenko自由表面多次波壓制方法對于復雜模型具有較強的適應性，通過進一步研究有望應用于實際數(shù)據(jù).

圖7 Marchenko方法多次波壓制前后零偏移距剖面(a) 多次波壓制前； (b) 預測的多次波； (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 結論

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

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