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

1、時(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和WignerVille分布和S變換的特點(diǎn)，并利用Matlab對(duì)一個(gè)地震信號(hào)進(jìn)行時(shí)頻分析的實(shí)現(xiàn)。關(guān)鍵詞：時(shí)頻分析；短時(shí)傅立葉變換（STFT）；S變換；WignerVille分布引言在傳統(tǒng)的信號(hào)處理領(lǐng)域，基于 Fourier 變換的信號(hào)頻域表示及其能量的頻域分布揭示了信號(hào)在頻域的特征，它們?cè)趥鹘y(tǒng)的信號(hào)分析與處理的發(fā)展史上發(fā)揮了極其重要的作用

2、。但是，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í)間的變化特性。這使我們不但能夠同

3、時(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(t，f)，使p(t，f)=在時(shí)間t和頻率f的強(qiáng)度

4、，或者P(t，f)在時(shí)間t和頻率f，在時(shí)一頻單元內(nèi)的部分能量。地震信號(hào)的時(shí)頻分析一般將時(shí)頻分析方法分為線性和非線性兩種。典型的線性時(shí)頻表示有短時(shí)傅葉變換(簡(jiǎn)記為STFT)、Gabor展開和小波變換(Wavelet Transformation，簡(jiǎn)記為WT)等。非線性時(shí)頻方法是一種二次時(shí)頻表示方法(也稱為雙線性)，最典型的是WVD(Wigner-Ville Distribution)和Cohen類。 采用的地震信號(hào)如下圖所示：圖1 時(shí)間域的信號(hào)圖2 用FFT方法求取的能量譜1、短時(shí)傅里葉變換（STFT）傳統(tǒng)的傅立葉變換只在頻率域具備了局部分析的能力，而在時(shí)間域不具有這種能力。要得到滿足同時(shí)獲得時(shí)

5、間和頻率的局部分析能力的要求，一種最基本的方法就是:取出信號(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, 可以視為信號(hào)s(t)在該時(shí)刻的“局部頻譜”，從而整個(gè)變

6、換的結(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變換。圖3 用STFT做時(shí)頻分析2、S變換主要對(duì)S變換的定義、推導(dǎo)和特性進(jìn)行詳細(xì)闡述，S變換綜合短時(shí)傅立葉變換和小波變換的優(yōu)點(diǎn)，又避免了它們的不足:它與傅立葉變換有著直接的聯(lián)系，具有無損可逆性;與短時(shí)傅立葉變換和

7、小波變換一樣，也是一種線性時(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ù)則

8、進(jìn)行伸縮和平移。ST也可以認(rèn)為是CWT的“相位校正”。函數(shù)的S變換表示為：其中：可以用表示的傅里葉變換，其中和是相同的意義。圖4 S變換時(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)的Wigner-Ville變換用公式表示為式中是的解析信號(hào)，即是實(shí)信號(hào)的Hilbert變換。Wigner-Ville變換也可以用解析信號(hào)的頻譜來表示從上面兩種不同形式的Wi

9、gner-Ville變換表達(dá)式中可以看出，式中不包含任何窗函數(shù)，從而避免了線性時(shí)頻變換中實(shí)踐分辨率和頻率分辨率的相互牽制，難以兼顧的問題。因?yàn)閃igner-Ville變換的時(shí)間-帶寬積可達(dá)到Heisenberg測(cè)不準(zhǔn)原理(不確定原理給出的下界，故可以證明，沒有任何一種時(shí)頻變換方法的時(shí)間-頻率分辨率及聚集性能出其右。在Wigner-Ville變換表達(dá)式中，信號(hào)出現(xiàn)了兩次，故稱之為雙線性時(shí)頻變換。Wigner-Ville變換不是線性的，即兩信號(hào)之和的Wigner-Ville變換不等于每一個(gè)信號(hào)的Wigner-Ville變換之和，其中多出了一個(gè)附加項(xiàng)。令則有其中式中前兩項(xiàng)是自由項(xiàng)(auto term

10、s)，第三項(xiàng)是交叉項(xiàng)(cross terms)。交叉項(xiàng)常常導(dǎo)致時(shí)頻平面上出現(xiàn)偽影現(xiàn)象：交叉項(xiàng)是實(shí)的，混雜于自由項(xiàng)成分當(dāng)中，并且其幅度是自由項(xiàng)成分的兩倍；交叉項(xiàng)是震蕩型的，每?jī)蓚€(gè)信號(hào)分量就會(huì)產(chǎn)生一個(gè)交叉項(xiàng)。交叉項(xiàng)的存在嚴(yán)重地干擾著對(duì)Wigner-Ville分布的解釋，當(dāng)信號(hào)變得非常復(fù)雜時(shí)，Wigner-Ville分布甚至變得毫無意義。圖5 WVD時(shí)頻分析4、小波變換 在短時(shí)傅立葉變換和Gabor展開中我們都使用了固定的時(shí)間窗函數(shù)，這就引出了時(shí)間分辨率和頻率分辨率的概念，時(shí)間分辨率和頻率分辨率是一對(duì)矛盾。根據(jù)海森堡的測(cè)不準(zhǔn)原理，即時(shí)間窗函數(shù)的長度越長，頻率分辨率就越高，而對(duì)于時(shí)間分辨率則越差。為了

11、平衡時(shí)間分辨率和頻率分辨率這個(gè)矛盾，可以采取對(duì)存在高頻分量的部分采用高的時(shí)間分辨率和低的頻率分辨率，而對(duì)于低頻分量則采用高的頻率分辨率和低的時(shí)間分辨率的方法，這就是多分辨分析的思想。 小波變換是一種在時(shí)間-尺度平面內(nèi)，利用多分辨率分析思想分析非平穩(wěn)信號(hào)的方法。所謂小波，就是一個(gè)滿足容許條件的一個(gè)函數(shù)族 可以看出函數(shù)族是由窗函數(shù)在時(shí)間上平移b，在尺度上伸縮a，再乘上歸一化因子后的結(jié)果，所以非平穩(wěn)信號(hào)的連續(xù)小波變換定義為 :其中是小波基函數(shù)的共軛。 將小波變換和短時(shí)傅立葉變換兩者的基函數(shù)相比較，可以看出，小波變換基函數(shù)的尺度參數(shù)決定了小波變換的多分辨分析特性，即利用時(shí)間-尺度聯(lián)合函數(shù)來分析非平穩(wěn)信

12、號(hào)的“變焦距”法，以達(dá)到分析信號(hào)局部特性的目的。 小波變換由于其本身分辨力的優(yōu)良性能，因此一經(jīng)提出，很快就成了非平穩(wěn)信號(hào)分析和處理的一大熱點(diǎn)，經(jīng)過近 20 年的發(fā)展，小波變換取得了突破性的發(fā)展，形成了多分辨分析，框架和濾波器組三大完整豐富的小波變換理論體系?，F(xiàn)在小波變換己經(jīng)被廣泛地應(yīng)用在信號(hào)的奇異性檢測(cè)、計(jì)算機(jī)視覺、圖像處理、語音分析與合成等等諸多領(lǐng)域、在分形和混沌理論中也有了很多的應(yīng)用。 聲音信號(hào)分析：1. 女聲版江南style：圖6 女聲版聲音信號(hào)圖7 FFT所得振幅信息圖8 短時(shí)傅里葉變換時(shí)頻分析圖9 S變換時(shí)頻分析 圖10 WVD時(shí)頻分析圖11 小波變換時(shí)頻分析圖12 Hilbert變

13、換后小波變換時(shí)頻分析圖13 Hilbert變換后WVD時(shí)頻分析2. 男聲版江南style：圖12 男聲版聲音信號(hào)圖13 FFT所得振幅信息圖14 短時(shí)傅里葉變換時(shí)頻分析圖15 S變換時(shí)頻分析圖16 WVD時(shí)頻分析圖17 小波變換時(shí)頻分析圖18 加入Hilbert變換小波變換時(shí)頻分析圖19 加入Hilbert變換WVD時(shí)頻分析對(duì)比男生與女生版兩段相似的聲音信號(hào)發(fā)現(xiàn)女生的頻率段基本在0.35kHz左右，而男生版的聲音信號(hào)則有一段明顯的低頻段。再者由于聲音信號(hào)內(nèi)存在的其他干擾，可能是的信號(hào)的頻率與能量等有差異，但是最基本的頻率信息還是可以看出來的。另外，交叉項(xiàng)是震蕩型的，每?jī)蓚€(gè)信號(hào)分量就會(huì)產(chǎn)生一個(gè)交

14、叉項(xiàng)。交叉項(xiàng)的存在嚴(yán)重地干擾著對(duì)Wigner-Ville分布的解釋，當(dāng)信號(hào)變得非常復(fù)雜時(shí)，Wigner-Ville分布甚至變得毫無意義。就像上面分析的聲音信號(hào)，除了看出簡(jiǎn)單的對(duì)稱性之外，其他的分析意義基本沒有。幾種不同的Wigner-Ville變換原始信號(hào)： x1=cos(2*pi*200*t).*(t>=0 & t<=0.8); x2=cos(2*pi*50*t).*(t>=0.2 & t<=1); x=x1+x2;圖20 原始信號(hào)圖21 三種WVD時(shí)頻分析圖22 WVD時(shí)頻分析三維顯示結(jié)論由于地下介質(zhì)的復(fù)雜性，地震信號(hào)是一種非平穩(wěn)隨機(jī)信號(hào)，時(shí)頻分析是

15、分析處理地震信號(hào)的一種非常有效的方法。前人在此領(lǐng)域已做了大量工作，將很多時(shí)頻分析的方法應(yīng)用在地震數(shù)據(jù)處理解釋中。本文對(duì)各種時(shí)頻分析方法的比較研究和分析其各自特點(diǎn)。并將時(shí)頻分析的方法用于地震信號(hào)分析處理中，結(jié)果表明:通過時(shí)頻分析方法對(duì)地震剖面進(jìn)行處理，可以得到瞬時(shí)地震屬性參數(shù)(包括頻率類瞬時(shí)屬性參數(shù)、瞬時(shí)振幅譜)，這些新的地震剖面對(duì)地質(zhì)解釋有一定的輔助作用。另外后面聲音信號(hào)的分析說明，時(shí)頻分析對(duì)于頻率有很好的敏感性，對(duì)于區(qū)分不同頻帶的信息有著非常好的分辨率。像不同的樂器，相同物體不同的激發(fā)聲音方式等，都可以通過時(shí)頻分析進(jìn)行解釋。時(shí)頻分析以聯(lián)合時(shí)頻分布的形式來表示信號(hào)的特性，具有很多的優(yōu)點(diǎn): 1）

16、它克服了傅里葉分析時(shí)域和頻域完全分離的缺陷，將時(shí)頻兩域聯(lián)合起來對(duì)信號(hào)進(jìn)行分析，能同時(shí)考慮到兩個(gè)方面的性能。 2）彌補(bǔ)了信號(hào)的時(shí)間能量密度和頻譜能量密度不能充分描述信號(hào)的物理特性的缺陷。 3）在時(shí)頻相平面上，可以精確地定位在某一時(shí)刻出現(xiàn)了哪些頻率分量，以及某一分量出現(xiàn)在哪些時(shí)刻。 4）對(duì)于不同情況的信號(hào)，通過對(duì)核函數(shù)施加一些約束條件，就可以設(shè)計(jì)符合期望性能的時(shí)頻分布，來滿足處理信號(hào)的要求。 5）為非平穩(wěn)信號(hào)的分析提供了有效的工具，為信號(hào)的分析開辟了新的途徑。時(shí)頻分析雖然具有很多優(yōu)點(diǎn)，但同時(shí)也具有不少缺點(diǎn): 1）由于雙線性形式的時(shí)頻分布是非線性的，使得兩個(gè)信號(hào)和的時(shí)頻分布已不再是兩個(gè)信號(hào)各自分布的

17、和，即存在交叉項(xiàng)。 2）在時(shí)頻域進(jìn)行去噪時(shí)計(jì)算復(fù)雜，不易實(shí)現(xiàn)。 3）時(shí)頻分布雖然反映了信號(hào)的能量分布，但不能用信號(hào)的“瞬時(shí)能量”來解釋某一時(shí)刻或某一頻率處的時(shí)頻分布。 現(xiàn)實(shí)中很多信號(hào)，比如語音信號(hào)，都是時(shí)變非平穩(wěn)的，時(shí)變非平穩(wěn)特性是現(xiàn)實(shí)信號(hào)的普遍規(guī)律，聯(lián)合時(shí)頻分析技術(shù)正是應(yīng)現(xiàn)實(shí)的科學(xué)和工程應(yīng)用需求而產(chǎn)生和發(fā)展起來的。對(duì)于許多信號(hào)，僅用時(shí)域或頻域里的各種方法去分析往往不能揭示信號(hào)內(nèi)部的局部特征和信息，而時(shí)頻分析作為一種能將頻譜隨時(shí)間的演變關(guān)系明確表現(xiàn)出來的新手段，自然更符合實(shí)際應(yīng)用的需要。參考文獻(xiàn)1. 鄒紅星,周小波,李衍達(dá). 時(shí)頻分析:回溯與前瞻J. 電子學(xué)報(bào),2000,09:78-84.2.

18、 劉葵,劉招君,朱建偉,張峰. 時(shí)頻分析在石油地球物理勘探中的應(yīng)用J. 世界地質(zhì),2000,03:282-285.3. 劉麗娟. 時(shí)頻分析技術(shù)及其應(yīng)用D.成都理工大學(xué),2008.4. 陳斌. 時(shí)頻分析及在地震信號(hào)分析中的應(yīng)用研究D.成都理工大學(xué),2007. 5. 于振江. 一種基于Matlab的語音信號(hào)采集與分析系統(tǒng)設(shè)計(jì)J. 科技情報(bào)開發(fā)與經(jīng)濟(jì),2012,12:112-114.6. 陳家焱,陳冬嬌,張達(dá)響. 基于Matlab的聲音信號(hào)采集與分析處理J. 計(jì)算機(jī)與現(xiàn)代化,2005,06:91-92+96.7. 題原,張勁松. 基于MATLAB的語音信號(hào)采集和分析系統(tǒng)的可視化設(shè)計(jì)J. 齊齊哈爾大學(xué)

19、學(xué)報(bào),2006,06:43-46.8. 蔣濟(jì)同,唐世振. 基于MATLAB的振動(dòng)信號(hào)采集與分析系統(tǒng)的研究A. 中國力學(xué)學(xué)會(huì)、中國振動(dòng)工程學(xué)會(huì)、中國航空學(xué)會(huì)、中國機(jī)械工程學(xué)會(huì)、中國宇航學(xué)會(huì).第九屆全國振動(dòng)理論及應(yīng)用學(xué)術(shù)會(huì)議論文集C.中國力學(xué)學(xué)會(huì)、中國振動(dòng)工程學(xué)會(huì)、中國航空學(xué)會(huì)、中國機(jī)械工程學(xué)會(huì)、中國宇航學(xué)會(huì):,2007:5.9. 宋楊潔. 基于LabVIEW與MATLAB的語言信號(hào)的采集與分析D.武漢理工大學(xué),2012.10. 周淵,王炳和,劉斌勝. 基于MATLAB的噪聲信號(hào)采集與分析系統(tǒng)的設(shè)計(jì)J. 電聲技術(shù),2004,07:52-54.附：1、時(shí)頻分析基于Matlab的實(shí)現(xiàn)clc;clear

20、;nfft=512;fid=fopen('xinhao1.txt','r');fclose(fid);%x=u(:,200);figure(1);tfrstft(u);figure(2);tfrsp(u);figure(3);tfrwv(u);2、 不同的Wigner-Ville變換clear all;N=2048;fs=2048;n=0:N;t=n*1/2000;x1=cos(2*pi*200*t).*(t>=0 & t<=0.8);x2=cos(2*pi*50*t).*(t>=0.2 & t<=1);x=x1+x2;figure(1);plot(t,x);title('原始信號(hào)');xlabel('時(shí)間 t');ylabel('幅值 A');grid on;figure(2);subplot(311);tfr,t,f=tfrwv(x',1:N,N);cont