時(shí)頻分析在地震數(shù)據(jù)處理中的應(yīng)用

時(shí)頻分析在地震數(shù)據(jù)處理中的應(yīng)用摘要：經(jīng)典的傅立葉分析只適用于分析平穩(wěn)信號(hào)，而不適用于非平穩(wěn)信號(hào)。為了分析非平穩(wěn)信號(hào)，我們采用時(shí)頻分析方法。時(shí)頻分析能夠清楚的揭示信號(hào)的時(shí)變頻譜特征，是對(duì)時(shí)變、非平穩(wěn)信號(hào)進(jìn)行分析與處理的有力工具。本文介紹了時(shí)頻分析中常用的STFT和Wigner-切但分布和S變換的特點(diǎn)，并利用Matlab對(duì)一個(gè)地震信號(hào)進(jìn)行時(shí)頻分析的實(shí)現(xiàn)。關(guān)鍵詞：時(shí)頻分析；短時(shí)傅立葉變換（STFT）;S變換；Wigner-3帕分布引言在傳統(tǒng)的信號(hào)處理領(lǐng)域，基于Fourier變換的信號(hào)頻域表示及其能量的頻域分布揭示了信號(hào)在頻域的特征，它們?cè)趥鹘y(tǒng)的信號(hào)分析與處理的發(fā)展史上發(fā)揮了極其重要的作用。但是，F(xiàn)ourier變換是一種整體變換，即對(duì)信號(hào)的表征要么完全在時(shí)域，要么完全在頻域，作為頻域表示的功率譜并不能告訴我們其中某種頻率分量出現(xiàn)在什么時(shí)候及其變化情況。然而，在許多實(shí)際應(yīng)用場(chǎng)合，信號(hào)是非平穩(wěn)的，其統(tǒng)計(jì)量（如相關(guān)函數(shù)、功率譜等）是時(shí)變函數(shù)。這時(shí)，只了解信號(hào)在時(shí)域或頻域的全局特性是遠(yuǎn)遠(yuǎn)不夠的，最希望得到的乃是信號(hào)頻譜隨時(shí)間變化的情況。為此，需要使用時(shí)間和頻率的聯(lián)合函數(shù)來表示信號(hào)，這種表示簡(jiǎn)稱為信號(hào)的時(shí)頻表示。時(shí)頻分析方法旨在通過構(gòu)造一種時(shí)間和頻率的密度函數(shù)，將一個(gè)一維的時(shí)間信號(hào)以二維的時(shí)間一頻率函數(shù)形式表示出來，以揭示信號(hào)中所包含的頻率分量及其隨時(shí)間的變化特性。這使我們不但能夠同時(shí)掌握非平穩(wěn)信號(hào)的時(shí)域及頻域信息，而且可以清楚地了解非平穩(wěn)信號(hào)的頻率是如何隨時(shí)間變化的。通過時(shí)頻分析方法技術(shù)對(duì)地震信號(hào)進(jìn)行分析處理可在獲得地震信號(hào)的瞬時(shí)頻率、瞬時(shí)相位、瞬時(shí)振幅等瞬時(shí)參數(shù)的同時(shí)獲得時(shí)頻譜圖等重要時(shí)頻域信息，實(shí)現(xiàn)對(duì)地震信號(hào)的邊緣檢測(cè)、屬性提取等。時(shí)頻分析是非平穩(wěn)信號(hào)分析處理領(lǐng)域的重要方法，時(shí)頻分布的基本任務(wù)是建立一個(gè)函數(shù)，要求這個(gè)函數(shù)能夠同時(shí)用時(shí)間和頻率來描述信號(hào)的能量密度。如果有了這樣的一個(gè)分布，就可以計(jì)算某一確定的頻率和時(shí)間范圍內(nèi)能量的百分率、計(jì)算某一特定時(shí)刻的頻率密度、計(jì)算該分布的整體和局部的各階矩。即尋找一個(gè)聯(lián)合密度函數(shù)P（tf）,使「（匕力=在時(shí)間t和頻率十的強(qiáng)度，或者P（t,f）tf在時(shí)間t和頻率f,在時(shí)一頻單元t。內(nèi)的部分能量。地震信號(hào)的時(shí)頻分析一般將時(shí)頻分析方法分為線性和非線性兩種。典型的線性時(shí)頻表示有短時(shí)傅葉變換（簡(jiǎn)記為STFT）、Gabor展開和小波變換（WaveletTransformation，簡(jiǎn)記為WT）等。非線性時(shí)頻方法是一種二次時(shí)頻表示方法也稱為雙線性），最典型的是WVD（Wigner-Vi^eDistribution）和Cohen類。采用的地震信號(hào)如下圖所示：圖1時(shí)間域的信號(hào)圖2用FFT方法求取的能量譜1、短時(shí)傅里葉變換(STFT)傳統(tǒng)的傅立葉變換只在頻率域具備了局部分析的能力，而在時(shí)間域不具有這種能力。要得到滿足同時(shí)獲得時(shí)間和頻率的局部分析能力的要求，一種最基本的方法就是:取出信號(hào)在所關(guān)心時(shí)刻附近的一小段，而忽略信號(hào)的其它部分，對(duì)其作傅立葉變換，即可得到這一特定時(shí)刻的頻率分量。因?yàn)樗〉臅r(shí)間長度與整個(gè)信號(hào)相比很短，所以將這種方法稱為短時(shí)傅立葉變換(STFT)，它是時(shí)頻分析中最簡(jiǎn)單的形式。短時(shí)傅立葉變換(STFT)的基本思想:用窗函數(shù)來截取信號(hào)，假定信號(hào)在窗內(nèi)是平穩(wěn)的，采用傅立葉變換來分析窗內(nèi)信號(hào)，以便確定在那個(gè)時(shí)間存在的頻率，然后沿著信號(hào)移動(dòng)窗函數(shù)，得到信號(hào)頻率隨時(shí)間的變化關(guān)系，這就得到了我們所需要的時(shí)頻分布。STFT的物理意義在于，對(duì)于一定的分析時(shí)刻t,STFT(t,/)可以視為信號(hào)s(t)在該時(shí)刻的“局s部頻譜”，從而整個(gè)變換的結(jié)果也就能揭示信號(hào)頻譜的變化特性。短時(shí)Fourier變換的時(shí)頻分辨率受制于窗函數(shù)的形狀和寬度短時(shí)傅立葉變換的時(shí)間分辨率與分析窗函數(shù)的時(shí)間域?qū)挾瘸烧?，而其頻率分辨率與分析窗的頻寬成正比。從而，一個(gè)好的時(shí)間分辨率需要一個(gè)短的窗函數(shù)，而一個(gè)好的頻率分辨率需要一個(gè)長的窗函數(shù)。因此，短時(shí)傅立葉變換不能同時(shí)兼顧時(shí)間分辨率和頻率分辨率。當(dāng)選用的窗函數(shù)為Gaussian函數(shù)時(shí)，該變換為Gabor變換。交交一Hu當(dāng)bEg圖3用STFT做時(shí)頻分析2、S變換主要對(duì)S變換的定義、推導(dǎo)和特性進(jìn)行詳細(xì)闡述，S變換綜合短時(shí)傅立葉變換和小波變換的優(yōu)點(diǎn)，又避免了它們的不足:它與傅立葉變換有著直接的聯(lián)系，具有無損可逆性;與短時(shí)傅立葉變換和小波變換一樣，也是一種線性時(shí)頻表示，因此不存在交義項(xiàng)的干擾;S變換具有多種分辨率，克服了短時(shí)傅立葉變換固定分辯率的不足;S變換中含有相位因子，這是小波變換所不具備的特性?？傮w來說，S變換是近幾年發(fā)展起來的一種新的時(shí)頻分析方法.S變換結(jié)合了短時(shí)傅里葉變換和小波變換的優(yōu)點(diǎn),具有相位信息,同時(shí)該變換與小波變換一樣,其時(shí)頻窗可以調(diào)節(jié)大小以適應(yīng)非平穩(wěn)信號(hào)的特點(diǎn).S變換的這些優(yōu)點(diǎn),使得它在地球物理方面得到廣泛的應(yīng)用。S變換首先是由Stockwell等人提出的，是以Morlet小波為基本小波的連續(xù)小波變換的延展。在S變換中，簡(jiǎn)諧波與高斯函數(shù)的乘機(jī)構(gòu)成了基本小波，因?yàn)楹?jiǎn)諧波在時(shí)間域可以作伸縮變換，而高斯函數(shù)則進(jìn)行伸縮和平移。ST也可以認(rèn)為是CWT的“相位校正”。函數(shù)h(t)的S變換表示為：ST(t,f)=』+sh(t)to(t-1)e-j2兀ftdt-8其中：TOC\o"1-5"\h\z，、小 1 z t2、3(t)=—e=exp(———)o%：2兀 2o2f (T—t)2f2 -八，ST(t,f)=j+sh(t)dexp((一4)exp(—j2kft)dt—6 、/：2k 2可以用w()=exp(-壬y2)表示①(t)的傅里葉變換，其中丫和f是相同的意義。f2SPrLh=62,N1^249.5,lin.scale,imiagiesc,Threshold=5%50 100 1S0 20Ci 250 300 350 400 450Timis[e]圖4S變換時(shí)頻分析3、WVDWigner-Ville變換是1932年由Wigner首次提出的，并應(yīng)用于量子力學(xué)領(lǐng)域，后來Ville等人將其引入到信號(hào)分析處理領(lǐng)域。 1966年，Cohen發(fā)現(xiàn)各種發(fā)現(xiàn)各種時(shí)頻分析只是Wigner-Ville變換的不同形式，可以統(tǒng)一起來，成為Cohen類雙線性時(shí)頻分析。信號(hào)s(t)的Wigner-Ville變換用公式表示為w(t,f)」+6z(t+T).z*(t—T).exp(—j2兀Tf)dTz —6 2 2式中z(t)是s(t)的解析信號(hào)，即z(t)=s(t)+js(t)笆h(t)=s(t)+j/^^dN+jH[s(t)]

t—|L1—6h[s(t)]是實(shí)信號(hào)s(t)的Hilbert變換。Wigner-Ville變換也可以用解析信號(hào)的頻譜來表示(t,f)」+6z(f―-),Z*(f+-).exp(—j2kt)d—6 乙 乙(完整word版)時(shí)頻分析在地震數(shù)據(jù)處理中的應(yīng)用從上面兩種不同形式的Wigner-Ville變換表達(dá)式中可以看出，式中不包含任何窗函數(shù)，從而避免了線性時(shí)頻變換中實(shí)踐分辨率和頻率分辨率的相互牽制，難以兼顧的問題。因?yàn)閃igner-Ville變換的時(shí)間-帶寬積可達(dá)到Heisenberg測(cè)不準(zhǔn)原理(不確定原理給出的下界，故可以證明，沒有任何一種時(shí)頻變換方法的時(shí)間-頻率分辨率及聚集性能出其右。在Wigner-Ville變換表達(dá)式中，信號(hào)s(t)出現(xiàn)了兩次，故稱之為雙線性時(shí)頻變換。Wigner-Ville變換不是線性的，即兩信號(hào)之和的Wigner-Ville變換不等于每一個(gè)信號(hào)的Wigner-Ville變換之和，其中多出了一個(gè)附加項(xiàng)。令s(t)=s(t)+s(t)則有W(t,f)=W(t,f)+W(t,f)+2Re{Ws(t,f)}其中卡t TW(t,f)=]s(t+—)s*(t——)e-j2kftdtsis2 1 22 2-8式中前兩項(xiàng)是自由項(xiàng)(autoterms)，第三項(xiàng)是交叉項(xiàng)(crossterms)。交叉項(xiàng)常常導(dǎo)致時(shí)頻平面上出現(xiàn)偽影現(xiàn)象：交叉項(xiàng)是實(shí)