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

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

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

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

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

1 方法原理

1.1 自適應稀疏反演多次波壓制原理

為了便于理解和推導,EPSI頻率域方程表示為:

P=GQ+GRP

(1)

式中:P為頻率域數(shù)據(jù)矩陣;G為頻率域一次波脈沖響應;R為頻率域自由表面反射算子;Q為頻率域震源矩陣。如果用q(t)表示對應于時間域數(shù)據(jù)矩陣p的震源特性,并假定q(t)對所有的炮和道是不變的,則Q=w(ω)I。其中,w(ω)是頻率的標量函數(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)的時移。當然,也可以施加其它需要的震源特性。

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

(9)

(10)

(11)

(12)

調整因子sk為:

(13)

1.2 初始震源子波的確定

實際中的地震子波是一個很復雜的問題,因為地震子波與地層巖石性質有關,地層巖石性質本身就是一個復雜體。為了研究方便,對地震子波進行模擬,這就是理論子波。目前普遍認為Ricker提出的地震子波數(shù)學模型具有廣泛的代表性,因此,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 自適應逼近實際地震子波

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

(21)

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

(22)

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

1) 輸入地震記錄p,設定殘差值σ(‖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 應用算例

2.1 簡單巖丘模型

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

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

2.1.1. 應用效果測試分析

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

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

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

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

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

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

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

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

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

圖6 不同相位子波輸入時的子波迭代更新結果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為采用本文方法進行多次波壓制后的單炮記錄,可以看出,不同階次的多次波都得到很好的壓制,一次波得到恢復。

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

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

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

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

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

3 結論

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