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

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

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

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

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

1 基于時頻域譜模擬的時變子波估計方法

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

圖1 基于時頻域譜模擬的時變子波估計方法流程

1.1 地震記錄的時頻分析

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

常用的時頻分析方法有短時傅里葉變換、Gabor變換、小波變換和S變換等。其中，S變換是在繼承短時傅里葉變換和小波變換優(yōu)點(diǎn)的同時，克服了短時傅里葉變換時頻分辨率不變以及小波變換的尺度因子與頻率無關(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變換采用時窗寬度隨頻率f呈反比例變化的高斯窗函數(shù)。在高頻段時窗較窄，可獲得較高的時間分辨率;在低頻段時窗較寬，可獲得較高的頻率分辨率。這些廣義S變換對于大多數(shù)非平穩(wěn)信號能獲得較好的時頻效果，但仍不能很好滿足非平穩(wěn)地震記錄處理所要求的分辨率。齊春燕等[11]在廣義S變換的基礎(chǔ)上提出一種改進(jìn)的廣義S變換，其窗函數(shù)的表達(dá)式為

(2)

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

(3)

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

1.2 地震記錄時頻譜的時頻濾波

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

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

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

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

(4)

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

圖2 時頻濾波器的三維曲面表示

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

(5)

f≠0

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

1.3 基于時頻域譜模擬的時變子波提取

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

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

(6)

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

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

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

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

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

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

2 數(shù)值模擬

我們采用含有5個反射系數(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個相鄰的同向反射系數(shù)和2個相鄰的異向反射系數(shù)。采用主頻分別為60,50，40，35，30Hz且振幅逐漸衰減的零相位雷克子波與反射系數(shù)序列褶積得到逐漸衰減的合成地震記錄，如圖3b所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

常用的時變子波提取方法大都采用分段提取，當(dāng)相鄰層段的子波屬性存在較大變化時，分段方法不能有效地反應(yīng)這一變化。為此，我們研究提出了一種基于時頻域譜模擬的時變子波提取方法。理論分析和數(shù)值模擬結(jié)果表明，該方法有效刻畫了地震記錄中的非平穩(wěn)成分，且不需要對反射系數(shù)序列進(jìn)行白噪假設(shè);時頻濾波器的引入大大提高了時間分辨率，能夠更精細(xì)地提取出時變子波;在提高地震資料分辨率方面，該方法不需要研究地震記錄的動態(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] 彭才,朱仕軍,黃中玉,等.基于動態(tài)褶積模型的動態(tài)子波估計[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] 彭才,曾濤,常智,朱仕軍.混合相位動態(tài)子波估計[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é)報,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] 馮晅,劉財,楊寶俊,等.分時窗提取地震子波及在合成地震記錄中的應(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é)報,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)代信號處理[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ù)估算[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