自適應(yīng)稀疏反演多次波壓制方法

楊旭明,王 麗,陶長(zhǎng)江,鮑 偉

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

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

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

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

1 方法原理

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

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

P=GQ+GRP

(1)

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

p=M(g,q;p)

(2)

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

p=M(g,q)

(3)

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

(4)

(5)

(6)

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

(7a)

(7b)

(8)

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

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

(9)

(10)

(11)

(12)

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

(13)

1.2 初始震源子波的確定

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

(14)

式中:t為時(shí)間;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)

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

(19)

(20)

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

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

(21)

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

(22)

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

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

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

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

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

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

12)k=k+1。

2 應(yīng)用算例

2.1 簡(jiǎn)單巖丘模型

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

圖2 頂部有水層的簡(jiǎn)單巖丘模型

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

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

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

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

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

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

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

不同相位子波輸入時(shí)的子波迭代更新結(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)用本文方法進(jìn)行處理,最終子波都收斂到正演模擬數(shù)據(jù)的實(shí)際理論子波,而多次波壓制效果沒有明顯變化,說明子波迭代更新結(jié)果受產(chǎn)生數(shù)據(jù)的實(shí)際子波控制,本文方法對(duì)子波誤差具有較高的容忍度。

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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