基于時(shí)頻域譜模擬的時(shí)變子波估計(jì)方法

張漫漫，戴永壽，張亞南，丁進(jìn)杰，王蓉蓉

(中國石油大學(xué)(華東)信息與控制工程學(xué)院,山東東營257061)

地震波在地下介質(zhì)傳播過程中，由于地層的吸收和濾波作用，子波能量逐漸被吸收，高頻成分逐漸衰減，因此待估計(jì)的子波具有動(dòng)態(tài)衰減性。為了滿足油氣勘探的需要，得到高分辨率的地震資料，時(shí)變子波的提取成為地球物理工作者一個(gè)重要的研究方向。

時(shí)變子波估計(jì)方法很多，如胡啟宇[1]和Liang等[2]使用同態(tài)理論和高階統(tǒng)計(jì)量估算時(shí)變和空變地震子波;彭才等[3-4]提出了基于動(dòng)態(tài)褶積模型的子波估計(jì)方法;Baan[5]提出了基于峰值最大化的時(shí)變地震子波估算方法;高靜懷等[6]基于反射地震記錄變子波模型提取地震子波進(jìn)而提高地震記錄分辨率;馮晅等[7]提出分時(shí)窗提取地震子波的方法并應(yīng)用于合成地震記錄中。這些方法大都采用分段處理，即首先對地震記錄在時(shí)間上進(jìn)行分段，假設(shè)每段內(nèi)的地震子波是時(shí)不變的，進(jìn)而在每段內(nèi)提取出一個(gè)時(shí)不變的地震子波。受分段長度的影響，每段提取的平均意義子波與實(shí)際子波勢必存在一定誤差，使得后續(xù)的地震資料處理和解釋結(jié)果不準(zhǔn)確。當(dāng)相鄰地層的地震反射振幅和頻率成分差異比較大時(shí)，分段提取的方法將不能很好的反映相鄰層段中子波的變化。

為了突破地震記錄平穩(wěn)性假設(shè)，解決上述常規(guī)分段子波提取方法存在的局限性和不足，本文考慮到地震記錄中由于地層的吸收濾波作用導(dǎo)致的地震波振幅和高頻成分衰減的情況，提出了一種基于時(shí)頻譜模擬的時(shí)變子波估計(jì)方法。

1 基于時(shí)頻域譜模擬的時(shí)變子波估計(jì)方法

首先對地震記錄做時(shí)頻分析，在二維時(shí)間-頻率域內(nèi)進(jìn)行濾波處理，然后采用譜模擬的方法擬合出每一反射點(diǎn)處的子波振幅譜(這里假設(shè)地震記錄已經(jīng)過零相位化處理)，從而實(shí)現(xiàn)時(shí)變子波的提取，其處理流程如圖1所示。該方法不需要考慮對子波進(jìn)行衰減補(bǔ)償，與分段子波提取方法相比，摒棄了分段平穩(wěn)的假設(shè)，更精細(xì)地刻畫了地震記錄中的非平穩(wěn)成分。

圖1 基于時(shí)頻域譜模擬的時(shí)變子波估計(jì)方法流程

1.1 地震記錄的時(shí)頻分析

實(shí)際地下介質(zhì)是粘彈性的，地震子波在地下介質(zhì)中傳播時(shí)由于波前擴(kuò)散和吸收衰減等效應(yīng)會(huì)使子波能量衰減、頻帶變窄，所以地震子波隨時(shí)間變化，記錄的地震資料是一個(gè)非平穩(wěn)信號(hào)。而時(shí)頻分析是研究非平穩(wěn)信號(hào)的有力工具，它以聯(lián)合時(shí)頻分布的形式來表示信號(hào)的特性，可以較為準(zhǔn)確地定位某一時(shí)刻出現(xiàn)的頻率分量。所以，對地震記錄進(jìn)行時(shí)頻分析能夠充分顯示出其中的非平穩(wěn)特性。

常用的時(shí)頻分析方法有短時(shí)傅里葉變換、Gabor變換、小波變換和S變換等。其中，S變換是在繼承短時(shí)傅里葉變換和小波變換優(yōu)點(diǎn)的同時(shí)，克服了短時(shí)傅里葉變換時(shí)頻分辨率不變以及小波變換的尺度因子與頻率無關(guān)的缺點(diǎn)。S變換采用與頻率有關(guān)的可變高斯窗函數(shù)，其窗函數(shù)定義如下：

(1)

式中：f為頻率?？梢钥闯?，S變換窗函數(shù)的變化趨勢固定不變，不能根據(jù)實(shí)際需要靈活調(diào)整，為此很多學(xué)者對S變換的窗函數(shù)加以改進(jìn)提出了多種廣義S變換[8-10]。常用的廣義S變換采用時(shí)窗寬度隨頻率f呈反比例變化的高斯窗函數(shù)。在高頻段時(shí)窗較窄，可獲得較高的時(shí)間分辨率;在低頻段時(shí)窗較寬，可獲得較高的頻率分辨率。這些廣義S變換對于大多數(shù)非平穩(wěn)信號(hào)能獲得較好的時(shí)頻效果，但仍不能很好滿足非平穩(wěn)地震記錄處理所要求的分辨率。齊春燕等[11]在廣義S變換的基礎(chǔ)上提出一種改進(jìn)的廣義S變換，其窗函數(shù)的表達(dá)式為

(2)

式中：q,p為大于零的調(diào)節(jié)因子。從(2)式可以看出窗函數(shù)的寬度與頻率呈正比例變化，即在低頻處可獲得較高的時(shí)間分辨率，在高頻處獲得較高的頻率分辨率，此性質(zhì)符合地震記錄動(dòng)態(tài)衰減特性。信號(hào)s(t)的時(shí)頻譜為

(3)

利用該方法對地震記錄進(jìn)行時(shí)頻分析，可獲得更高的時(shí)頻分辨率，且具有很好的時(shí)頻聚焦性，能夠有效地分辨出地震資料中頻率成分的變化，有利于更精細(xì)地分析地層結(jié)構(gòu)。此外，該方法還可進(jìn)行無能量損失的反變換，從而準(zhǔn)確重構(gòu)時(shí)間域地震信號(hào)。

1.2 地震記錄時(shí)頻譜的時(shí)頻濾波

由上述可知，改進(jìn)的廣義S變換具有良好的時(shí)頻聚焦性，能夠獲得高質(zhì)量的時(shí)頻譜信息，在壓制噪聲、判斷反射層的位置以及油氣的直接顯示等方面可獲得較好的應(yīng)用。然而受測不準(zhǔn)原理的限制，我們不可能同時(shí)獲得較高的時(shí)間分辨率和頻率分辨率[12]。由文獻(xiàn)[13]可知，對所獲得的非平穩(wěn)地震記錄的時(shí)頻譜直接處理，有一定效果，但不是很理想。為了更好的適應(yīng)薄互層的地震分析和解釋，在保持現(xiàn)有頻率分辨率的同時(shí)盡可能提高時(shí)間分辨率，避免相鄰反射地層之間的影響，引入一種時(shí)頻濾波器[14]對地震記錄的時(shí)頻譜進(jìn)行濾波處理。該濾波器應(yīng)是時(shí)間和頻率的光滑函數(shù)，且具有以下特性：

1) 對某一時(shí)間τ，當(dāng)t=τ時(shí)，H(t,ω)=1;當(dāng)|t-τ|→∞時(shí)，H(t,ω)=0;

2) 對某一時(shí)間ξ，當(dāng)ω=ξ時(shí)，H(t,ω)=1;當(dāng)|ω-ξ|→∞時(shí)，H(t,ω)=0。

即H(t,ω)對空間點(diǎn)(τ,ξ)具有時(shí)頻局域化性質(zhì)。據(jù)此定義如下形式的3參數(shù)(λ,p,v)時(shí)頻濾波器：

(4)

函數(shù)的后一部分對前一部分起衰減控制作用，不同的λ,p對應(yīng)不同的衰減。λ越小，p越大，H(t,ω)越尖銳，支集越小，實(shí)際應(yīng)用中一般λ取負(fù)值，p取正值。取λ=-1,p=1,v=2時(shí)的時(shí)頻濾波器的三維曲面圖如圖2所示。

圖2 時(shí)頻濾波器的三維曲面表示

利用上述濾波器對地震記錄的時(shí)頻譜進(jìn)行濾波：

(5)

f≠0

式中：S(τ,f)為處理后的地震記錄時(shí)頻譜;T為時(shí)頻濾波器的時(shí)窗寬度;tr為反射系數(shù)出現(xiàn)的時(shí)間點(diǎn)，且T/2小于相鄰反射系數(shù)最小間隔的一半。

1.3 基于時(shí)頻域譜模擬的時(shí)變子波提取

得到地震記錄的時(shí)頻譜后，需要從中擬合出子波的時(shí)變振幅譜。常用的方法是在反射系數(shù)是白噪序列的假設(shè)條件下，通過地震記錄的自相關(guān)得到子波振幅譜。然而實(shí)際的反射系數(shù)序列并不一定符合白噪假設(shè)，且自相關(guān)法不能用于時(shí)頻域提取時(shí)變子波，為此本文提出采用譜模擬[15]的方法估計(jì)時(shí)變子波。該方法一般用于反褶積處理中[13,16-17]，是一種比較有效且穩(wěn)健性較強(qiáng)的地震子波模擬方法。

譜模擬假設(shè)地震子波的振幅譜是平滑的，認(rèn)為子波振幅譜類似于雷克子波譜的單峰光滑曲線，對其在頻率域進(jìn)行數(shù)學(xué)建模，數(shù)學(xué)表達(dá)式為

(6)

式中：k是一個(gè)常數(shù);N為階數(shù);an是關(guān)于f的多項(xiàng)式系數(shù)，一般0

基于時(shí)頻域譜模擬的時(shí)變子波估計(jì)方法具體實(shí)現(xiàn)過程：

1) 對地震記錄x(t)進(jìn)行改進(jìn)的廣義S變換，得到相應(yīng)的時(shí)頻譜SNGST(τ,f);

2) 對地震記錄的時(shí)頻譜進(jìn)行時(shí)頻濾波處理，得到處理后的時(shí)頻譜S(τ,f)。

3) 固定時(shí)間T，提取相應(yīng)時(shí)刻的地震記錄頻譜S(T,f)，通過譜模擬的方法擬合出該時(shí)刻的地震子波振幅譜W(T,f)，根據(jù)實(shí)際需要對相位進(jìn)行處理(本文假設(shè)子波已通過零相位化處理)，反傅里葉變換即可得到該時(shí)刻的地震子波。

根據(jù)地震記錄中反射點(diǎn)出現(xiàn)的不同位置，改變時(shí)間τ的取值，可實(shí)現(xiàn)時(shí)變子波的精細(xì)提取。子波提取的最終目的是進(jìn)行反褶積處理，在去除子波影響時(shí)，采用時(shí)變維納濾波法[5]進(jìn)行反褶積處理，從而提高了地震記錄的處理精度。

2 數(shù)值模擬

我們采用含有5個(gè)反射系數(shù)脈沖的稀疏模型進(jìn)行試驗(yàn)，其長度為1000ms，采樣頻率為1ms，如圖3a所示。各反射系數(shù)脈沖的位置和大小分別為(200ms,1.0)，(240ms,0.8)，(520ms,-1.0)，(570ms,0.9)和(840ms,-0.8)，其中包括2個(gè)相鄰的同向反射系數(shù)和2個(gè)相鄰的異向反射系數(shù)。采用主頻分別為60,50，40，35，30Hz且振幅逐漸衰減的零相位雷克子波與反射系數(shù)序列褶積得到逐漸衰減的合成地震記錄，如圖3b所示。

圖3 反射系數(shù)序列(a)和逐漸衰減的合成地震記錄(b)

圖4是對上述合成地震記錄進(jìn)行多種時(shí)頻分析的結(jié)果。從圖4可以看出，與S變換和廣義S變換相比，改進(jìn)的廣義S變換具有更好的時(shí)頻分辨率，且能準(zhǔn)確地反應(yīng)出相鄰2個(gè)地層的頻率成分及其差異。然而，從改進(jìn)的廣義S變換譜可看出，雖然原反射系數(shù)出現(xiàn)的時(shí)間點(diǎn)200，240，520，570，840ms處的能量最大，但在鄰域內(nèi)仍存在頻率成分，這會(huì)對后續(xù)子波估計(jì)進(jìn)而去除子波成分有一定的影響。圖4d是對改進(jìn)的廣義S變換譜濾波后的結(jié)果，可以看出，在保持頻率分辨率的同時(shí)，獲得了更好的時(shí)間分辨率。

從處理后的地震記錄時(shí)頻譜中，分別提取每個(gè)時(shí)刻的地震子波。圖5分別為200，520，570，650ms時(shí)的原始理論子波的振幅譜、地震記錄的振幅譜以及模擬出的子波振幅譜對比結(jié)果。圖6分別為200，520，570ms處的原始理論子波與估計(jì)子波時(shí)域?qū)Ρ冉Y(jié)果。

圖4 S變換(a)、廣義S變換(b)、改進(jìn)的廣義S變換(c)以及濾波后的時(shí)頻譜(d)

由圖5和圖6可以看出：①模擬出的子波振幅譜與原始子波振幅譜基本吻合，反變換到時(shí)域亦是如此，從而證明譜模擬提取時(shí)變子波振幅譜的準(zhǔn)確性;②對比520ms和570ms處的模擬結(jié)果，相鄰地層之間雖然存在相互影響，但本文提出的方法仍能準(zhǔn)確的提取出每一反射點(diǎn)處的子波;③隨著時(shí)間的增加，擬合出的子波其主頻逐漸降低且能量減少，能夠反應(yīng)出子波在傳播過程中的動(dòng)態(tài)衰減特性;④提取反射點(diǎn)以外時(shí)間點(diǎn)處的子波(如圖5d)，雖然能得到一個(gè)振幅譜，但其主頻較原始子波有一定誤差，且能量很低。

用得到的時(shí)變子波對地震記錄進(jìn)行動(dòng)態(tài)反褶積處理，結(jié)果如圖7所示。可以看出，用時(shí)頻域譜模擬提取的時(shí)變子波進(jìn)行反褶積處理后的地震記錄，消弱了子波的影響。對比圖7中濾波前、后的處理效果，可以明顯看出，濾波處理后的反褶積結(jié)果在反射點(diǎn)處更加接近脈沖。

圖5 200ms(a)，520ms(b)，570ms(c)，650ms(d)時(shí)理論子波、地震記錄和擬合子波的振幅譜對比

圖6 200ms(a)，520ms(b)，570ms(c)時(shí)的原始子波與估計(jì)子波時(shí)域?qū)Ρ?/p>

圖7 用估計(jì)子波作動(dòng)態(tài)反褶積處理后的地震記錄a 原始合成記錄; b 處理后記錄(無濾波); c 處理后記錄(加濾波)

若采用分段處理的方法提取子波，200和240ms,520和570ms很可能被分到一個(gè)時(shí)段內(nèi)。圖8是200和240ms所在段內(nèi)提取的平均子波，可以看出，該子波與兩個(gè)時(shí)刻的原始子波均有一定的誤差，而用此子波進(jìn)行反褶積處理結(jié)果必然不準(zhǔn)確。

為了進(jìn)一步說明問題，合成一個(gè)比較接近實(shí)際的非平穩(wěn)地震記錄。隨機(jī)生成一個(gè)稀疏反射系數(shù)利用本文方法處理該非平穩(wěn)地震記錄，圖10a顯示了對地震記錄進(jìn)行改進(jìn)的廣義S變換得到的時(shí)頻譜，圖10b是時(shí)頻濾波處理后的結(jié)果，可見，處理后的時(shí)頻譜在頻率分辨率不變的基礎(chǔ)上，其時(shí)間分辨率大大提高，這有利于后續(xù)的地震資料處理。僅以176，378，577ms為例，圖11a為擬合出的子波振幅譜，圖11b為反變換回時(shí)域后的歸一化子波。從3個(gè)時(shí)刻的圖形對比不難看出，隨著時(shí)間的增加，子波主頻逐漸降低，頻帶變窄，相應(yīng)時(shí)域波形變寬，說明本文方法能夠很好地反映地震記錄中的非平穩(wěn)成分，準(zhǔn)確精細(xì)地提取出時(shí)變子波。由于378ms處的反射系數(shù)較大，模擬的378ms處的子波振幅譜的幅值比176ms處的大，但其主頻仍是降低的。

序列(圖9a)，對應(yīng)的非平穩(wěn)地震記錄為圖9b。

圖8 分段提取子波對比

圖9 反射系數(shù)序列(a)和非平穩(wěn)合成地震記錄(b)

圖10 對非平穩(wěn)合成地震記錄進(jìn)行廣義S變換(a)及濾波后的結(jié)果(b)

圖11 不同時(shí)刻子波提取對比a 子波振幅譜; b 時(shí)域子波

提取出每一反射點(diǎn)處的子波后，將其用于反褶積處理，得到提高分辨率后的地震記錄(圖12)，從處理前、后地震記錄的對比可以很明顯的看出，反褶積處理后的地震記錄上在反射點(diǎn)處更接近脈沖序列，從而說明本文方法能夠準(zhǔn)確提取時(shí)變子波，有效提高地震資料分辨率。

圖12 用提取子波作反褶積處理前(a)、后(b)地震記錄對比

3 結(jié)束語

常用的時(shí)變子波提取方法大都采用分段提取，當(dāng)相鄰層段的子波屬性存在較大變化時(shí)，分段方法不能有效地反應(yīng)這一變化。為此，我們研究提出了一種基于時(shí)頻域譜模擬的時(shí)變子波提取方法。理論分析和數(shù)值模擬結(jié)果表明，該方法有效刻畫了地震記錄中的非平穩(wěn)成分，且不需要對反射系數(shù)序列進(jìn)行白噪假設(shè);時(shí)頻濾波器的引入大大提高了時(shí)間分辨率，能夠更精細(xì)地提取出時(shí)變子波;在提高地震資料分辨率方面，該方法不需要研究地震記錄的動(dòng)態(tài)衰減機(jī)制，直接對地震資料進(jìn)行處理，避免了誤差累積。

本文提出的方法在實(shí)際應(yīng)用中的穩(wěn)定性和普適性等問題，有待后續(xù)在對實(shí)際地震資料的處理中作進(jìn)一步探討。

參 考 文 獻(xiàn)

[1] 胡啟宇.同態(tài)反褶積的一種可能途徑[J].石油物探,1984,23(2):109-111

Hu Q Y.A possible way for hommorphic deconvolution [J].Geophysical Prospecting for Petroleum,1984,23(2):109-111

[2] Liang G H,Cai X P,Li Q Y.Using high-order cumulants to extrapolate spatially variant seismic wavelets [J].Geophysics,2002,67(6):1869-1876

[3] 彭才,朱仕軍,黃中玉,等.基于動(dòng)態(tài)褶積模型的動(dòng)態(tài)子波估計(jì)[J].石油物探,2007,46(4):224-328

Peng C,Zhu S J,Huang Z Y,et al.Dynamic wavelet estimation based on the dynamic convolution model [J].Geophysical Prospecting for Petroleum,2007,46(4):224-328

[4] 彭才,曾濤,常智,朱仕軍.混合相位動(dòng)態(tài)子波估計(jì)[J].油氣地球物理,2008,6(1):21-24

Peng C,Zeng T,Chang Z,et al.Estimation of mixed-phase nonstationary wavelet [J].Petroleum Geophysics,2008,6(1):21-24

[5] Mirko van der B.Time-varying wavelet estimation and deconvolution by kurtosis maximization [J].Geophysics,2008,73(2):11-18

[6] 高靜懷,汪玲玲,趙偉.基于反射地震記錄變子波模型提高地震記錄分辨率[J].地球物理學(xué)報(bào),2009,52(5):1289-1300

Gao J H,Wang L L,Zhao W.Enhancing resolution of seismic trances based on the changing wavelet model of the seismogram [J].Chinese Journal of Geophysics,2009,52(5):1289-1300

[7] 馮晅,劉財(cái),楊寶俊,等.分時(shí)窗提取地震子波及在合成地震記錄中的應(yīng)用[J].地球物理學(xué)進(jìn)展,2002,17(1):71-77

Feng X,Liu C,Yang B J.The extractive method of seismic wavelet in different time window and the application in synthetic seismogram [J].Progress in Geophysics,2002,17(1):71-77

[8] 高靜懷,陳文超,李幼銘，等.廣義S變換與薄互層地震響應(yīng)分析[J].地球物理學(xué)報(bào),2003,46(4):526-532

Gao J H,Chen W C,Li Y M,et al.Generalized S transform and seismic response analysis of thin interbeds [J].Chinese Journal of Geophysics,2003,46(4):526-532

[9] Stockwell R G.A basis for ef digital signal processing of the S-transform [J].Digital Signal Processing,2006:1-22

[10] Ervin S, Djurovi L, Jin J. A window width optimized S-transform[J]. EURASIP Journal on Advances in Signal Processing,28(1):1-13

[11] 齊春艷,李彥鵬,彭繼新,等.一種改進(jìn)的廣義S變換[J].石油地球物理勘探，2010,45(2):215-218

Qi C Y,Li Y P,Peng J X,et al.An improved generalized S-transform [J].Oil Geophysical Prospecting,2010,45(2):215-218

[12] 張賢達(dá).現(xiàn)代信號(hào)處理[M].第2版.北京：清華大學(xué)出版社,2002:360-362

Zhang X D.Modern signal processing [M].2nd ed.Beijing:Tsinghua University Press,2002:360-362

[13] Tian J,Song W,Yang F.Enhancing the resolution of seismic data based on the generalized S-transform [J].Petroleum Science,2009,6:153-157

[14] 高靜懷,張兵.基于時(shí)頻濾波的吸收衰減參數(shù)估算[J].石油地球物理勘探,2012,47(6):931-936

Gao J H,Zhang B.Seismic attenuation parameter estimation based on the time-frequency filtering [J].Oil Geophysical Prospecting,2012,47(6):931-936

[15] Rosa A L，Ulrych T J.Processing via spectral modeling [J].Geophysics.1991,56(8):1244-1251

[16] 孫成禹.譜模擬方法及其在提高地震資料分辨率中的應(yīng)用[J].石油地球物理勘探,2000,35(1):27-35

Sun C Y.Spectrum modeling method and its application to seismic resolution improvement [J].Oil Geophysical Prospecting,2000,35(1):27-35

[17] 唐博文,趙波,吳艷輝,等.一種實(shí)現(xiàn)譜模擬反褶積的新途徑[J].石油地球物理勘探,2010,45(增刊1):66-70

Tang B W,Zhao B,Wu Y H,et al.A new method for spectral-modeled deconvolution [J].Oil Geophysical Prospecting,2010,45(S1):66-70