楊旭明,王 麗,陶長江,鮑 偉

(中國石油化工股份有限公司江蘇油田分公司物探研究院,江蘇南京210046)

地震勘探中多次波問題普遍存在,多次波的存在使得一次反射波的相位、振幅和頻率等動力學(xué)特征失真,對構(gòu)造解釋和地震反演都造成了不利影響。當(dāng)一次波和多次波相互重疊時,會更加難以識別和處理多次波。因此,在地震資料處理中,有效壓制多次波一直是地震處理研究的重點和難點之一。

BERKHOUT[1]提出了描述復(fù)雜多次波系統(tǒng)的反饋理論框架,奠定了反饋迭代多次波壓制方法的數(shù)學(xué)物理基礎(chǔ),提出的多維反演算法通過數(shù)據(jù)矩陣相乘預(yù)測多次波,能夠適應(yīng)任意復(fù)雜的地下結(jié)構(gòu)。但消除多次波時需要已知震源子波。VERSCHUUR[2]引入基于一次波能量最小假設(shè)的自適應(yīng)相減方法,成功地從實際地震數(shù)據(jù)中估計出震源子波,極大地推動了表面相關(guān)多次波壓制方法(surface-related multiple elimination,SRME)的發(fā)展。為了避免因求解震源子波而產(chǎn)生的非線性優(yōu)化問題,VERSCHUUR等[3]提出了迭代SRME方法,在每次迭代過程中采用最小二乘匹配多次波模型與原始數(shù)據(jù),確定震源子波,將原來的非線性問題轉(zhuǎn)化為線性問題,增強了方法的實用性。為了避免自適應(yīng)相減可能損傷有效信號,SRME被重新定義為反演一次波脈沖響應(yīng)和震源子波的全波形反演問題,形成了基于稀疏反演的一次波估計方法(estimation of primaries by sparse inversion,EPSI)[4-5]。

傳統(tǒng)的EPSI方法需要精確地包含有一次波的時窗和反演參數(shù)設(shè)定,存在一定的不穩(wěn)定性。因此,需要對EPSI方法進行改進,將其轉(zhuǎn)化成為雙凸優(yōu)化問題[6],即反演一次波脈沖響應(yīng)過程中采用的目標(biāo)函數(shù)是一個凸函數(shù),應(yīng)用的約束集合同樣也是一個凸函數(shù),通過交替優(yōu)化迭代過程,對一次波脈沖響應(yīng)采用譜梯度投影(SPGL1)[7]算法求解,同時引入帕累托優(yōu)化(pareto)曲線。改進的EPSI方法[8-10]在求解過程中,需要在輸入記錄的多維互相關(guān)數(shù)據(jù)中拾取與最強反射一次波相一致的同相軸來估計初始的震源子波。該工作不僅費時費事,還可能由于識別錯誤和拾取參數(shù)設(shè)置誤差,得不到合適的震源子波。馮飛[11]使用常規(guī)最小平方正交三角(QR)分解算法求解震源子波。宋家文等[12]和白蘭淑等[13]通過匹配濾波方法求取地震子波。本文提出的自適應(yīng)稀疏反演多次波壓制方法(adaptive estimation of primaries by sparse inversion,AEPSI),不需要人工拾取和時窗設(shè)置,完全由數(shù)據(jù)驅(qū)動解決初始震源子波及震源子波迭代更新問題。

1 方法原理

1.1 自適應(yīng)稀疏反演多次波壓制原理

為了便于理解和推導(dǎo),EPSI頻率域方程表示為:

P=GQ+GRP

(1)

式中:P為頻率域數(shù)據(jù)矩陣;G為頻率域一次波脈沖響應(yīng);R為頻率域自由表面反射算子;Q為頻率域震源矩陣。如果用q(t)表示對應(yīng)于時間域數(shù)據(jù)矩陣p的震源特性,并假定q(t)對所有的炮和道是不變的,則Q=w(ω)I。其中,w(ω)是頻率的標(biāo)量函數(shù),是q(t)的頻率域表示;I是單位矩陣。假定自由表面反射算子R=-I。方程通過逆傅里葉變換返回到時間域,G的時間域表示為g,P的時間域表示為p,時間域子波q(t)簡記為q,則時間域數(shù)據(jù)矩陣p可表示為:

p=M(g,q;p)

(2)

由于稀疏約束是施加在時間域,所以用作用在時間域的M(g,q;p)函數(shù)討論問題會更方便。如果基于數(shù)據(jù)矩陣p的方程(2)用隱式表示,則有:

p=M(g,q)

(3)

對于方程(3)來說,如果給定變量g、q中的一個,則M(g,q)變成了相對另一個變量的線性算子。因此,M(g,q)也叫做雙線性算子,可以寫成:

(4)

(5)

(6)

根據(jù)(5)式和(6)式,定義:

(7a)

(7b)

(8)

式中:q∈Λ表示滿足任何所需的或先驗約束施加在震源上。這里假定q在時間域是短時窗的震源子波,并允許非因果的描述記錄數(shù)據(jù)中可能出現(xiàn)的時移。當(dāng)然,也可以施加其它需要的震源特性。

根據(jù)迭代求解約束優(yōu)化問題的方法,則可以通過下列兩個子問題的迭代求解EPSI約束優(yōu)化問題:

(9)

(10)

(11)

(12)

調(diào)整因子sk為:

(13)

1.2 初始震源子波的確定

實際中的地震子波是一個很復(fù)雜的問題,因為地震子波與地層巖石性質(zhì)有關(guān),地層巖石性質(zhì)本身就是一個復(fù)雜體。為了研究方便,對地震子波進行模擬,這就是理論子波。目前普遍認(rèn)為Ricker提出的地震子波數(shù)學(xué)模型具有廣泛的代表性,因此,Ricker子波可作為初始震源子波。Ricker子波公式如下:

(14)

式中:t為時間;A為振幅;fm為主頻。圖1展示了主頻從15Hz到45Hz,振幅A為1.0的Ricker子波。

圖1 不同主頻Ricker子波

假定地震子波已知為q(t),且振幅A=1.0,由(4)式有:

M(g,q):=

=Mqg=Qqg+Pqg=p

(15)

這里定義算子:

(16)

(17)

(18)

同理,要求多次波的估計值pm和一次波估計值pp,與實際地震記錄p也在能量水平上接近,分別用sm和sp表示比例系數(shù),則有:

(19)

(20)

1.3 自適應(yīng)逼近實際地震子波

通常在EPSI算法中,使用SPGL1方法求解,其中子波求解采用常規(guī)的最小平方QR分解算法或匹配濾波算法。這里采用無約束非線性優(yōu)化算法。定義目標(biāo)函數(shù)為:

(21)

則無約束非線性優(yōu)化問題為:

(22)

實現(xiàn)流程如下。

1) 輸入地震記錄p,設(shè)定殘差值σ(‖p‖2的5%～10%)。

2) 根據(jù)地震記錄的主頻fm計算振幅A=1的Ricker子波q(t)。

5) 迭代計數(shù)k=0,L1范數(shù)約束τ0=0。

6) 循環(huán)開始。

7) 由σ和τk,通過SPGL1確定τk+1。

12)k=k+1。

2 應(yīng)用算例

2.1 簡單巖丘模型

設(shè)計一個如圖2所示的頂部有水層的巖丘地質(zhì)模型。該模型的水層深度在200m左右,橫向變化的鹽丘層大致位于深度400～800m。通過時間二階、空間四階有限差分聲波程序進行正演模擬,使用主頻為25Hz的零相位Ricker子波。震源線與接收線精確重合,道間距為15m。

圖2 頂部有水層的簡單巖丘模型

2.1.1. 應(yīng)用效果測試分析

抽取排列在巖丘高點附近的150道數(shù)據(jù)進行處理。圖3a是測線中部相鄰兩炮原始數(shù)據(jù),水層底部和鹽丘的頂部產(chǎn)生的一階及高階多次波(圖中白色箭頭指示處)嚴(yán)重影響了一次波同相軸。圖3b是采用本文方法壓制多次波后的對應(yīng)單炮記錄,圖中白色箭頭指示處的多次波已得到有效壓制,一次波得到很好的恢復(fù),1.6s以下空白,沒有信息,這與模型1700m以下沒有反射層相對應(yīng)。圖3c是采用EPSI方法壓制多次波后的對應(yīng)兩炮記錄,白色箭頭處的多次波也得到壓制。對比圖3b和圖3c可以看出,圖3b中的繞射波更完整,同時最右邊的箭頭和最下邊兩個箭頭處,原來被多次波掩蓋的有效波得到更好的恢復(fù),說明采用本文方法壓制多次波時,有效波得到更好的保護和恢復(fù)。

圖4a為原始數(shù)據(jù)的零炮檢距剖面,多次波嚴(yán)重影響資料品質(zhì)。圖4a中存在兩個產(chǎn)生較強多次波的界面,一個是水層頂面,一個是水層底面,圖中白色箭頭指示了表面相關(guān)一階和高階多次波。水層頂面產(chǎn)生表面相關(guān)多次波,水層底面產(chǎn)生層間多次波。圖4b 為采用本文方法進行多次波壓制后的零炮檢距剖面,可以看出,白色箭頭所指水層表面產(chǎn)生的各階多次波得到大幅壓制,被掩蓋的一次波同相軸得到了很好的恢復(fù),多次波壓制效果良好。在0.8s左右存在由水底產(chǎn)生的層間多次波。在1.2s左右處,存在由水底產(chǎn)生的高階層間多次波。對于這些層間多次波,可采用其它層間多次波壓制方法進一步壓制。圖4c 為采用EPSI進行多次波壓制后的零炮檢距剖面?？梢钥闯?圖4c整體效果與圖4b相當(dāng),但圖4b中的繞射波信息更完整,同時最下面的兩個白色箭頭處的弱反射也完整地顯現(xiàn)出來了。

圖4 采用不同方法進行多次波壓制后的零炮檢距剖面a 原始數(shù)據(jù); b 本文方法;c EPSI

2.1.2 初始子波誤差容忍度分析

為了分析本文方法對初始子波誤差的容忍度,基于同樣的數(shù)據(jù)集,分別采用主頻為20Hz的零相位、最小相位和混合相位Ricker子波作為初始輸入子波,應(yīng)用本文方法進行處理和分析。

不同相位子波輸入的多次波壓制結(jié)果如圖5所示。圖5a為原始單炮記錄(局部),含有嚴(yán)重的表面相關(guān)多次波,圖中白色箭頭指示了主要多次波。圖5b 至圖5d分別為零相位子波、最小相位子波和混合相位Ricker子波輸入時的多次波壓制結(jié)果。從圖5 可以看出,表面相關(guān)多次波得到有效壓制(如圖中白色箭頭所示),一次波得到較好恢復(fù),不同相位子波輸入對多次波壓制效果影響較小,表明本文方法不受輸入子波類型影響,對初始子波誤差的容忍度較高。

不同相位子波輸入時的子波迭代更新結(jié)果如圖6 所示。圖6a為零相位子波迭代更新結(jié)果,初始子波(k=1)經(jīng)過4次迭代,子波的主頻和振幅趨于穩(wěn)定,并且向?qū)嶋H理論子波逼近。求取的子波(k=4)和模擬正演使用的理論子波(主頻為25Hz的零相位Ricker子波)一致。圖6b為最小相位子波迭代更新結(jié)果,初始子波(k=1)經(jīng)過5次迭代,求取的子波(k=5)和模擬正演使用的理論子波一致。圖6c為混合相位子波迭代更新結(jié)果,初始子波(k=1)經(jīng)過5次迭代,求取的子波(k=5)和模擬正演使用的理論子波一致。

采用主頻為20Hz的零相位、最小相位和混合相位Ricker子波作為初始輸入子波,應(yīng)用本文方法進行處理,最終子波都收斂到正演模擬數(shù)據(jù)的實際理論子波,而多次波壓制效果沒有明顯變化,說明子波迭代更新結(jié)果受產(chǎn)生數(shù)據(jù)的實際子波控制,本文方法對子波誤差具有較高的容忍度。

圖5 輸入不同相位子波的多次波壓制結(jié)果a 原始數(shù)據(jù); b 零相位子波; c 最小相位子波; d 混合相位子波

圖6 不同相位子波輸入時的子波迭代更新結(jié)果a 零相位子波; b 最小相位子波; c 混合相位子波

2.2 實際數(shù)據(jù)算例

實際地震數(shù)據(jù)來自蘇伊士灣海上地區(qū)(Gulf of Suez region),淺海底產(chǎn)生了強烈的多次波,每炮178道,采樣間隔為4ms,記錄長度為4s。初始子波使用20Hz零相位Ricker子波。

圖7a為包含多次波的原始單炮記錄,整個數(shù)據(jù)體被強烈的多次波覆蓋。圖7b為采用本文方法進行多次波壓制后的單炮記錄,可以看出,不同階次的多次波都得到很好的壓制,一次波得到恢復(fù)。

圖7 采用本文方法進行多次波壓制前(a)、后(b)的單炮記錄

圖8a為原始數(shù)據(jù)疊加剖面,多次波在剖面中普遍存在。圖8b為采用本文方法進行多次波壓制后的疊加剖面,可以看出,海洋表面相關(guān)各階多次波得到很好的壓制。

采用本文方法迭代求解時的子波變化情況如圖9 所示?？梢钥闯?初始子波(k=1)經(jīng)過5次迭代,子波趨于穩(wěn)定收斂,第4次迭代(k=4)和第5次迭代(k=5)的子波幾乎重合。

圖8 采用本文方法進行多次波壓制前(a)、后(b)的疊加剖面

圖9 采用本文方法迭代求解時子波變化情況

3 結(jié)論

本文以EPSI多次波壓制方法頻域理論公式為基礎(chǔ),將EPSI過程表示為雙線性算子一次波脈沖響應(yīng)和震源子波的約束優(yōu)化問題,提出了一種自適應(yīng)稀疏反演多次波壓制方法(AEPSI),在迭代求解的過程中,震源子波的主頻、振幅、相位等特性,自適應(yīng)逐步逼近實際地震子波,避免了通過輸入記錄的多維互相關(guān)數(shù)據(jù)來估計初始震源子波的不確定性以及求取一次脈沖響應(yīng)時的人工拾取和時窗設(shè)置。模型數(shù)據(jù)和實際資料處理結(jié)果表明,本文方法正確有效,穩(wěn)定收斂,對子波誤差具有較高的容忍度,整個處理過程完全由數(shù)據(jù)驅(qū)動,無需人工干預(yù),具有較好的應(yīng)用前景。