基于時(shí)變子波的基追蹤地震反演

周玉毅, 周懷來,2, 廖璐瑤, 王元君,2, 張 舜

(1.成都理工大學(xué) 地球物理學(xué)院，成都 610059；2.“油氣藏地質(zhì)及開發(fā)工程”國家重點(diǎn)實(shí)驗(yàn)室，成都 610059)

0 引言

地震子波提取的研究由來已久，許多研究者相繼提出和發(fā)展了多種提取地震子波方法[1-3]。地震子波估計(jì)準(zhǔn)確程度直接影響著反演的精度與可靠性，傳統(tǒng)的地震子波估計(jì)方法應(yīng)用前提是,地震波在傳播的過程中其波形以及能量保持不變，然而由于地下介質(zhì)對(duì)地震波能量的衰減吸收，導(dǎo)致地震子波在地層傳播過程中將產(chǎn)生高頻成分缺失、相位畸變以及頻散現(xiàn)象，即實(shí)際地震記錄中的地震子波具有時(shí)變性，因此基于Robinson平穩(wěn)褶積模型估計(jì)子波的理論是不完善的。

針對(duì)地震子波的非平穩(wěn)特性，Clarke[4]首次提出非平穩(wěn)褶積模型;Margrave等[5]在非平穩(wěn)褶積模型基礎(chǔ)上研究了地震子波具有時(shí)變性的原因;彭才等[6]建立符合地震子波衰減特性的動(dòng)態(tài)褶積模型，在隨機(jī)白噪反射系數(shù)和地震子波最小相位的假設(shè)前提下，通過對(duì)非平穩(wěn)地震道進(jìn)行加高斯窗傅里葉變換后，可以獲得地震道振幅譜，再經(jīng)方形窗平滑處理得到動(dòng)態(tài)子波；Margrave等[7]在Gabor變換的基礎(chǔ)上，模擬各個(gè)時(shí)窗內(nèi)地震記錄振幅譜提取動(dòng)態(tài)子波。上述方法應(yīng)用的前提是反射系數(shù)序列為白噪聲的假設(shè)，后續(xù)研究者采用譜模擬技術(shù)對(duì)其進(jìn)行規(guī)避。Rosa等[8]結(jié)合Ricker的研究結(jié)果，認(rèn)為子波在頻率域內(nèi)的曲線趨勢是光滑連續(xù)的，與之相對(duì)應(yīng)的是反射系數(shù)在頻率域內(nèi)是振蕩的，地震記錄的振幅譜主要受地震子波影響，而反射系數(shù)只是對(duì)細(xì)節(jié)部分有所影響，因此可以通過模擬地震記錄振幅譜來提取子波。由于譜模擬方法對(duì)于反射系數(shù)序列有很好的適應(yīng)性，因此，諸多學(xué)者將該方法拓展到高精度時(shí)頻分析方法中來提高時(shí)頻聚焦性。李振春等[9]利用S變換將地震記錄轉(zhuǎn)換到頻率域，在頻率域內(nèi)用數(shù)學(xué)方法的多項(xiàng)式擬合出時(shí)變子波的振幅譜，在地震子波相位為零的假設(shè)下，得到頻率域中的時(shí)變子波;張漫漫等[10]利用廣義S變換結(jié)合時(shí)頻濾波求取時(shí)頻譜，對(duì)時(shí)頻譜進(jìn)行多項(xiàng)式譜模擬來提取時(shí)變子波;Zhou等[11]對(duì)譜模擬反褶積方法進(jìn)行推廣，提出基于廣義S變換的動(dòng)態(tài)反褶積算法，估算出的子波準(zhǔn)確性有所提高。在考慮實(shí)際地層吸收衰減基礎(chǔ)上，時(shí)變子波譜可近似表示為震源子波譜與衰減函數(shù)的乘積形式，因此可以從實(shí)際地震記錄中提取時(shí)變子波進(jìn)一步提高反射系數(shù)反演精度。姚振岸等[12]對(duì)非平穩(wěn)地震道進(jìn)行基追蹤譜分解獲得各時(shí)刻的局部功率譜，采用廣義地震子波構(gòu)建方法提取時(shí)變子波，最終實(shí)現(xiàn)非平穩(wěn)地震道的基追蹤反演;李婧等[13]基于改進(jìn)廣義S變換提取符合地震數(shù)據(jù)非平穩(wěn)特性的時(shí)變子波，建立時(shí)變子波矩陣，通過求解對(duì)應(yīng)的稀疏約束正則化目標(biāo)函數(shù)獲得與實(shí)際接近的反射系數(shù)。

筆者以前人的研究為基礎(chǔ)，利用廣義S變換譜模擬方法從地震資料中提取時(shí)變子波，用估算的時(shí)變子波逐點(diǎn)替換時(shí)不變子波構(gòu)建時(shí)變子波庫，對(duì)地震記錄采用時(shí)變子波進(jìn)行基追蹤反演。理論模型和實(shí)際工區(qū)地震資料測試結(jié)果表明，基于本文方法提取的時(shí)變子波基追蹤地震反演的結(jié)果,與傳統(tǒng)時(shí)不變子波基追蹤反演的結(jié)果相比，井震標(biāo)定吻合度更高、橫向連續(xù)性更好以及細(xì)節(jié)展示更為豐富，說明對(duì)非平穩(wěn)地震信號(hào)進(jìn)行基追蹤反演時(shí)采用時(shí)變子波庫的有效性。

1 方法原理

1.1 基于廣義S變換的時(shí)變子波提取

Rosa結(jié)合Ricker的研究結(jié)果分析認(rèn)為：子波振幅譜的曲線趨勢是光滑連續(xù)的，與之相對(duì)應(yīng)，反射系數(shù)振幅譜的曲線趨勢是波動(dòng)的。因此，只要設(shè)置適當(dāng)?shù)臄M合參數(shù)，就能采用數(shù)學(xué)手段對(duì)地震記錄的振幅譜利用多項(xiàng)式擬合進(jìn)行處理，最后得到子波振幅譜[8]。現(xiàn)在假設(shè)一個(gè)地震子波的振幅譜W(t,f)曲線趨勢與雷克子波單峰光滑連續(xù)的振幅譜相似，則可以將其擴(kuò)展為基于時(shí)頻變換的時(shí)變子波，從而在頻率范圍內(nèi)構(gòu)成一種可用的子波數(shù)學(xué)分析模式：

(1)

其中:f為頻率；k、n為常數(shù)；an為關(guān)于f的多項(xiàng)式的系數(shù)。通常情況下，0≤k≤3，4≤n≤7[9]。

描述信號(hào)時(shí)間和頻率聯(lián)合分布的函數(shù)有很多種，如STFT、CWT、Hilbert-Huang Transform等，為了實(shí)現(xiàn)局部時(shí)頻分析，Stockwell等[14]繼承和擴(kuò)展了STFT和CWT的優(yōu)勢方面，提出了一種與其他時(shí)頻分析方法相同的無損可逆S變換，該方法對(duì)信號(hào)進(jìn)行處理后，可以觀察到信號(hào)在頻域和時(shí)域的能量分布。設(shè)x(t)為平方可積的函數(shù)，則x(t)的S變換可以定義為：

(2)

式中：τ表示時(shí)間;f表示頻率;ST(x(τ,f))表示x(t)的S變換。

其中S變換的基本小波函數(shù)定義為：

gf(t)exp(-i2πft)

(3)

式中g(shù)f(t)為高斯窗函數(shù)，可定義為式(4)。

(4)

由于基本小波窗函數(shù)形態(tài)固定，S變換不能夠依據(jù)實(shí)際情況來調(diào)整窗口的大小，在高頻部分分辨率低，導(dǎo)致在實(shí)際應(yīng)用中不能得到較好的效果。針對(duì)S變換的問題，陳學(xué)華等[15]在Stockwell的S變換基礎(chǔ)上提出了廣義S變換(Generalized S Transform)，通過加入兩個(gè)參數(shù)rgs和ρ來調(diào)節(jié)高斯窗函數(shù)的大小，其數(shù)學(xué)表達(dá)式為式(5)。

exp(-i2πft)dt

(5)

通過調(diào)整控制rgs和ρ的大小，該方法能夠?qū)λ械臅r(shí)頻譜具有較高的時(shí)頻聚焦性，可以用來提取不同時(shí)間點(diǎn)的時(shí)變子波，利用廣義S變換時(shí)頻譜提取子波過程(圖1)如下：

圖1 時(shí)變子波提取流程圖

1)對(duì)沿t方向的非平穩(wěn)地震記錄進(jìn)行廣義S變換，得到時(shí)頻譜|GST(x(τ,f))|。

2)固定某一個(gè)時(shí)間T，通過設(shè)置合理的參數(shù)k和n，采用數(shù)學(xué)手段對(duì)該時(shí)間點(diǎn)地震記錄的振幅譜|GST(x(T,f))|利用多項(xiàng)式擬合進(jìn)行處理，擬合得到相應(yīng)的子波振幅譜|GST(W(T,f))|。

3)改變?chǔ)拥娜≈?，逐點(diǎn)計(jì)算子波振幅譜，從而實(shí)現(xiàn)“時(shí)變子波”的時(shí)頻譜提取。

4)在子波相位為零的假設(shè)前提下[9]，對(duì)各時(shí)刻的子波振幅譜進(jìn)行傅里葉逆變換后，可求得各時(shí)刻零相位時(shí)變子波時(shí)域表達(dá)式。

1.2 時(shí)變子波基追蹤反演

地震記錄通常是由一個(gè)地震子波和地層內(nèi)部反射界面的反射系數(shù)進(jìn)行褶積而得到的，用公式表示如下：

s(t)=w(t)*r(t)

(6)

其中：w(t)為一個(gè)地震子波；r(t)為一個(gè)反射系數(shù)序列；“*”為褶積運(yùn)算的符號(hào)。

將式(6)寫為矩陣形式:

s=Wr

(7)

其中：假設(shè)地震數(shù)據(jù)的采樣點(diǎn)為M；s是由地震數(shù)據(jù)所構(gòu)成的M×1矩陣;W是由地震子波w(t)經(jīng)過逐點(diǎn)移動(dòng)構(gòu)成的M×M對(duì)角矩陣；r是由地下反射系數(shù)構(gòu)成的M×1矩陣。雖然式(6)和式(7)一個(gè)是褶積運(yùn)算，一個(gè)是矩陣運(yùn)算，但是它們的運(yùn)算結(jié)果是完全相等的。

基追蹤是一種基于范數(shù)約束的最優(yōu)化準(zhǔn)則，基追蹤反演結(jié)果與實(shí)際地層相接近，分辨率較高，近年來逐漸被應(yīng)用于地震反演，其中目標(biāo)函數(shù)被定義為式(8)[16]。

J=argmin||s-Wr||2+λ||r||1

(8)

其中：J為目標(biāo)函數(shù)；下標(biāo)1和2為目標(biāo)向量，分別用于表示L1和L2范數(shù)；λ是影響調(diào)整值的因子。滿足J→min即為一個(gè)解函數(shù)向量r的最終函數(shù)值。

為了表征地下介質(zhì)的稀疏性，而不是反射系數(shù)的稀疏，根據(jù)Helholm定理，將原始反射系數(shù)對(duì)分解成如圖2所示的偶脈沖對(duì)與奇脈沖對(duì)的線性組合，則任意反射系數(shù)序列可表示為式(9)。

bm,nro(t,m,n,Δt)}

(9)

其中：ro(t,m,n, Δt)為奇楔反射系數(shù)對(duì)分量；re(t,m,n, Δt)為偶楔反射系數(shù)對(duì)分量；am,n和bm,n為對(duì)應(yīng)奇、偶反射系數(shù)對(duì)對(duì)應(yīng)的系數(shù)；Δt為采樣間隔；m、n分別為地下介質(zhì)層頂、底界面所對(duì)應(yīng)的采樣點(diǎn)位置；由于不同反射對(duì)出現(xiàn)的位置不固定，頂?shù)追瓷湎禂?shù)以平移量mΔt沿著圖2中t的方向進(jìn)行逐點(diǎn)移動(dòng)，其中m的取值范圍為“0”到M，M、N由實(shí)際地震資料給定。

根據(jù)褶積定理，將反射系數(shù)對(duì)與地震子波進(jìn)行褶積可得到地震信號(hào):

re(t,m,n,Δt)}+bm,n[w(t,m,n)*

ro(t,m,n,Δt)]}

(10)

用矩陣的形式表示為式(11)。

s=Dwr

(11)

也就是說，地震響應(yīng)是通過奇、偶楔形反射系數(shù)模型和地震子波褶積獲得的奇、偶地震響應(yīng)的線性組合。式(11)中，Dw是奇、偶楔形模型地震響應(yīng)。傳統(tǒng)的基追蹤反演是在褶積模型假設(shè)的前提下，用統(tǒng)計(jì)方法或井震標(biāo)定方法從實(shí)際地震數(shù)據(jù)提取出地震子波w(t)。由于隨機(jī)噪聲干擾、地下介質(zhì)吸收衰減作用的影響，從實(shí)際地層采集的地震數(shù)據(jù)呈現(xiàn)非平穩(wěn)特性，運(yùn)用本文方法逐道提取各時(shí)刻的子波，逐列替代時(shí)不變子波構(gòu)建時(shí)變子波核矩陣，可得：

re(t,m,n,Δt)]+bm,n[w(t,m,n)*

ro(t,m,n,Δt)]}

(12)

反演過程中，目標(biāo)函數(shù)變?yōu)椋?/p>

(13)

采用基追蹤最優(yōu)化準(zhǔn)則從式(12)中求得參數(shù)am,n和bm,n，因?yàn)槭?9)和式(12)采用相同的系數(shù)，因此將am,n、bm,n代入式(9)通過求和重構(gòu)可得到有效的反射系數(shù)。在準(zhǔn)確提取地震時(shí)變子波的前提下，通過基追蹤最優(yōu)化準(zhǔn)則可求得最優(yōu)的反射系數(shù)值。

2 模型測試

2.1 不含噪聲模型

為了檢驗(yàn)時(shí)變子波庫的建立對(duì)于基追蹤反演的重要性，設(shè)計(jì)如圖3(a)所示以1 ms為采樣間隔且時(shí)間長度為1 000 ms的反射系數(shù)序列，其中各反射系數(shù)脈沖的值分別設(shè)定為(120 ms，-0.15)、(165 ms，0.09)、(380 ms，-0.06)、(540 ms，-0.06)、(760 ms，0.13)、(900 ms，0.18)，非平穩(wěn)地震記錄由主頻從60 Hz單調(diào)下降到40 Hz的Ricker子波與設(shè)計(jì)反射系數(shù)褶積合成，如圖3(b)所示。分別采用60 Hz子波字典、50 Hz子波字典、40 Hz子波字典以及由本文方法提取的隨時(shí)間變化子波構(gòu)建的子波字典，對(duì)合成的非平穩(wěn)地震記錄進(jìn)行基追蹤反演，其反演結(jié)果如圖4所示。由圖4(a)～圖4(c)可知，反演所用子波的主頻越接近合成地震記錄所用的子波主頻時(shí)，反演出的反射系數(shù)才更接近實(shí)際，否則反演出的反射系數(shù)的值與實(shí)際相差較遠(yuǎn)，甚至出現(xiàn)虛假的反射系數(shù)。由此可知，針對(duì)非平穩(wěn)地震記錄而言，若采用固定主導(dǎo)頻率的時(shí)不變子波字典用于基追蹤反演，則難以確保反演結(jié)果的準(zhǔn)確度和可靠性，而采用時(shí)變子波字典進(jìn)行基追蹤反演時(shí)能獲得準(zhǔn)確的反射系數(shù)系列(圖4(d))，說明時(shí)變子波的估計(jì)準(zhǔn)確程度對(duì)于基追蹤反演精度有重要影響。

圖3 反射系數(shù)序列及合成地震記錄

圖4 基追蹤反演結(jié)果

2.2 含噪聲模型

為檢驗(yàn)本文方法的抗噪性，對(duì)非平穩(wěn)地震道分別添加SNR=15 dB和SNR=25 dB的隨機(jī)噪聲后,得到如圖5(a)與圖5(d)所示的含噪非平穩(wěn)地震記錄，對(duì)加噪后的非平穩(wěn)地震道采用上述方法處理,可得到含噪非平穩(wěn)地震道的時(shí)頻譜和“時(shí)變子波”時(shí)頻譜。對(duì)圖5(b)和圖5(e)進(jìn)行分析可知，廣義S變換在含噪的情況下依舊具有較好的時(shí)頻聚焦性。利用本文方法分別對(duì)含噪地震數(shù)據(jù)提取時(shí)變子波，進(jìn)而基追蹤反演，反演結(jié)果如圖6所示，將反演出的反射系數(shù)與原始反射系數(shù)作對(duì)比可知，時(shí)變子波基追蹤反演結(jié)果與原始反射系數(shù)結(jié)果近似，由此可以得出本文子波提取及反演方法具有較好的抗噪性。

圖5 含噪地震記錄及時(shí)頻譜

圖6 含噪地震記錄基追蹤反演結(jié)果

3 實(shí)際應(yīng)用

這里所用數(shù)據(jù)為南海某實(shí)際工區(qū)的疊后過井地震數(shù)據(jù)剖面如圖7所示，該地震剖面共有201道，每道采樣點(diǎn)數(shù)為256，截取的時(shí)窗范圍為2 401 ms～2 912 ms，其中測井曲線為密度曲線。

圖7 實(shí)際地震數(shù)據(jù)剖面

選擇圖7中紫色虛線標(biāo)注的第86道地震記錄(圖8(a))進(jìn)行廣義S變換處理后，可以獲得如圖8(b)所示的時(shí)頻譜圖，對(duì)時(shí)頻譜進(jìn)行譜模擬后得到如圖8(c)所示的 “時(shí)變子波”時(shí)頻譜圖。對(duì)比圖8(b)和圖8(c)可知，時(shí)變子波時(shí)頻譜消除了地層反射系數(shù)的影響，與地震波形吻合度較高。

圖8 第86道地震記錄及時(shí)頻譜圖

利用廣義S變換譜模擬方法對(duì)該地震剖面逐道提取時(shí)變子波，其不僅沿時(shí)間軸方向變化，而且還沿著聯(lián)絡(luò)測線的方向變化。在整個(gè)地震剖面中，對(duì)每道地震記錄的所有時(shí)間點(diǎn)提取隨時(shí)間變化的地震子波，而對(duì)于地震剖面而言，提取的地震子波隨時(shí)間變化的同時(shí)，還隨空間變化?？紤]到地震子波在不同道之間的差異，以提取的時(shí)變子波為基礎(chǔ)，對(duì)所估計(jì)的地震子波采用多地震道加權(quán)空間變化處理[17]。為了研究時(shí)不變子波和時(shí)變子波對(duì)反演結(jié)果精度的影響和差異，在此基礎(chǔ)上采用本文方法提取的時(shí)變子波和時(shí)不變子波,用于實(shí)際工區(qū)地震資料的基追蹤反演，反演結(jié)果如圖9所示。對(duì)比圖9(a)與圖9(b)可以看出，與時(shí)變子波基追蹤結(jié)果圖中的綠色方框部分相比，時(shí)變子波基追蹤地震反演結(jié)果的分辨率要高于時(shí)不變子波基追蹤反演結(jié)果的分辨率；時(shí)不變子波基追蹤反演出的反射系數(shù)出現(xiàn)虛假反射，時(shí)變子波基追蹤反演出的反射系數(shù)與原始剖面更為接近；同時(shí)，與時(shí)不變子波基追蹤結(jié)果圖中的綠色方框部分相比，時(shí)變子波基追蹤反演出的反射系數(shù)橫向連續(xù)性更好，究其原因，是因?yàn)楸疚姆椒ㄔ谇笕r(shí)變子波時(shí)還考慮了地震數(shù)據(jù)縱橫向的變化特性。

圖9 不同子波字典基追蹤反演剖面對(duì)比圖

4 結(jié)論

1)本文方法提取的時(shí)變子波符合地震子波在地下傳播時(shí)的頻率衰減特性，利用時(shí)變子波建立的子波庫具有較高的精度，為后續(xù)的基追蹤反演提供了保證。

2)筆者提取時(shí)變子波過程中，沒有考慮子波的相位影響，后續(xù)可以結(jié)合子波相位的估計(jì)來提取時(shí)變子波。

3)實(shí)際地震剖面的反演結(jié)果表明，基于本文方法提取的時(shí)變子波基追蹤地震反演的結(jié)果，與傳統(tǒng)時(shí)不變子波基追蹤反演得到的反演結(jié)果相比，能進(jìn)一步提高反演精度，可用于后期儲(chǔ)層預(yù)測與精細(xì)刻畫。