楊學(xué)亭,劉 財(cái),劉 洋,羅 騰,周 寅,張 鵬,李繼龍

(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春130026)

地震波在地下介質(zhì)中傳播時(shí)，由于地層的非彈性性質(zhì)，使其能量有一定的衰減損耗。國(guó)內(nèi)外專家學(xué)者已經(jīng)對(duì)地震波衰減機(jī)理和補(bǔ)償方法進(jìn)行了大量的研究。應(yīng)用比較早的有譜白化處理和反Q濾波。譜白化是通過(guò)展平信號(hào)的振幅譜來(lái)實(shí)現(xiàn)對(duì)衰減頻率的補(bǔ)償，可分別在時(shí)間域和頻率域進(jìn)行。譜白化處理的特點(diǎn)是只對(duì)信號(hào)的振幅譜進(jìn)行處理而不改變信號(hào)的相位譜。反Q濾波最早由Hale提出[1]。Hargreave等[2]提出與Stolt偏移原理類似的相移校正方法，采用級(jí)聯(lián)式常Q補(bǔ)償?shù)姆绞介_(kāi)時(shí)窗逐級(jí)向下計(jì)算。Wang提出一種基于波場(chǎng)延拓理論的反Q濾波方法[3]，并于2006年對(duì)該方法作了改進(jìn)，推廣到Q隨著地層深度連續(xù)變化的情況[4]。

大地濾波算子是時(shí)間、頻率和品質(zhì)因子的函數(shù)，因此也可以在時(shí)頻域內(nèi)進(jìn)行地震波能量衰減補(bǔ)償。白樺等[5]提出了時(shí)頻域的衰減補(bǔ)償方法，將短時(shí)傅里葉變換(STFT)應(yīng)用到衰減補(bǔ)償中。針對(duì)STFT時(shí)窗寬度固定，不能根據(jù)非平穩(wěn)地震信號(hào)在不同時(shí)刻的不同特征調(diào)節(jié)時(shí)窗寬度的問(wèn)題，李鯤鵬等[6]提出基于小波包分解的地震波能量衰減補(bǔ)償方法。劉喜武等[7]提出了基于廣義S變換的地震波能量衰減補(bǔ)償方法。Braga等[8]提出基于一維連續(xù)小波變換的小波域內(nèi)地震波能量衰減補(bǔ)償方法，取得了很好的效果，文中將尺度和頻率的關(guān)系定義為s=f0/f(f0是小波的中心頻率)。小波變換是一種時(shí)間-尺度分析方法，每一個(gè)尺度對(duì)應(yīng)著一個(gè)頻帶而非單一的頻率，尺度和頻率之間沒(méi)有嚴(yán)格的數(shù)學(xué)關(guān)系。

本文利用Sinha等[9]提出的基于一維連續(xù)小波變換的時(shí)頻譜(TFCWT)和基于Kolsky衰減模型的大地濾波算子，在時(shí)頻域內(nèi)實(shí)現(xiàn)地震波能量的衰減補(bǔ)償。對(duì)比文獻(xiàn)[8]提出的小波域(時(shí)間-尺度域)內(nèi)的衰減補(bǔ)償方法，基于TFCWT的地震波能量衰減補(bǔ)償方法能夠更好地恢復(fù)地震波的振幅和相位，提高地震資料的分辨率。

1 理論基礎(chǔ)

1.1 基于一維連續(xù)小波變換的時(shí)頻譜

短時(shí)傅里葉變換(STFT)和Gabor變換都是時(shí)頻分析常用的方法，但是這些方法的時(shí)頻窗口是固定的，窗口的大小直接決定了時(shí)頻分辨率。一旦窗口函數(shù)選定，其時(shí)頻分辨率就已經(jīng)確定，不會(huì)隨著時(shí)間和頻率而改變，因此不能滿足對(duì)時(shí)頻分辨率的要求。

相反，小波變換不需要預(yù)先設(shè)定窗口的寬度，它利用平移參數(shù)τ和尺度參數(shù)s得到時(shí)間-尺度譜。小波變換具有很好的時(shí)間和頻率局部分析的能力。一維連續(xù)小波變換公式如下:

設(shè)信號(hào)f(t)∈L2(R)，ψ(t)為小波基函數(shù)

(1)

式中:s為尺度參數(shù);τ為平移參數(shù)。且小波基函數(shù)需要滿足可容性條件

(2)

式中:Ψ(ω)是ψ(t)的傅里葉變換;Cψ對(duì)每一個(gè)小波基函數(shù)是一個(gè)常數(shù)。則滿足相容性條件的一維連續(xù)小波變換為

(3)

小波函數(shù)的選取是應(yīng)用小波變換的一個(gè)重要問(wèn)題，在此選用的是Morlet小波

(4)

Morlet小波是一種復(fù)小波基函數(shù)，利用Morlet小波函數(shù)可以將信號(hào)的振幅和相位分開(kāi)，有利于非平穩(wěn)信號(hào)的分析處理[11]。Morlet小波作為一種分析小波，已經(jīng)被認(rèn)為是分析非彈性介質(zhì)中非平穩(wěn)地震信號(hào)最合適的小波。Kulesh等[12]將Morlet小波應(yīng)用于天然地震學(xué)中確定表面波的傳播參數(shù)，群速度、相速度和衰減系數(shù)等。Reine等[13]將Morlet小波應(yīng)用于反射波地震學(xué)。

根據(jù)Addison[11]，一維連續(xù)小波變換在頻率域內(nèi)的計(jì)算公式為

(5)

對(duì)尺度參數(shù)s的選擇，采用Christopher等[14]提出的方式，經(jīng)過(guò)頻率域內(nèi)的小波變換得到時(shí)間-尺度譜。

Hlawatsch等[15]提出利用關(guān)系式s=f0/f(f0是小波中心頻率)可以將時(shí)間-尺度譜轉(zhuǎn)化為時(shí)頻譜。Sinha等[9]利用該數(shù)學(xué)關(guān)系式實(shí)現(xiàn)了時(shí)間-尺度譜到時(shí)頻譜的轉(zhuǎn)換，其結(jié)果顯示在高頻部分，時(shí)頻分析并未取得很好效果，因此提出一種新的方法，對(duì)一維連續(xù)小波變換的反變換做傅里葉變換，即可得到時(shí)頻譜，具體公式推導(dǎo)如下:

首先，根據(jù)Daubechies[16]給出的反小波變換的公式

(6)

對(duì)式(6)進(jìn)行傅里葉變換得到:

(7)

一維連續(xù)小波變換中，在任何尺度s和時(shí)間t上，窗口面積ΔtΔω保持不變，即時(shí)頻分辨率相互制約，但是尺度參數(shù)s和平移參數(shù)τ是相互獨(dú)立的[9]，因此可以去掉對(duì)平移參數(shù)τ的積分，然后把F(ω)替換為F(ω,τ)，得到:

(8)

由方程(8)，就可以利用一維連續(xù)小波變換得到時(shí)頻譜，稱之為T(mén)FCWT。對(duì)信號(hào)的重構(gòu)可以分兩步進(jìn)行:首先對(duì)F(ω,τ)按照時(shí)間方向進(jìn)行疊加求和，然后對(duì)疊加和做反傅里葉變換。

這里通過(guò)對(duì)模擬信號(hào)的處理，實(shí)現(xiàn)基于TFCWT的時(shí)頻分析，并與Gabor變換方法對(duì)比分析。選取chirp信號(hào)，信號(hào)的頻率以對(duì)數(shù)的方式增加，頻率范圍是20～120Hz，得到如圖1所示的時(shí)頻譜。

圖1 chirp信號(hào)的時(shí)頻譜a chirp信號(hào); b 基于Gabor變換的時(shí)頻譜; c基于TFCWT的時(shí)頻譜

從圖1中可以看出，在高頻部分，基于TFCWT方法更有優(yōu)勢(shì)。由于一維連續(xù)小波變換可以利用尺度參數(shù)和平移參數(shù)的相互調(diào)節(jié)，使其具有更好的局部時(shí)頻分析能力。

通過(guò)對(duì)非平穩(wěn)信號(hào)做基于TFCWT的時(shí)頻譜，然后做反變換重構(gòu)原信號(hào)，最后做誤差分析(圖2)，以驗(yàn)證方法的正確性。

分析圖2可以看出,TFCWT時(shí)頻分析方法對(duì)非平穩(wěn)信號(hào)具有很好的分析能力。通過(guò)求取原信號(hào)和重構(gòu)信號(hào)的誤差，驗(yàn)證了TFCWT時(shí)頻分析方法以及信號(hào)重構(gòu)方法的正確性。

1.2 線性衰減模型和衰減地震道的構(gòu)建

在勘探地震學(xué)中，常用由Kolsky[17]提出的衰減模型，假設(shè)衰減因子λ是頻率f的線性函數(shù)(Q≥1)。根據(jù)Kolsky的模型，衰減因子和相速度分別為

(9)

式中：vr和fr分別表示參考相速度和參考頻率。

Wang[18]提出了一種衰減模型，其基本理論是在均勻粘彈性介質(zhì)中，波數(shù)k變成了一個(gè)關(guān)于相速度v，衰減系數(shù)λ和頻率f的復(fù)數(shù)k(f)=2πf/v+iλ。經(jīng)過(guò)測(cè)試發(fā)現(xiàn)，把fr設(shè)置為大值(比如500Hz)更有利于對(duì)擴(kuò)散的恢復(fù)。將(9)代入到復(fù)波數(shù)k(f)

(10)

通過(guò)在頻率域內(nèi)求解Helmoltz方程的平面波解，可以推導(dǎo)出頻率域內(nèi)衰減信號(hào)S(f)和原信號(hào)S0(f)之間的關(guān)系為

(11)

式中:ζ表示波前面沿著射線路徑從震源到檢波器的距離。將(10)代入(11)中，因?yàn)樵诳碧降卣饘W(xué)中常用的Q值范圍是20～200，因此可以忽略含有1/Q2的項(xiàng)，化簡(jiǎn)得到

(12)

圖2 非平穩(wěn)信號(hào)的時(shí)頻譜及其信號(hào)重構(gòu)a 非平穩(wěn)信號(hào); b 基于TFCWT的時(shí)頻譜; c 重構(gòu)信號(hào); d 誤差分析

(12)式中并對(duì)所有的頻率進(jìn)行疊加，得到時(shí)間域內(nèi)的衰減地震道為

(13)

這就得到了基于Kolsky衰減模型的衰減地震道。濾波算子G是時(shí)間和頻率的函數(shù)，由此可以實(shí)現(xiàn)在時(shí)頻域內(nèi)對(duì)衰減信號(hào)進(jìn)行補(bǔ)償。

1.3 時(shí)頻域內(nèi)地層衰減補(bǔ)償

定義反Q濾波算子為

(14)

(15)

2 理論模型測(cè)試分析

2.1 合成地震信號(hào)的衰減補(bǔ)償

采用主頻f0=50Hz的Ricker子波合成地震信號(hào)(圖3a)，時(shí)間采樣間隔為dt=0.004s，最大延續(xù)時(shí)間tmax=1.2s;品質(zhì)因子Q=80，采用Kolsky衰減模型合成衰減地震信號(hào)(圖3b)。分別在基于一維連續(xù)小波變換的小波域和TFCWT時(shí)頻域內(nèi)做地震波能量衰減補(bǔ)償(圖4)。圖4a，圖4c和圖4e 是合成地震信號(hào)在小波域(時(shí)間-尺度域)內(nèi)的衰減補(bǔ)償，可見(jiàn)，補(bǔ)償信號(hào)在振幅和相位上都得到較好恢復(fù);圖4b，圖4d和圖4f是合成地震信號(hào)在TFCWT時(shí)頻域內(nèi)進(jìn)行的能量衰減補(bǔ)償，可見(jiàn)，振幅和相位的補(bǔ)償效果都比在小波域內(nèi)好。原因可能是小波域內(nèi)衰減補(bǔ)償定義尺度和頻率的關(guān)系為s=f0/f，這并不是它們之間嚴(yán)格的數(shù)學(xué)關(guān)系。

圖3 合成地震信號(hào)(a)及其基于Kolsky模型的衰減地震信號(hào)(b)

圖4 合成地震信號(hào)在小波域和TFCWT時(shí)頻域內(nèi)的衰減補(bǔ)償a 衰減信號(hào)的時(shí)間-尺度譜; b 衰減信號(hào)基于TFCWT的時(shí)頻譜; c 衰減補(bǔ)償后信號(hào)的時(shí)間-尺度譜; d 衰減補(bǔ)償后信號(hào)的時(shí)頻譜; e 小波域內(nèi)衰減補(bǔ)償后的信號(hào); f 時(shí)頻域內(nèi)衰減補(bǔ)償后的信號(hào)

2.2 模擬地震記錄的衰減補(bǔ)償

采用主頻f0=50Hz的最小相位子波，白噪的反射系數(shù)序列，通過(guò)正演獲得模擬地震記錄，時(shí)間采樣間隔為dt=0.004s，最大延續(xù)時(shí)間tmax=1.2s;品質(zhì)因子Q=80，基于Kolsky模型合成衰減記錄，在TFCWT時(shí)頻域內(nèi)做衰減補(bǔ)償(圖5)。

對(duì)比圖5a至圖5c可以看出，基于一維連續(xù)小波變換的時(shí)頻分析方法能夠分析非平穩(wěn)的地震信號(hào)，在TFCWT時(shí)頻域內(nèi)做地層衰減補(bǔ)償能夠很好地補(bǔ)償深層信號(hào)的振幅和相位，大的反射層位都得到了很好的恢復(fù)，分辨率明顯提高。

圖5 模擬地震記錄在TFCWT時(shí)頻域內(nèi)的衰減補(bǔ)償a 模擬地震記錄及其時(shí)頻譜; b 衰減地震記錄及其時(shí)頻譜; c衰減補(bǔ)償后的時(shí)頻譜及重構(gòu)信號(hào)

3 實(shí)際地震資料試處理

利用某地區(qū)的實(shí)際二維疊后地震資料進(jìn)行本文方法的試處理。圖6a是選取的其中第212～301道(1.6～2.2s)實(shí)際地震數(shù)據(jù)，該數(shù)據(jù)的時(shí)間采樣間隔為1ms;圖6b是利用基于TFCWT的地震波能量補(bǔ)償方法進(jìn)行補(bǔ)償處理后的結(jié)果。

對(duì)比圖6a和圖6b可以看出，地震波能量得到了很好的補(bǔ)償，補(bǔ)償后的剖面上反射波同相軸更加清晰，深層分辨率得到了明顯的改善。實(shí)際資料試處理結(jié)果說(shuō)明基于TFCWT的地震波能量補(bǔ)償方法對(duì)大地濾波衰減具有很好的補(bǔ)償和恢復(fù)作用，能夠較好地提高地震資料的分辨率。

圖6 某地區(qū)二維疊后地震資料(a)及其基于TFCWT時(shí)頻域內(nèi)地震波能量補(bǔ)償后的結(jié)果(b)

4 結(jié)束語(yǔ)

一維連續(xù)小波變換具有很好的局部分析能力，利用一維連續(xù)小波變換可以得到時(shí)間-尺度譜，每個(gè)尺度表示的是一個(gè)頻帶而不是單個(gè)頻率。通過(guò)對(duì)一維連續(xù)小波變換的反變換做傅里葉變換得到時(shí)頻譜，就是基于一維連續(xù)小波變換的時(shí)頻分析方法。Kolsky模型的大地濾波算子是時(shí)間和頻率的函數(shù)，因此可以實(shí)現(xiàn)在時(shí)頻域內(nèi)做地震波能量的衰減補(bǔ)償。理論模型和實(shí)際地震資料試算的結(jié)果表明，基于TFCWT的地震波能量補(bǔ)償方法能夠很好地補(bǔ)償深層衰減地震信號(hào)，提高地震資料分辨率，其應(yīng)用效果優(yōu)于小波域衰減補(bǔ)償方法。

參 考 文 獻(xiàn)

[1] Hale D.An inverse Q filter[J].SEP,1981,26(22):231-243

[2] Hargreaves N D,Calvert A J.Inverse Q filtering by Fourier transform[J].Geophysics,1991,56(4):519-527

[3] Wang Y H.A stable and efficient approach of inverse Q filtering[J].Geophysics,2002,67(2):657-664

[4] Wang Y H.Inverse Q-filter for seismic resolution enhancement[J].Geophysics,2006,71(3):51-61

[5] 白樺,李鯤鵬.基于時(shí)頻分析的地層吸收補(bǔ)償[J].石油地球物理勘探,1999,34(6):642-648

Bai H,Li K P.Stratigraphic absorption compensation based on the time-frequency analysis[J].Oil Geophysical prospecting,1999,34(6):642-648

[6] 李鯤鵬,李衍達(dá),張學(xué)工.基于小波包分解的地層吸收補(bǔ)償[J].地球物理學(xué)報(bào),2000,43(4):542-549

Li K P,Li Y D,Zhang X G.A method to compensate earth filtering based on wavelet packet[J].Chinese Journal of Geophysics,2000,43(4):542-549

[7] 劉喜武,年靜波,劉洪.基于廣義S變換的吸收衰減補(bǔ)償方法[J].石油物探,2006,45(1):9-14

Liu X W,Nian J B,Liu H.The absorption attenuation compensation method based on the generalized S transform[J].Geophysical Prospecting for Petroleum,2006,45(1):9-14

[8] Braga I L S,Moraes F S.High-resolution gathers by inverse Q filtering in the wavelet domain[J].Geophysics,2013,78(2):53-61

[9] Sinha S,Routh P S,Anno P D,et al.Spectral decomposition of seismic data with continuous-wavelet transform[J].Geophysics,2005,70(6):19-25

[10] 許惠平,周云軒，孫運(yùn)生，等.小波和球面小波技術(shù)及其在位場(chǎng)分析中的應(yīng)用[M].北京:科學(xué)出版社，2004:3-40

Xu H P，Zhou Y X,Sun Y S,et al.Wavelet and spherical wavelet theories and their applications in potential field[M].Beijing:Science Press,2004:3-40

[11] Addison P S.The illustrated wavelet transform handbook:Introductory theory and applications in science,engineering,medicine and finance[M].London:Institute of Physics Publishing,2002:6-63

[12] Kulesh M,Holschneider M,Diallo M S,et al.Modeling of wave dispersion using continuous wavelet transforms[J].Pure and Applied Geophysics,2005,162(5):843-855

[13] Reine C,Baan M van der,Clark R.The robustness of seismic attenuation measurements using fixed and variable window time-frequency transforms[J].Geophysics,2009,74(2):123-135

[14] Torrence C,Compo G P.A practical guide to wavelet analysis[J].Bulletin of the American Meteorological Society,1998,79(1):61-78

[15] Hlawatsch F,Boudreaux-Bartels G F.Linear and quadratic time-frequency signal representations[J].IEEE Signal Processing,1992,9(2):21-67

[16] Daubechies I.Ten lectures on wavelets[M].SIAM:Society of Industrial and Applied Mathematics,1992:1-355

[17] Kolsky H.The propagation of stress pulses in viscoelastic solids[J].Philosophical Magazine,1956,1(8):693-710

[18] Wang Y H.Seismic inverse Q filtering[M].New York:Wiley-Black well,2008:1-248