深度域地震子波提取方法綜述

張 杰,陳學華,蔣 偉

(1.成都理工大學油氣藏地質(zhì)及開發(fā)工程國家重點實驗室,四川成都610059;2.成都理工大學地球勘探與信息技術(shù)教育部重點實驗室,四川成都610059)

深度偏移在復雜地質(zhì)構(gòu)造和地質(zhì)體的高精度成像方面具有突出的優(yōu)勢[1-3]。隨著疊前深度偏移處理技術(shù)的廣泛應用,業(yè)界對系統(tǒng)的真深度域地震反演和解釋方法技術(shù)提出了日益迫切的需求。深度域地震子波提取是進行深度域地震反演和解釋等工作的基礎。由于深度域地震子波是與地下介質(zhì)速度有關的函數(shù)[4],其波形隨著介質(zhì)速度的增加而拉伸、隨著介質(zhì)速度的降低而壓縮,且在不同介質(zhì)層界面處的波形呈現(xiàn)非對稱特征,因此,“線性時不變”條件在深度域中并不成立,基于褶積模型的地震子波提取方法對深度域地震數(shù)據(jù)并不是直接適用的。如何準確提取深度域地震子波,從而為深度域地震資料反演和解釋等工作奠定可靠的基礎具有重要意義。

目前,常規(guī)深度域地震反演和解釋通常是先將深度域測井數(shù)據(jù)和地震數(shù)據(jù)根據(jù)時深關系轉(zhuǎn)換至時間域,完成反演和解釋后,再將結(jié)果轉(zhuǎn)換回深度域。這樣的流程存在兩個主要缺點:一是在將測井的各物性參數(shù)曲線從深度域轉(zhuǎn)換至時間域的過程中,會損失其中的重要高頻特征信息;二是在不同域之間的轉(zhuǎn)換過程中,會在有效數(shù)據(jù)中引入累積誤差。因此,常規(guī)深度域地震反演和解釋流程,不能充分利用深度域地震數(shù)據(jù)中反映儲層結(jié)構(gòu)、巖性和流體情況等的有效信息,不利于挖掘深度域地震數(shù)據(jù)的應用價值和潛力,也不利于實現(xiàn)儲層的精細描述。為了避免常規(guī)深度域地震反演和解釋中的深-時和時-深轉(zhuǎn)換,需要利用好深度域地震子波的特性。如果介質(zhì)速度是一個常數(shù),那么在該介質(zhì)中的地震子波波形就能在不同深度位置處保持不變,可以滿足“線性時不變”條件。因此,可以采用一定的速度變換方法,將深度域測井數(shù)據(jù)和深度域地震數(shù)據(jù)無損地轉(zhuǎn)換到一個介質(zhì)速度處處相等的等效介質(zhì)(即常速度深度域)中[4-7],從而可以基于褶積模型提取深度域地震子波,進而充分利用深度域測井數(shù)據(jù)和深度域地震數(shù)據(jù)中的有效信息,實現(xiàn)更高質(zhì)量和更加精細的深度域地震反演和解釋[8-12]。

本文首先對深度域地震子波提取中存在的主要問題進行了分析;然后,對目前已有的深度域地震子波提取方法進行了梳理,并對每種方法的優(yōu)缺點及適用條件進行了總結(jié);最后,對深度域地震子波提取方法的未來發(fā)展提出了一些建議。

1 深度域地震子波提取存在的主要問題

1.1 深度域地震子波對介質(zhì)速度的依賴特性

一個波動過程,既可以用時間域函數(shù)表征,也可以用空間域函數(shù)表征。時間域地震子波是在波動延續(xù)時間范圍內(nèi),某一質(zhì)點的各時刻振動幅度形成的振動曲線?？臻g(深度)域地震子波是在某一時刻波動延續(xù)空間范圍內(nèi),一系列質(zhì)點的振動幅度形成的波形曲線。地震子波可以被表示為一系列諧波(一組不同振幅、頻率和相位的正弦曲線)疊加而成的信號[4,13],在時間域和深度域分別表示為:

w(t)=∑w(f)ei2πft

(1)

w(d)=∑w(k)eikd

(2)

式中:w(t)和w(d)分別表示時間域和深度域的地震子波;t和d分別表示時間和深度;w(f)和w(k)分別表示時間域子波的頻譜和深度域子波的波數(shù)譜;k和f分別表示波數(shù)和頻率。波數(shù)與頻率之間的關系為:

(3)

式中:v是地下介質(zhì)的速度,隨深度變化。

為直觀展示時間域地震子波與深度域地震子波的異同,圖1a和圖1b分別顯示了在2層不同介質(zhì)中的時間域地震子波和深度域地震子波以及各自的前20個主要正弦分量。從圖1可以看出:零相位的時間域地震子波波形在不同介質(zhì)的分界面處是對稱的,對應的主要正弦分量在分界面處也是對稱的;深度域地震子波波形在不同介質(zhì)的分界面處是不對稱的,主要的正弦分量在分界面處也是不對稱的,且各正弦分量之間存在不同的相位差。圖2顯示了3層介質(zhì)模型中的深度域合成地震記錄及其在速度為2000,4000,6000m/s的常速度介質(zhì)中的深度域地震子波。圖2通過一個3層水平地質(zhì)模型和與之對應的深度域合成地震記錄,以及每一層常速度介質(zhì)中的深度域地震子波,展示了深度域地震子波波形隨介質(zhì)速度變化的特性。由圖2可以看出,介質(zhì)速度越大,深度域地震子波的波形越寬(對應的波數(shù)越小)。圖1和圖2 所展示的深度域地震子波對介質(zhì)速度的依賴特征在鹽丘體附近的深度域地震數(shù)據(jù)中很常見[14]。因此,與時間域地震子波相比,深度域地震子波除了受到地震震源特征、地下介質(zhì)的吸收和頻散等因素影響外[15-16],還受到介質(zhì)速度的影響。

圖1 在2層不同介質(zhì)中的時間域地震子波和深度域地震子波a 時間域零相位地震子波及構(gòu)成這條振動曲線的前20個時間域正弦分量; b 某一時刻(波形曲線振幅峰值恰好位于介質(zhì)分界面)的深度域地震子波(由101個質(zhì)點組成)及構(gòu)成這條波形曲線的前20個空間域正弦分量

圖2 3層介質(zhì)模型中的深度域合成地震記錄(a)及其在速度為2000,4000,6000m/s的常速度介質(zhì)中的深度域地震子波(都與主頻為30Hz的時間域Ricker子波相對應)(b)

時間域中常用的褶積運算的假設前提是“線性時不變”。對于一個線性時不變系統(tǒng),當輸入一個諧波時,其對應的輸出是一個同頻率但振幅和相位不同的諧波。而(2)式和(3)式表明,在深度域中,“線性時不變條件”不成立。因此,要對深度域地震數(shù)據(jù)進行反演、屬性分析等處理,需要采用一定的變換方法將數(shù)據(jù)轉(zhuǎn)換至一個能夠滿足“線性時不變條件”的轉(zhuǎn)換域中,目前主要有以下2種途徑:①將深度域地震數(shù)據(jù)根據(jù)時-深關系轉(zhuǎn)換到時間域;②將深度域地震數(shù)據(jù)進行速度變換,轉(zhuǎn)換到一種介質(zhì)速度處處相等的均勻介質(zhì)中,即常速度深度域中。這2種途徑的基本思想都是在轉(zhuǎn)換域(時間域或常速度深度域)中完成反演、屬性分析等處理之后,再將處理結(jié)果反變換回深度域。

1.2 深度域測井數(shù)據(jù)與地震數(shù)據(jù)的采樣間隔非一致性

通常,深度域測井數(shù)據(jù)的深度采樣間隔小于深度域地震數(shù)據(jù)的深度采樣間隔。在利用地震數(shù)據(jù)和測井數(shù)據(jù)提取地震子波時,測井數(shù)據(jù)和地震數(shù)據(jù)的采樣間隔須要保持一致。使測井數(shù)據(jù)與地震數(shù)據(jù)的采樣間隔保持一致的方法主要有以下3種:①將深度域測井數(shù)據(jù)按照時-深關系轉(zhuǎn)換到時間域,使測井數(shù)據(jù)與地震數(shù)據(jù)的時間采樣間隔保持一致;②對深度域測井數(shù)據(jù)進行降采樣,使其深度采樣間隔與深度域地震數(shù)據(jù)的采樣間隔保持一致;③對深度域地震數(shù)據(jù)進行插值,使其深度采樣間隔與深度域測井數(shù)據(jù)的采樣間隔保持一致。方法①和方法②都是數(shù)據(jù)的降采樣過程,常規(guī)的降采樣處理方法并未考慮地下反射界面的信息,從而導致有用的層位信息丟失。因此,在對深度域測井數(shù)據(jù)進行降采樣時,應使用能夠盡可能保留反射界面信息的降采樣處理方法[17]。方法③是數(shù)據(jù)的插值過程,該過程不會增加深度域地震數(shù)據(jù)本身的信息,其缺點是會增加數(shù)據(jù)量。但目前的計算機軟硬件條件足以應對這個缺點,故我們推薦使用該方法解決深度域測井數(shù)據(jù)與地震數(shù)據(jù)的采樣間隔非一致性問題。

2 深度域地震子波提取方法

若要基于褶積模型進行深度域地震子波提取,需要先將深度域地震數(shù)據(jù)和深度域測井數(shù)據(jù)轉(zhuǎn)換至時間域或常速度深度域。由于將深度域地震數(shù)據(jù)和測井數(shù)據(jù)轉(zhuǎn)換到時間域(降采樣處理)會丟失數(shù)據(jù)信息,所以,更佳的做法是將數(shù)據(jù)轉(zhuǎn)換至常速度深度域。具體的常速度深度域轉(zhuǎn)換方法流程可以參考文獻[4]和文獻[7]。需要說明的是,深度域轉(zhuǎn)換至常速度深度域是一個數(shù)據(jù)插值過程,不會增加深度域地震數(shù)據(jù)本身的信息。我們對目前已有的深度域地震子波提取方法以及可以借鑒的時間域地震子波提取方法進行了梳理,并根據(jù)在方法中是否用到測井信息,將方法分為確定性地震子波提取方法和統(tǒng)計性地震子波提取方法2大類。

2.1 確定性地震子波提取方法

當線性時不變系統(tǒng)條件得以滿足時,合成地震記錄s可由反射系數(shù)向量r和地震子波向量w進行褶積運算(線性褶積)再加上噪聲項n獲得,其具體形式為:

s=r*w+n

(4)

其中,“*”表示一維褶積運算,噪聲項通常被假設為高斯隨機白噪聲。從(4)式可以看出,當已知地震記錄和反射系數(shù)時,地震子波成了唯一要確定的量,這也是確定性地震子波提取方法名稱的由來。在時間域或常速度深度域中,利用測井數(shù)據(jù)和井旁地震記錄進行地震子波提取的方法就屬于這一類?；隈薹e模型的不同形式,目前的確定性地震子波提取方法主要有譜除法和最小二乘法。除了最小二乘法已經(jīng)用于深度域地震子波提取外[12,18],其余的確定性方法均還未見到在這方面的應用,但這些方法的思想是值得借鑒的。

2.1.1 譜除法

兩個向量的褶積運算在傅里葉變換域(頻域或波數(shù)域)中是兩個向量傅里葉變換結(jié)果的乘積運算:

s′=r′w′+n′

(5)

式中:s′,r′,w′和n′分別表示合成地震記錄、反射系數(shù)、地震子波和噪聲項的傅里葉變換結(jié)果。若忽略噪聲項,地震子波可由地震記錄與反射系數(shù)的傅里葉變換結(jié)果相除獲得[15-16,19]:

(6)

式中:IFT(·)表示反傅里葉變換。從(6)式可以看出,反射系數(shù)的傅里葉變換結(jié)果中若存在零值,將導致結(jié)果異常。可以通過在反射系數(shù)的傅里葉變換結(jié)果中添加一個非常小的常數(shù)以避免這個問題[16]。此外,由于數(shù)量級不同的小數(shù)在進行除法運算后會出現(xiàn)異常大值(主要出現(xiàn)在高頻或高波數(shù)部分),所以還需要對譜除結(jié)果用帶通濾波器作進一步的處理。譜除法對噪聲非常敏感,為了減少噪聲的影響,在實際應用中,更多的是使用反射系數(shù)與地震記錄的互相關結(jié)果和反射系數(shù)的自相關結(jié)果進行譜除[20-21]。由譜除法得到的結(jié)果通常作為迭代優(yōu)化地震子波提取方法的初始子波模型[22]。

2.1.2 最小二乘法

最小二乘法是要使得合成數(shù)據(jù)與已知數(shù)據(jù)之間的殘差能量達到最小,具體的目標函數(shù)可表示為如下形式:

(7)

式中:d為包含M個采樣點的已知地震記錄向量;R為由反射系數(shù)r構(gòu)成的M階Toeplitz方陣[23-24];‖·‖2表示向量的L2范數(shù)。(7)式的解析解形式為:

(8)

一般而言,地震子波的延續(xù)范圍有限[25-26],所以,在提取地震子波的過程中,可以用這一點作為約束,以便在地震數(shù)據(jù)信噪比不高的情況下獲得合理的結(jié)果。已有許多學者提出了不同的地震子波長度確定方法:梁光河[27]將褶積模型看作一個移動平均(moving average,MA)過程,從而將子波長度確定問題轉(zhuǎn)換為MA過程的定階問題;BUNCH等[28]基于合成地震記錄與已知地震記錄之間的均方根誤差,利用試錯法確定子波長度;BULAND等[29]與GUNNING等[30]指出,可以將子波長度作為一個未知參數(shù)包括在子波提取過程中;RIETSCH[31-32]用多道地震記錄的自相關和互相關所構(gòu)建矩陣的零特征值個數(shù)作為所提取地震子波的長度。將地震子波長度作為約束與最小二乘法相結(jié)合后,(8)式改寫為[24]:

(9)

式中:RP=RP,P是根據(jù)地震子波長度L(L

(10)

式中:I為L階單位方陣。由于求解(9)式比較容易且快速,因此,可以通過試錯法確定最合適的地震子波長度L。通過地震子波長度信息的約束,不但可以降低過擬合,還可以降低運算過程中的矩陣維度,提高計算效率。

在確定性地震子波提取方法的實際應用中,由于測井資料的深度范圍有限,與之對應的井旁地震記錄也只有有限的一段數(shù)據(jù),在這段數(shù)據(jù)中包含著在此深度范圍之外的信息,相對而言,這部分信息是要處理的“噪聲”,故還需要考慮在(9)式中加入正則約束項以形成更加穩(wěn)健的算法。一般而言,正則約束項多采用L1范數(shù),但是L1范數(shù)不可求導,對應的目標函數(shù)沒有解析解,只能通過迭代求解,且求解的結(jié)果是稀疏的,因此,最小二乘法中多用L2范數(shù)作為正則約束項[18,24,33],相應的目標函數(shù)及其解析解如下:

(11)

(12)

式中:λ(λ>0)是正則化參數(shù),可以通過試錯法或廣義交叉驗證(generalized cross-validation,GCV)方法確定[34-35]。由于L1范數(shù)對異常值的處理效果優(yōu)于L2范數(shù)[36],所以還可以進一步考慮將L1范數(shù)和L2范數(shù)相結(jié)合的處理方案,例如Hybird、Huber、Biweight范數(shù)等。張杰等[18]采用Huber范數(shù)約束最小二乘法在較短的深度窗中獲得了優(yōu)于L2范數(shù)約束最小二乘法的深度域地震子波提取結(jié)果。此外,BULAND等[29]基于貝葉斯理論提出了一種貝葉斯最小二乘法。貝葉斯理論重視利用先驗信息,貝葉斯最小二乘法是在(12)式的基礎上使用了多層先驗,例如:不僅假設子波滿足高斯分布,而且還在此基礎上,進一步假設該高斯分布的均值和方差分別滿足高斯分布和逆伽馬分布(inverse gamma distribution)。貝葉斯最小二乘法的優(yōu)點是能給出所提取子波結(jié)果的不確定性評價,但是,由于其沒有明確的后驗概率分布表達式,需要通過馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)算法進行大量抽樣,算法耗時長。

2.1.3 其它確定性地震子波提取方法

LIU等[37]使用測井數(shù)據(jù)作為輸入,井旁地震記錄作為輸出,利用誤差反向傳播(back propagation,BP)神經(jīng)網(wǎng)絡方法來提取地震子波。由于訓練好的網(wǎng)絡結(jié)果是一系列的網(wǎng)絡連接權(quán)重系數(shù),并不是一個子波波形,還需要將訓練好的網(wǎng)絡作用于單反射界面模型上,輸出的結(jié)果才是最終提取的地震子波。該方法的缺點是如果井震之間不匹配會導致算法收斂較慢,優(yōu)點是不需要依賴初始子波模型,也不需要對地震子波和反射系數(shù)進行任何假設。

2.2 統(tǒng)計性地震子波提取方法

當工區(qū)無測井資料可用時,就無法獲取反射系數(shù)信息。公式(4)存在兩個未知數(shù),解的多解性是必然的。故想要得到合理的反射系數(shù)和地震子波估計,就需要對二者進行一定的假設[21,38]:①反射系數(shù)是統(tǒng)計獨立的平穩(wěn)隨機過程,其對應的概率密度函數(shù)是對稱的非高斯函數(shù);②地震子波是最小相位或常相位的;③噪聲是統(tǒng)計獨立的平穩(wěn)隨機過程,其均值為0,滿足高斯分布,且與地震記錄和反射系數(shù)是相互獨立的統(tǒng)計量。常用的統(tǒng)計性地震子波提取方法有同態(tài)法、交替迭代法、二階統(tǒng)計量法和高階統(tǒng)計量法。目前,除了二階統(tǒng)計量法已用于提取深變地震子波外[39-40],其余的統(tǒng)計性方法均還未見到在深度域地震子波提取中的應用。此外,在缺少測井速度信息的情況下,如何進行速度變換也制約著統(tǒng)計性方法在深度域地震子波提取中的應用。

2.2.1 同態(tài)法

忽略噪聲項,將(5)式兩邊分別取自然對數(shù),可以得到如下線性關系:

ln(s′)=ln(r′)+ln(w′)

(13)

對(13)式中各變量進行反傅里葉變換,可以得到褶積運算在復倒譜域(complex cepstral domain)中[41]的對應形式:

IFT[ln(s′)]=IFT[ln(r′)]+IFT[ln(w′)]

(14)

(13)式和(14)式構(gòu)成一個同態(tài)系統(tǒng),可以將褶積運算轉(zhuǎn)變?yōu)榧臃ㄟ\算,基于該思想的統(tǒng)計性地震子波提取方法稱為同態(tài)法。與反射系數(shù)相比,地震子波波形較為平滑且其延續(xù)范圍有限,因此,同態(tài)法假設地震子波的復倒譜通常位于原點附近,反射系數(shù)的復倒譜遠離原點,且二者可分離。通過設計相應的低通濾波器,可以從地震數(shù)據(jù)的復倒譜中分離地震子波和反射系數(shù)。同態(tài)法不需要對地震子波的相位和反射系數(shù)的概率分布函數(shù)做出假設,但該方法涉及到相位解纏和移除解纏相位線性趨勢等復雜處理。ULRYCH等[42]使用微分復倒譜來避免相位解纏和線性趨勢移除,并采用多道地震數(shù)據(jù)的復倒譜疊加來突出地震子波的復倒譜。在實際地震數(shù)據(jù)中,地震子波和反射系數(shù)在復倒譜域并不是完全可分的,需要用多道地震記錄的統(tǒng)計平均[43-44]以及隨機時窗來壓制反射系數(shù)的影響[42,45]。DE MACEDO等[46]利用測井反射系數(shù)的復倒譜信息,構(gòu)建了穩(wěn)健的分離濾波器來提取子波,但這實際上屬于一種確定性地震子波提取方法。

2.2.2 交替松弛迭代法

交替松弛迭代法求解兩個未知變量的基本思想是:先給定其中一個變量的初始解,基于此初始解求解另一個變量;然后固定求解后的變量,再求解另一個變量,如此反復交替求解兩個變量,直至滿足停機準則。對于地震子波提取問題而言,由于有很多成熟的方法(如高階統(tǒng)計量法)可以僅利用地震數(shù)據(jù)信息獲得一個較為可靠的地震子波估計[47],所以,通常先給定子波的初始解。交替松弛迭代法的目標函數(shù)為:

(15)

式中:w0是給定的子波初始解;Reg(r)是反射系數(shù)約束項;λ1和λ2分別是地震子波約束項和反射系數(shù)約束項的正則參數(shù)。當固定反射系數(shù)求解地震子波時,(15)式退化為(11)式,可以采用(12)式的解析解形式直接求解子波。交替松弛迭代法的重點是求解反射系數(shù)。曹靜杰[48]在反射系數(shù)滿足廣義高斯分布的假設前提下,采用Lp(0

2.2.3 二階統(tǒng)計量法

二階統(tǒng)計量,即信號的自相關,其結(jié)果并不包含相位信息,但由于計算簡單,通常主要用于估計地震子波的振幅譜。一個零均值、統(tǒng)計獨立的平穩(wěn)隨機過程y(t),其二階統(tǒng)計量定義為[50]:

(16)

式中:N是信號的采樣點個數(shù);τ是時間延遲量?；诜瓷湎禂?shù)的非高斯白噪假設,其自相關是一脈沖函數(shù),可以近似當作一個常數(shù)處理,(4)式所示的褶積模型對應的二階統(tǒng)計量形式為:

c2s(τ)=σc2w(τ)

(17)

式中:σ是反射系數(shù)自相關的近似;c2s(τ)和c2w(τ)分別為地震記錄和子波的自相關。由于自相關的傅里葉變換就是功率譜,因此,可以利用地震記錄的功率譜估計地震子波的功率譜。ZHANG等[51]利用基于局部屬性的時頻分解方法對時間域地震數(shù)據(jù)進行時頻分解,然后基于瞬時振幅譜提取了零相位時變地震子波;ZHANG等[39]利用S變換對疊后深度域地震數(shù)據(jù)進行深度-波數(shù)分解,然后基于瞬時波數(shù)振幅譜提取了零相位深變地震子波;之后,ZHANG等[40]又將該方法擴展到疊前深度域地震數(shù)據(jù),提取了隨深度和入射角度變化的地震子波集。然而,在對深度域地震數(shù)據(jù)使用深度-波數(shù)分解方法時,需要考慮深度域地震子波隨地下介質(zhì)速度變化的特征,否則在利用深度域地震數(shù)據(jù)的深度-波數(shù)分解結(jié)果進行衰減頻散屬性分析時,就有可能出現(xiàn)誤判。此外,ROSA等[52]提出一種譜模擬方法:

|w(f)|=|f|neH(f)

(18)

式中:|w(f)|是時間域地震子波的振幅譜;n是大于等于0的實數(shù);H(f)是和頻率有關的多項式。利用(18)式,通過擬合地震記錄的功率譜,可以得到地震子波的振幅譜。與時頻分析方法(如S變換)相結(jié)合,譜模擬方法可用于提取時變地震子波的振幅譜。譜模擬方法的前提假設是,地震子波的振幅譜是一光滑的單峰曲線,若地震數(shù)據(jù)的振幅譜是多峰曲線(如短時窗中地震數(shù)據(jù)的振幅譜),不建議使用該方法。

2.2.4 高階統(tǒng)計量法

通常,將二階以上的統(tǒng)計量稱為高階統(tǒng)計量。常用的高階統(tǒng)計量主要是三階統(tǒng)計量和四階統(tǒng)計量。高階統(tǒng)計量包含高階累積量和高階矩兩個概念,高階累積量是一階至高階矩的和[50]。由于三階和四階累積量的維度分別是二維和三維,在計算過程中意味著將單道一維地震數(shù)據(jù)擴展到二維和三維,導致算法復雜、計算量大,限制了其在實際中的應用。但由于高階統(tǒng)計量包含了數(shù)據(jù)的相位信息,因此多將其用于估計地震子波的相位。高階統(tǒng)計量法實際上又分為兩類:①利用地震數(shù)據(jù)的高階統(tǒng)計量提取子波,主要是高階累積量擬合法;②利用地震數(shù)據(jù)的高階統(tǒng)計量提取子波的相位譜(子波的振幅譜通過二階統(tǒng)計量法獲得),常用方法主要是Kurtosis最大化準則法和高階累計量譜估算法。

一個零均值、統(tǒng)計獨立的平穩(wěn)隨機過程y(t),其三階和四階累積量定義為[50]:

(19)

y(t+τ3)-m4g(τ1,τ2,τ3)

(20)

式中:m4g(τ1,τ2,τ3)是一個與y(t)有著相同均值和自相關的高斯過程的四階矩,其定義為:

m4g(τ1,τ2,τ3)=c2y(τ1)c2y(τ3-τ2)+c2y(τ2)·

c2y(τ3-τ1)+c2y(τ3)c2y(τ2-τ1)

(21)

可以看出,一個隨機過程的概率密度函數(shù)若是高斯函數(shù),那么其三階和四階累積量均為0。(4)式中的噪聲項通常被假設為滿足高斯分布,故利用高階累積量可以壓制隨機白噪聲的影響。

2.2.4.1 高階累積量擬合法

(4)式褶積模型對應的四階累積量形式為[38]:

c4s(τ1,τ2,τ3)=c4r(τ1,τ2,τ3)*m4w(τ1,τ2,τ3)

(22)

式中:c4s(τ1,τ2,τ3)和c4r(τ1,τ2,τ3)分別為地震記錄和反射系數(shù)的四階累積量;m4w(τ1,τ2,τ3)為地震子波的四階矩。(22)式中的“*”對應三維褶積運算。(4)式褶積模型對應的三階累積量形式與(22)式類似。由于統(tǒng)計性地震子波提取方法對反射系數(shù)的非高斯白噪假設,其三階累積量為0,而其四階累積量不為0,因此更常采用四階累積量[38,53]。此外,反射系數(shù)序列長度越長,其四階累積量就越接近一個脈沖函數(shù)(近似于一個非零常數(shù)),因此(22)式進一步改寫為:

c4s(τ1,τ2,τ3)=σm4w(τ1,τ2,τ3)

(23)

式中:σ是反射系數(shù)四階累積量的近似?？梢钥闯?地震子波的四階矩與地震記錄的四階累積量之間只相差一個比例系數(shù),因此,可利用多道地震記錄的四階累積量疊加提取地震子波,相應的目標函數(shù)為:

(24)

LAZEAR[38]利用最速梯度下降法,通過擬合地震記錄的四階累積量提取地震子波。該方法的一個缺點是容易陷入初始模型附近的局部最優(yōu)值。此外,該方法還受限于地震數(shù)據(jù)的信噪比和數(shù)據(jù)量:所用地震數(shù)據(jù)的信噪比越高、數(shù)據(jù)量越大,反射系數(shù)的四階累積量才會越接近一個非零常數(shù),該方法才會越有效。HARGREAVES[54]進一步指出,若是地震數(shù)據(jù)的帶寬太窄,會導致數(shù)據(jù)的高階累積量中沒有足夠的相位信息,進而導致該方法失效,但目前的寬頻采集技術(shù)已能夠避免這個問題。為避免梯度法陷入局部最優(yōu),VELIS等[55]使用了混合方法:先使用模擬退火優(yōu)化算法找到初始子波,再用梯度下降法求得最優(yōu)子波。戴永壽等[56]基于(23)式,用自回歸滑動平均(autoregressive moving average,ARMA)方法對地震子波進行建模,只需要求解少量模型參數(shù)即可提取子波?？紤]到實際地震數(shù)據(jù)中噪聲的概率密度函數(shù)不一定是高斯函數(shù),故其四階累積量不一定為0,還需要對地震數(shù)據(jù)的高階累積量進行加窗平滑,以壓制噪聲的影響[55]。

2.2.4.2 峰度(Kurtosis)最大化準則法

Kurtosis反映了隨機變量概率密度分布曲線的峰部尖度:Kurtosis值越大,峰的形狀就越尖銳。Kurtosis是數(shù)據(jù)四階累積量在原點處的值,歸一化的Kurtosis定義為[57]:

(25)

對于(25)式,若忽略等號右邊的常數(shù)項,就是最大方差模[58]。反射系數(shù)的四階累積量可以近似為其峰度,且峰度值越大,反射系數(shù)越稀疏[59]。

陸文凱等[60]先將多道地震數(shù)據(jù)的平均對數(shù)振幅譜作為估計子波的振幅譜,然后給定常相位掃描范圍和步長,通過常相位掃描,分別找出使得每一道地震記錄的方差模達到最大的常相位,再將多道地震記錄的常相位求平均后作為子波的常相位,最后基于振幅譜和常相位構(gòu)建最終的子波提取結(jié)果。李鯤鵬等[61]先采用譜模擬方法[52]獲得子波的振幅譜,然后找出使得反褶積結(jié)果的Kurtosis值達到最大時的子波Z變換在單位圓內(nèi)的零點個數(shù),最后基于單位圓內(nèi)零點個數(shù)和振幅譜構(gòu)建最終的子波提取結(jié)果。VAN DER BAAN[62]基于分段平穩(wěn)假設,先利用相互重疊的時窗將地震數(shù)據(jù)分段,并將每段時窗中地震數(shù)據(jù)振幅譜的疊加平均結(jié)果作為子波的振幅譜,然后用該振幅譜的反傅里葉變換結(jié)果作為初始子波,基于Kurtosis最大準則法和常相位掃描,確定子波的常相位,進而獲得該時窗對應的地震子波。對每個時窗進行相同的操作后,就得到了時變地震子波。但由于對地震記錄的分段數(shù)量有限,故提取的時變地震子波相位是分段(非連續(xù))變化的。VAN DER BAAN等[63]進一步使用負熵(Negentropy,即廣義的Kurtosis)以提高子波相位的估計精度。VAN DER BAAN等[64]將單個時窗內(nèi)的Kurtosis最大化問題轉(zhuǎn)化為整個地震剖面上的正則化最小二乘優(yōu)化問題,以消除分段平穩(wěn)假設,使所提取的時變地震子波相位是連續(xù)變化的。Kurtosis最大化準則法的優(yōu)點是計算方便且對高斯噪聲不敏感,其缺點是不適用于當反射系數(shù)的概率密度函數(shù)是弱非高斯函數(shù)以及地震數(shù)據(jù)量少的情況。

2.2.4.3 高階累積量譜估算法

對(19)式的三階累積量進行傅里葉變換,可得到:

=|By(ω1,ω2)|exp[iφ(ω1,ω2)]

(26)

其中,By(ω1,ω2)稱為雙譜,其振幅譜|By(ω1,ω2)|和相位譜φ(ω1,ω2)分別為:

|By(ω1,ω2)|=|y(ω1)||y(ω2)||y(ω1+ω2)|

(27a)

φ(ω1,ω2)=φ(ω1)+φ(ω2)+φ(ω1+ω2)

(27b)

式中:|y(ω)|和φ(ω)分別為信號y(t)的振幅譜和相位譜。

同樣,對(20)式的四階累積量進行傅里葉變換,可得:

e-i(ω1τ1+ω2τ2+ω3τ3)=|By(ω1,ω2,ω3)|·

exp[iφ(ω1,ω2,ω3)]

(28)

其中,By(ω1,ω2,ω3)稱為三譜,其振幅譜|By(ω1,ω2,ω3)|和相位譜φ(ω1,ω2,ω3)分別為:

|By(ω1,ω2,ω3)|=|y(ω1)|·

|y(ω2)||y(ω3)||y(ω1+ω2+ω3)|

(29a)

φ(ω1,ω2,ω3)=φ(ω1)+φ(ω2)+φ(ω3)+

φ(ω1+ω2+ω3)

(29b)

由(27b)式和(29b)式可以看出,雙譜和三譜的相位譜具有線性關系,基于該線性關系估算子波的相位相對容易。

褶積模型的雙譜、三譜頻域表達式形式與(5)式類似:

Bs(ω1,ω2)=Br(ω1,ω2)Bw(ω1,ω2)

(30)

Bs(ω1,ω2,ω3)=Br(ω1,ω2,ω3)Bw(ω1,ω2,ω3)

(31)

由于反射系數(shù)的高階累積量可以近似為高維度的脈沖函數(shù)[38],對應的高階累積量譜可以近似看作一個常數(shù),因此(30)式和(31)式可改寫為:

Bs(ω1,ω2)=σBw(ω1,ω2)

(32)

Bs(ω1,ω2,ω3)=σBw(ω1,ω2,ω3)

(33)

可以看出,地震數(shù)據(jù)與地震子波之間的雙譜和三譜均是成比例的,因此可利用地震數(shù)據(jù)的高階累積量譜信息提取地震子波。此外,也可以基于高階累積量相位譜的線性關系,通過構(gòu)建對應的矩陣形式,用最小二乘法求解子波的相位譜[65-69]。但由于存在相位纏繞問題,因此需要進行相位解纏,否則會影響子波的相位提取結(jié)果。YU等[69]采用保角變換解決相位纏繞問題;戴永壽等[70]和DAI等[71]將基于三階累積量相位譜的子波相位提取方法與基于時頻分解的譜模擬方法相結(jié)合,提取了時變地震子波。由于地震數(shù)據(jù)中噪聲不一定滿足高斯分布,也不一定與地震數(shù)據(jù)相互統(tǒng)計獨立,故在使用高階累積量譜估算法進行地震子波提取時,需要使用一些窗函數(shù)(如Parzen窗)壓制噪聲的影響[72]。

2.2.5 其它統(tǒng)計性地震子波提取方法

RIETSCH[31-32]基于求解2項或2項以上多項式公因式的思想,提出了一種從多道地震記錄中分離反射系數(shù)和地震子波的方法。該方法假設這些地震記錄含有相同的子波,但每道的反射系數(shù)不同。首先,利用多道地震記錄的Z變換系數(shù)構(gòu)成一個矩陣,計算該矩陣的特征值,并將特征值中的零特征值個數(shù)作為地震子波的長度;然后,用該長度作為約束,計算矩陣最小特征值對應的特征向量,該特征向量即為多道地震記錄所對應的反射系數(shù)拼接而成的向量;最后,用提取的反射系數(shù)和地震記錄采用最小二乘法獲得最終的地震子波。該方法不需要對子波的相位和反射系數(shù)的概率分布做出假設,同時對數(shù)據(jù)量要求較小。但由于地震子波通常是空變的,不符合該方法的假設,所以該方法不適用于從構(gòu)造復雜區(qū)域的地震數(shù)據(jù)中提取地震子波。王德營等[73]基于單輸入多輸出系統(tǒng)(single input multiple output,SIMO)的假設,即假設相鄰2道或多道地震記錄的反射系數(shù)相同,而地震子波不同,通過計算多道地震記錄的二階統(tǒng)計量(即自相關),先將反子波計算出來,再根據(jù)反子波求取子波并進行反褶積。該方法的優(yōu)點也是不需要對子波的相位和反射系數(shù)的概率分布做出假設,對數(shù)據(jù)量要求也不高;缺點是該方法并不能直接得到子波結(jié)果,而是得到反子波結(jié)果,需要進一步的計算才能得到地震子波,且其中涉及到多個大型矩陣的求逆運算。孔德輝等[74]將時變子波提取問題轉(zhuǎn)化為在線字典學習問題,認為字典中的每個原子(一個固定的子波)代表局部時窗內(nèi)時變子波的一個分量,通過原子的線性組合可以實現(xiàn)對時變子波的有效逼近。該方法依賴于對時窗的選取,對于復雜地震數(shù)據(jù),如何選取合適的窗函數(shù)以及窗函數(shù)如何隨深度變化還需進一步探索。

2.3 對深度域地震子波提取的幾點認識

1) 目前的深度域地震子波提取方法都需要通過速度轉(zhuǎn)換處理,將深度域的地震數(shù)據(jù)和測井數(shù)據(jù)轉(zhuǎn)換至常速度深度域以完成子波提取。因此,需要進一步發(fā)展不用速度轉(zhuǎn)換處理的深度域地震子波提取方法,例如:基于波動方程理論的子波提取方法,或是基于機器學習的子波提取方法等。此外,每一種地震子波提取方法都有各自的優(yōu)缺點和適用條件,并沒有絕對的優(yōu)劣之分。在具體應用中,需要結(jié)合工區(qū)的實際情況,以確定采用何種地震子波提取方法。

2) 大多數(shù)地震子波提取方法基于單道地震數(shù)據(jù)計算,不利于體現(xiàn)地震子波隨空間變化這一基本特征。LECOMTE[75]指出,地震子波提取方法需要考慮地震數(shù)據(jù)采集過程中“照明條件”的影響。雖然基于一維褶積模型的地震子波提取方法運算效率高,但是無法體現(xiàn)“照明條件”對深度域地震數(shù)據(jù)的影響,從而導致后續(xù)的地震反演結(jié)果出現(xiàn)錯誤。因此,一維深度域地震子波提取方法還需要擴展到二維和三維的深度域進行地震子波提取。

3) 如何恢復地震子波相位是地震子波提取中的關鍵。振幅譜相同的地震子波,其波形會由于相位譜的不同而不同。YUAN等[76]通過正演模擬指出,振幅譜相同而波形不同的地震子波,若以合成地震記錄與井旁地震記錄的相關性為判斷標準,那么這些子波都是“合格”的;但實際上只有相位正確的子波才能得到正確的波阻抗反演結(jié)果。換言之,對子波提取結(jié)果優(yōu)劣的評判還應該與測井的阻抗信息關聯(lián)起來。使用Kurtosis最大化準則法也會出現(xiàn)相似的問題:由于Kurtosis值越大,對應的反射系數(shù)反演結(jié)果越稀疏,直觀上,就是突出強振幅的反射系數(shù)而壓制弱振幅的反射系數(shù)。若將Kurtosis最大化準則用于迭代方法中,會導致最終的反褶積結(jié)果不一定是最優(yōu)的[45]。此外,也應注意到測井數(shù)據(jù)與地震數(shù)據(jù)之間的匹配不一定是最優(yōu)的,在子波提取過程中需要不斷地根據(jù)當前子波提取結(jié)果調(diào)整井-震之間的匹配。ZHANG等[77]考慮到測井數(shù)據(jù)與地震數(shù)據(jù)之間的不匹配,在子波提取過程中用動態(tài)時間規(guī)整(dynamic time warping,DTW)方法更新與當前提取子波對應的井-震時深關系,直至基于當前提取子波的阻抗反演結(jié)果與測井阻抗達到最佳匹配。故筆者建議使用交互式地震子波提取方法,并將阻抗反演結(jié)果與測井阻抗的相關性納入到子波提取結(jié)果可靠性的評判標準中。

4) 若使用交互式地震子波提取方法,那么對方法的計算效率方面就有較高要求。對地震子波用較少的參數(shù)進行建模是提高方法計算效率的有效途徑之一。在基于地震數(shù)據(jù)擬合的子波提取方法中(如最小二乘法、高階累積量擬合法、高階累積量譜估算法等),子波的每一個采樣點都是要求解的未知量,如果能用含有少量參數(shù)的子波模型來替代,那么就可以減少未知量的計算量,從而提高計算效率[56,78-80]。關于地震子波模型,可以進一步參考ALDRIDGE[81]、RYAN[82]、俞壽朋[83]、張海燕等[84]、WANG[85]、SKAUVOLD等[22]的工作。

3 結(jié)束語

隨著計算機計算性能的不斷提升以及疊前深度偏移處理方法的不斷改善,今后深度域地震數(shù)據(jù)將會越來越豐富。然而,目前的深度域地震數(shù)據(jù)反演和解釋幾乎都是先將深度域地震數(shù)據(jù)轉(zhuǎn)換至時間域,在時間域完成反演和解釋后,再將結(jié)果轉(zhuǎn)換回深度域。在時-深和深-時轉(zhuǎn)換過程中,不僅會丟失能夠反映儲層結(jié)構(gòu)、巖性和流體情況等的有效信息,還會產(chǎn)生轉(zhuǎn)換累計誤差,這無益于降低油氣勘探和開發(fā)中的鉆探風險和成本,因此,對深度域地震數(shù)據(jù)的反演和解釋今后必將是直接在真深度域進行,發(fā)展真深度域的地震反演和屬性分析等方法也勢在必行。

在不考慮震源特征以及地下介質(zhì)的吸收頻散等因素的影響下,深度域地震子波由于受到地下介質(zhì)速度的影響,其波形在不斷發(fā)生變化,且在不同介質(zhì)分界面處的波形呈現(xiàn)非對稱特征。因此,褶積模型以及基于褶積模型的方法對深度域地震數(shù)據(jù)并不是直接適用的。如何準確提取深度域地震子波,對真深度域的地震反演和解釋等工作具有重要意義。筆者總結(jié)了目前深度域地震子波提取中存在的問題,整理和概括了現(xiàn)有的深度域地震子波提取方法。雖然目前的深度域地震子波提取方法還比較少,但是有很多時間域地震子波提取方法值得借鑒,筆者對這些時間域地震子波提取方法也一并進行了整理和概括。最后,對深度域地震子波提取方法的未來發(fā)展提出以下建議。

1) 統(tǒng)計性深度域地震子波提取方法有待進一步發(fā)展。如何利用已知速度信息(如具有地質(zhì)意義的速度模型),并在考慮了深度域地震子波隨介質(zhì)速度變化特征的基礎上,直接利用深度域地震數(shù)據(jù)的統(tǒng)計信息提取深度域地震子波是值得探索的。

2) 對于確定性地震子波提取方法,合成地震記錄與井旁地震記錄之間相關性高低仍然是評價確定性地震子波提取結(jié)果優(yōu)劣的主要標準。但是,最高的相關性并不意味著最優(yōu)的地震子波提取結(jié)果,還需進一步考慮地震反演結(jié)果與測井資料之間的相關性。

3) 使用參數(shù)化的地震子波模型,以減少地震子波提取中需要求解的未知量個數(shù),提升方法的計算效率。

4) 考慮地震數(shù)據(jù)采集過程中“照明條件”因素對深度域地震數(shù)據(jù)的影響,將一維深度域地震子波提取方法擴展至提取二維和三維深度域地震子波。

5) 利用機器學習的優(yōu)勢,建立不依賴褶積模型的直接深度域地震子波提取方法,或是建立可以省去地震子波提取步驟,直接完成深度域地震反演的方法。

致謝:感謝中海油田服務股份有限公司對本文研究工作的資助與支持,同時也感謝中海油田服務股份有限公司物探事業(yè)部特普公司的但志偉、肖為等專家對本文研究工作提出的許多有益意見和建議。