基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤時(shí)頻分析*

范興利,成 谷,2

(1. 中山大學(xué)地球科學(xué)與地質(zhì)工程學(xué)院，廣東 廣州 510275；2. 廣東省地質(zhì)過(guò)程與礦產(chǎn)資源探查重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510275)

對(duì)于地震信號(hào)這種典型的非平穩(wěn)信號(hào)來(lái)說(shuō)，傳統(tǒng)的傅里葉分析方法不能完整地刻畫(huà)其時(shí)變特征，需采用時(shí)間-頻率的聯(lián)合分析(時(shí)頻分析)方式，將一維的時(shí)間信號(hào)映射到二維的時(shí)頻平面，得到全面反映被觀測(cè)信號(hào)的時(shí)間-頻率聯(lián)合特征。常用的時(shí)頻分析方法有短時(shí)傅里葉變換(Short-Time Fourier Transform, 簡(jiǎn)稱(chēng)STFT)、連續(xù)小波變換(Continuous Wavelet Transform，簡(jiǎn)稱(chēng)CWT)、S變換(S Transform，簡(jiǎn)稱(chēng)ST)、廣義S變換(Generalized S Transform，簡(jiǎn)稱(chēng)GST)、匹配追蹤(Matching Pursuit，簡(jiǎn)稱(chēng)MP)方法等。

短時(shí)傅里葉變換也稱(chēng)加窗傅里葉變換，通過(guò)窗函數(shù)的加入對(duì)信號(hào)的局部特征進(jìn)行分析。但由于采用固定的窗函數(shù)，存在兩個(gè)明顯的弊端：一是當(dāng)分析包含兩個(gè)以上的多分量信號(hào)時(shí)，很難使同一個(gè)窗函數(shù)同時(shí)滿(mǎn)足幾種不同的要求；二是窗函數(shù)的形狀在移動(dòng)過(guò)程中不能隨時(shí)間發(fā)生改變。其實(shí)質(zhì)是只具有單一分辨率，時(shí)頻分辨率的高低取決于所選窗函數(shù)[1-5]。

連續(xù)小波變換繼承和發(fā)展了短時(shí)傅里葉變換的局部化思想，引入了尺度因子，克服了短時(shí)傅里葉變換單一分辨率的缺點(diǎn)，窗口大小隨頻率變化，具有更高的時(shí)頻分辨率，對(duì)信號(hào)具有自適應(yīng)性[6]。但小波變換中尺度因子與頻率沒(méi)有直接關(guān)聯(lián)起來(lái)，實(shí)際上是利用時(shí)間-尺度函數(shù)來(lái)分析非平穩(wěn)信號(hào)，而不是嚴(yán)格的時(shí)間-頻率的方式，結(jié)果不是真正意義上的時(shí)頻譜；對(duì)于中等或者高頻信號(hào)要想獲得較高的頻率分辨率，小波變換的方式并不高效。此外利用不同母函數(shù)進(jìn)行小波變換分析信號(hào)時(shí)時(shí)常會(huì)得到不同的結(jié)論和效果，甚至出現(xiàn)無(wú)解的情況[6-9]。

S變換既可以看成是一種加窗的傅里葉變換，同時(shí)也可以看成是以Morlet小波為基本小波的連續(xù)小波變換。它結(jié)合了短時(shí)傅里葉變換和小波變換各自的優(yōu)點(diǎn)，解決了短時(shí)傅里葉變換中不能改變分析窗口頻率大小的問(wèn)題，又引入了小波變換中的多分辨分析思想，同時(shí)與傅里葉譜保持著直接的聯(lián)系。S變換較小波變換引入了相位因子，可無(wú)損恢復(fù)原始信號(hào)。但在S變換中，基本小波是固定的，這使其在應(yīng)用中受到一定限制，在有些情況下不適用。廣義S變換在S變換的基礎(chǔ)上引入了兩個(gè)調(diào)節(jié)參數(shù)對(duì)尺度與頻率的關(guān)系進(jìn)行調(diào)節(jié)，可靈活地調(diào)節(jié)高斯窗函數(shù)隨頻率尺度的變化趨勢(shì)，加快或減慢時(shí)窗寬度隨信號(hào)頻率變化的速度，更好地適應(yīng)具體信號(hào)的分析和處理[10-13]。

匹配追蹤算法是由Mallat與Zhang首次提出的，本質(zhì)是對(duì)信號(hào)進(jìn)行稀疏分解[14]。與其它的時(shí)頻分析方法相比，匹配追蹤方法的時(shí)頻分辨率最高，目前在噪聲壓制、油氣檢測(cè)及薄層識(shí)別等方面都有應(yīng)用[15-17]。但匹配追蹤傳統(tǒng)的貪婪迭代算法計(jì)算效率相對(duì)較低，很多人致力于匹配追蹤方法計(jì)算效率提高方面的研究。Liu等先后提出了基于Ricker小波和Morlet小波庫(kù)的匹配追蹤算法[18-19]，將地震信號(hào)的瞬時(shí)屬性(瞬時(shí)振幅、瞬時(shí)相位、瞬時(shí)頻率)特征引入到了匹配追蹤算法當(dāng)中，提高了算法的執(zhí)行效率。張繁昌等采用基于信號(hào)瞬時(shí)包絡(luò)和頻譜分布雙重控制的方法，得到了包含相位信息的高分辨率匹配追蹤時(shí)頻分布，后來(lái)提出了一種采用雙參數(shù)動(dòng)態(tài)掃描的快速匹配追蹤算法[20-22]，進(jìn)一步提高了匹配追蹤的計(jì)算效率。邵君對(duì)基于MP的信號(hào)稀疏分解算法進(jìn)行了研究[23]。Wang在Liu提出的基于Morlet小波庫(kù)匹配追蹤分解算法基礎(chǔ)上，通過(guò)引入控制Morlet小波時(shí)寬、帶寬的尺度參數(shù)，給出了三步法匹配追蹤算法[24]，提高了算法的效率。本文探討了以Morlet小波作為時(shí)頻原子的匹配追蹤算法中尺度因子對(duì)信號(hào)與時(shí)頻原子匹配特征的控制作用，在此基礎(chǔ)上進(jìn)行了基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤時(shí)頻分析。并利用模型數(shù)據(jù)測(cè)試了算法的準(zhǔn)確性及在噪聲壓制、薄層厚度求取等方面的應(yīng)用。

1 匹配追蹤算法思想與流程

1.1 匹配追蹤算法思想

匹配追蹤算法的基本思想是：將原始信號(hào)投影到一系列時(shí)頻原子上，即把原始信號(hào)表示為這些時(shí)頻原子的線性組合，利用這些時(shí)頻原子精確地表達(dá)原始信號(hào)。

時(shí)頻原子這一抽象的概念，在具體實(shí)現(xiàn)過(guò)程中通常具化為由頻率、相位、時(shí)延等參數(shù)確定的小波形態(tài)，Ricker小波和Morlet小波就是兩種通常采用的時(shí)頻原子形式。

從數(shù)學(xué)上表示，匹配追蹤算法思想如下：設(shè)H表示希爾伯特空間，定義H中的一個(gè)時(shí)頻原子庫(kù)D，gγ∈D為時(shí)頻原子庫(kù)中的一個(gè)時(shí)頻原子，是一單位長(zhǎng)度的矢量，即數(shù)學(xué)上滿(mǎn)足

‖gγ‖=1

(1)

對(duì)信號(hào)f(f∈H)，匹配追蹤時(shí)頻分析的思想是找到時(shí)頻原子逐步最佳地表征信號(hào)的某一特征，數(shù)學(xué)上表述為將信號(hào)進(jìn)行如下分解

f=〈f,gγ0〉gγ0+R1(f)

(2)

〈f，gγ0〉表示f與gγ0的內(nèi)積，R1(f)表示用時(shí)頻原子gγ0表示信號(hào)f所產(chǎn)生的誤差。顯然，gγ0與R1(f)是正交的，所以有

‖f‖2=‖〈f，gγ0〉‖2+‖R1(f)‖2

(3)

為了使殘差能量‖R1(f)‖最小，須選擇gγ0∈D使得‖〈f,gγ0〉‖最大。匹配追蹤采用迭代算法，開(kāi)始時(shí)設(shè)R0(f)=f，進(jìn)行到第k次迭代時(shí)，上一步迭代的殘差信號(hào)為Rk-1(f)，選擇gγk∈D，使gγk與Rk-1(f)最匹配，即Rk-1(f)被分解為

Rk-1(f)=〈Rk-1(f),gγk〉gγk+Rk(f)

(4)

重復(fù)此分解過(guò)程，直到殘差能量小于所設(shè)閾值，信號(hào)最終被分解為

(5)

將信號(hào)和時(shí)頻原子抽象表示為高維空間中的矢量，則匹配追蹤的算法思想可由圖1表示。其中，帶箭頭虛線表示當(dāng)前待分解信號(hào)，細(xì)實(shí)線表示針對(duì)當(dāng)前待分解信號(hào)在時(shí)頻原子庫(kù)中找到的最佳時(shí)頻原子，即在眾多的時(shí)頻原子中待分解信號(hào)在其上投影最大的一個(gè)(虛線所示矢量向細(xì)實(shí)線所示矢量投影值最大)，粗實(shí)線表示投影殘差。上一步的投影殘差作為下一步的信號(hào)(圖1左側(cè)的垂直虛線與圖1右側(cè)的垂直粗實(shí)線表示同一矢量)重復(fù)匹配過(guò)程。通常一個(gè)信號(hào)經(jīng)有限次分解以后殘差即可足夠小，即匹配追蹤算法的實(shí)質(zhì)為稀疏信號(hào)分解。

圖1 匹配追蹤分解算法圖示

1.2 匹配追蹤算法流程

在匹配追蹤算法中，時(shí)頻原子形式和參數(shù)的確定是關(guān)鍵。對(duì)于時(shí)頻原子中的相位、頻率和時(shí)延等信息，通常通過(guò)希爾伯特變換確定。

匹配追蹤算法具體流程如下：

1)將實(shí)地震道信號(hào)通過(guò)希爾伯特變換轉(zhuǎn)換成復(fù)地震道信號(hào)；

2)尋找復(fù)地震道信號(hào)包絡(luò)最大值處，此處的時(shí)間作為時(shí)間延遲u，瞬時(shí)頻率作為頻率參數(shù)ζ，瞬時(shí)相位作為相位參數(shù)φ；

3)以利用希爾伯特變換求得的時(shí)間延遲、瞬時(shí)頻率、瞬時(shí)相位等參數(shù)作為待尋找的時(shí)頻原子參數(shù)的初值，在其周?chē)欢ㄠ徲騼?nèi)進(jìn)行擾動(dòng)構(gòu)建多個(gè)時(shí)頻原子，并將待分解信號(hào)向時(shí)頻原子進(jìn)行投影，對(duì)應(yīng)投影值最大的參數(shù)即為尋找的最佳時(shí)頻原子的參數(shù)。

4)從地震信號(hào)中減去上面生成的最佳時(shí)頻原子的實(shí)部，并將剩余值作為新的地震信號(hào)；

5)重復(fù)步驟1-4，直到迭代誤差小于閾值。

上述流程的第3步驟中，在計(jì)算機(jī)實(shí)現(xiàn)時(shí)通常需采用多重(循環(huán)的重?cái)?shù)取決于時(shí)頻原子的參數(shù)個(gè)數(shù))循環(huán)的貪婪迭代尋找算法，因此計(jì)算效率相對(duì)較低。

2 基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤

在對(duì)地震信號(hào)進(jìn)行匹配追蹤時(shí)頻分析時(shí)，主要采用Ricker和Morlet兩種小波作為時(shí)頻原子。Morlet小波在時(shí)域中的表達(dá)式為

(6)

其中，ζ是頻率參數(shù)，u是時(shí)間延遲，φ是相位參數(shù)，σ是尺度參數(shù)。

以Morlet小波作為時(shí)頻原子，其時(shí)頻原子包括四個(gè)參變量，即γ=(σ,ζ,u,φ)，σ，ζ，u，φ分別表示時(shí)頻原子中的尺度因子、頻率因子、時(shí)間延遲和相位因子。

基于Morlet小波的匹配追蹤算法將地震信號(hào)分解表示為一系列Morlet小波的線性組合。經(jīng)過(guò)N次迭代，地震信號(hào)f(t)被分解成如下形式

(7)

其中，mn為在分解過(guò)程中與信號(hào)局部特征匹配最佳的Morlet小波時(shí)頻原子，對(duì)應(yīng)圖1中綠線所示矢量，an為分解過(guò)程中信號(hào)向最佳時(shí)頻原子mn投影值(圖1中虛線向細(xì)實(shí)線方向投影值)大小，RN(f)為剩余量(圖1中粗實(shí)線所示)，當(dāng)達(dá)到閾值后可認(rèn)為是噪聲 。

在匹配追蹤時(shí)頻分析方法中，選取的時(shí)頻原子應(yīng)能很好地匹配信號(hào)的局部特征，分析時(shí)頻原子中參數(shù)對(duì)信號(hào)特征的影響有重要的意義。在以Morlet小波為時(shí)頻原子的匹配追蹤算法中，Morlet小波中的尺度參數(shù)對(duì)于待分解信號(hào)向時(shí)頻原子的投影能量具有重要的影響，因此本文對(duì)Morlet小波時(shí)頻原子中的尺度參數(shù)對(duì)信號(hào)特征的影響進(jìn)行了分析。

2.1 Morlet小波時(shí)頻原子中尺度參數(shù)的作用

對(duì)一個(gè)連續(xù)的Morlet小波來(lái)說(shuō)，尺度σ是一個(gè)重要的自適應(yīng)參數(shù)，它控制著小波在時(shí)間域的寬度和頻譜的寬度。圖2所示是頻率ζ=50，u=0.2 s，φ=0，尺度σ取不同值(1、1.5、2)時(shí)所對(duì)應(yīng)的Morlet小波在時(shí)域和頻域的展布情況，從圖中可以看出，尺度取1、1.5、2時(shí)，信號(hào)時(shí)間域延續(xù)時(shí)間逐步加大，旁瓣個(gè)數(shù)增多幅度加大，頻帶寬度逐步變窄。即尺度取值較小時(shí)，Morlet小波具有相對(duì)較窄的時(shí)寬和較寬的頻寬；當(dāng)尺度取值較大時(shí)，Morlet小波具有相對(duì)較寬的時(shí)寬和較窄的頻寬。因此，尺度參數(shù)對(duì)時(shí)頻原子在時(shí)間域的形態(tài)具有較大的影響。

本文分析了匹配追蹤算法中時(shí)頻原子的尺度參數(shù)選取的準(zhǔn)確度對(duì)投影(信號(hào)對(duì)時(shí)頻原子的投影)能量的影響。在時(shí)頻原子其它參數(shù)相同(頻率為50 Hz、時(shí)延為0.2 s、相位為π/4)的情況下、以尺度參數(shù)為2的時(shí)頻原子作為假想的信號(hào)，其它尺度參數(shù)為1到3之間以0.1為間隔的共21個(gè)尺度參數(shù)組成的Morlet小波作為時(shí)頻原子與信號(hào)之間求取投影能量。發(fā)現(xiàn)尺度參數(shù)對(duì)投影能量的大小有較大的影響(見(jiàn)圖3)。從圖中可以看出，在其它參數(shù)均相同，僅尺度參數(shù)不同的情況下，不同的尺度參數(shù)得到不同的投影能量。在真實(shí)尺度參數(shù)2處，投影能量最強(qiáng)，越偏離真實(shí)尺度參數(shù)，投影能量越弱。由此可見(jiàn)：Morlet小波中的尺度因子對(duì)時(shí)頻原子與信號(hào)的匹配特性具有重要的控制作用。由于尺度因子的引入，Morlet小波比Ricker小波更適合作為匹配追蹤時(shí)頻分析方法的時(shí)頻原子。因此本文在進(jìn)行匹配追蹤算法設(shè)計(jì)時(shí)，采用Morlet小波作為時(shí)頻原子的基本形式，并考慮到尺度因子對(duì)時(shí)頻原子與信號(hào)局部特征匹配特性的控制作用，在對(duì)時(shí)頻原子最佳參數(shù)確定時(shí)對(duì)尺度參數(shù)進(jìn)行一維優(yōu)先尋優(yōu)迭代。

圖2 尺度對(duì)Morlet小波在時(shí)頻域展布的影響

2.2 基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤算法流程

基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤算法流程如下：

1)將實(shí)地震道信號(hào)通過(guò)希爾伯特變換轉(zhuǎn)換成復(fù)地震道信號(hào)；

2)尋找復(fù)地震道信號(hào)包絡(luò)最大值處，此處的時(shí)間作為時(shí)間延遲u，瞬時(shí)頻率作為頻率參數(shù)ζ，瞬時(shí)相位作為相位參數(shù)φ；

3)在固定其它參數(shù)情況下，對(duì)尺度參數(shù)進(jìn)行-維尋優(yōu)選代求取尺度參數(shù)σ，σ通過(guò)下式求得

(8)

4)求取最佳振幅an

(9)

其對(duì)應(yīng)的尺度參數(shù)即為最佳的尺度參數(shù)。

5)在求得的瞬時(shí)頻率、瞬時(shí)相位、時(shí)間延遲、最佳尺度參數(shù)的基礎(chǔ)上進(jìn)行局部微調(diào)得到當(dāng)前步驟最佳的時(shí)頻原子。從地震信號(hào)中減去上面生成的時(shí)頻原子的實(shí)部，并將剩余值作為新的地震信號(hào)；

6)重復(fù)步驟1-5，直到迭代誤差小于閾值；

7)最后采用Wigner-Ville分布公式分別計(jì)算生成小波的時(shí)頻譜，然后進(jìn)行疊加，得到總的時(shí)頻分布

(10)

式(10)中WVDmn(t,ω)是mn(t)的Wigner-Ville分布。

圖3 時(shí)頻原子的尺度參數(shù)與投影能量的關(guān)系

3 模型測(cè)試與分析

3.1 算法準(zhǔn)確性測(cè)試

本文采用模型數(shù)據(jù)對(duì)基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤算法進(jìn)行了測(cè)試。圖4是用Ricker小波合成的一段地震記錄。圖中橫向第1道數(shù)據(jù)包含兩個(gè)頻率為10 Hz、延時(shí)分別為0.2和0.9 s、相位分別為0和π/2的Ricker小波；第2道數(shù)據(jù)包含兩個(gè)頻率為20 Hz，延時(shí)分別為0.3和0.6 s，相位分別為0和π/2的Ricker小波；第3道數(shù)據(jù)包含3個(gè)頻率為30 Hz，延時(shí)分別為0.7、1.1、1.15 s，相位分別為π/2、π/2和0的Ricker小波；第4道是上述所有小波(即前3道數(shù)據(jù))疊加后的合成信號(hào)。圖5是對(duì)圖4中的合成信號(hào)采用基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤方法得到的幾個(gè)時(shí)頻原子，圖中第1道是合成信號(hào)，從第2道到第8道是依次迭代得到的時(shí)頻原子，第9道以后是7次迭代后的殘差能量。從圖中可以看出，組成合成信號(hào)的不同頻率、延時(shí)和相位的ricker子波(表示信號(hào)中的局部特征)能很好地被Morlet子波時(shí)頻原子匹配，充分體現(xiàn)了匹配追蹤算法的稀疏信號(hào)分解本質(zhì)。對(duì)圖4中的合成信號(hào)分別進(jìn)行短時(shí)傅里葉變換(STFT)、連續(xù)小波變換(CWT)、廣義S變換(GST)、匹配追蹤分解(MP)，得到四種時(shí)頻分布，見(jiàn)圖6所示。觀察圖中得到的四種時(shí)頻分布不難發(fā)現(xiàn)，短時(shí)傅里葉變換時(shí)頻分布效果較差，時(shí)間、頻率分辨率都不佳，很難區(qū)分出各個(gè)小波，分別出現(xiàn)在0.6和1.1 s附近的兩個(gè)小波不能被區(qū)分開(kāi)，出現(xiàn)了“同時(shí)不同頻率”的假象。相對(duì)短時(shí)傅里葉變換，小波變換和廣義S變換得到的時(shí)頻分布效果要好很多，時(shí)間和頻率分辨率都較高，體現(xiàn)了多分辨率分析的優(yōu)勢(shì)，基本能夠看出各個(gè)小波出現(xiàn)的時(shí)間以及頻率。時(shí)頻分布效果最好的是匹配追蹤得到的結(jié)果，從圖中可以看出，在0.2，0.6以及1.1 s附近，各自出現(xiàn)了兩個(gè)相隔很近的小波，在小波變換和廣義S變換圖中，均出現(xiàn)了能量團(tuán)“藕斷絲連”的現(xiàn)象，而在匹配追蹤時(shí)頻圖中可以看到這三個(gè)時(shí)刻出現(xiàn)的相隔較近小波都能夠被很好的區(qū)分開(kāi)，能量分布更為集中。由此可以看出利用匹配追蹤算法進(jìn)行時(shí)頻分析的效果要優(yōu)于短時(shí)傅里葉變換，連續(xù)小波變換以及廣義S變換等傳統(tǒng)時(shí)頻分析方法。

圖4 合成信號(hào)

圖5 匹配追蹤得到的時(shí)頻原子

3.2 算法應(yīng)用測(cè)試

采用基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤時(shí)頻分析方法，本文針對(duì)模型對(duì)算法在去噪和薄層厚度求取等方面的應(yīng)用進(jìn)行了測(cè)試。

圖7是利用基于Morlet小波尺度參數(shù)尋優(yōu)的匹配追蹤時(shí)頻分析方法進(jìn)行的去噪方面的測(cè)試。其中圖7(a)圖為不含噪音的合成地震記錄，圖7(b)圖為加上高斯白噪聲的含噪剖面。對(duì)含噪剖面中的每一道進(jìn)行匹配追蹤時(shí)頻分析，給定誤差閥值，得出對(duì)應(yīng)每一道信號(hào)的重建信號(hào)和誤差信號(hào)。將各道的重建信號(hào)排列在一起形成圖7(c)，即為利用匹配追蹤方法去除噪聲的剖面，將各道的誤差道排列在一起形成圖7(d)，即為去除的噪聲剖面。從圖中可以看出，去噪效果較好。

圖6 四種時(shí)頻分析方法比較

圖7 匹配追蹤方法去除噪聲

圖8(a)是一楔形反射系數(shù)模型，一層砂巖楔形侵入泥巖空間，砂巖速度取為3 700 m/s，砂巖最大厚度100 m，一共有30道，采樣間隔2 ms。利用反射系數(shù)模型生成合成地震記錄。對(duì)合成地震剖面利用匹配追蹤分解得到三維數(shù)據(jù)體，提取40 Hz的單頻剖面如圖8(b)所示。根據(jù)此楔形體模型薄層調(diào)諧厚度與振幅能量之間的關(guān)系，當(dāng)薄層厚度為地震子波波長(zhǎng)的1/4時(shí)，其振幅值最大，能量最強(qiáng)，這一厚度也稱(chēng)為調(diào)諧厚度。所以，利用單頻剖面中的能量極值可估算對(duì)應(yīng)該頻率的薄層調(diào)諧厚度。對(duì)于40 Hz單頻剖面而言，地震子波波長(zhǎng)的1/4(即調(diào)諧厚度)為

(11)

式中，v是薄層砂巖體速度，f是單頻剖面對(duì)應(yīng)的頻率，h是根據(jù)調(diào)諧理論計(jì)算出的對(duì)應(yīng)該頻率的調(diào)諧厚度。觀察40 Hz單頻剖面圖，可以看出最強(qiáng)能量振幅出現(xiàn)在第7道，根據(jù)模型設(shè)置參數(shù)，30道對(duì)應(yīng)厚度100 m，則第7道對(duì)應(yīng)的厚度為

(12)

式中，h′是根據(jù)單頻剖面計(jì)算出的振幅能量最強(qiáng)位置的薄層厚度，n是強(qiáng)能量振幅處道數(shù)，H是楔形體薄層最大厚度，N是總的地震道數(shù)。根據(jù)調(diào)諧頻率計(jì)算出的調(diào)諧厚度與真實(shí)薄層厚度之間的誤差量Δh為

Δh=|h-h′|=|23.12-23.33|=0.20 m

(13)

楔形體砂巖按道數(shù)的最小遞變厚度Δh′是

(14)

對(duì)比Δh與Δh′，可以看出Δh?Δh′，因此，估算誤差是可以接受的，利用40 Hz單頻剖面能夠達(dá)到對(duì)該楔形體薄層厚度進(jìn)行求取的目的。

圖8 匹配追蹤方法求取薄層厚度

4 結(jié) 語(yǔ)

與短時(shí)傅里葉變換、連續(xù)小波變換、廣義S變換等傳統(tǒng)的時(shí)頻分析方法比較，利用匹配追蹤分解算法得到的時(shí)頻分布具有更高的時(shí)頻分辨率。匹配追蹤時(shí)頻分析方法本質(zhì)上是對(duì)信號(hào)進(jìn)行稀疏分解，將信號(hào)分解為多個(gè)時(shí)頻原子的疊加。在算法實(shí)現(xiàn)過(guò)程中，通常采用Ricker子波和Morlet子波作為時(shí)頻原子的基本形式建立時(shí)頻原子庫(kù)，然后在時(shí)頻原子庫(kù)中尋找與當(dāng)前待分解信號(hào)局部特征匹配最佳的時(shí)頻原子。時(shí)頻原子庫(kù)的建立是匹配追蹤算法中關(guān)鍵的步驟，也是計(jì)算量消耗最多的步驟。傳統(tǒng)的匹配追蹤貪婪迭代算法在實(shí)現(xiàn)過(guò)程中需針對(duì)頻率、相位、時(shí)延、振幅等多個(gè)時(shí)頻原子參數(shù)采用多重循環(huán)迭代方式尋找最優(yōu)的時(shí)頻原子參數(shù)，算法效率較低。

本文對(duì)Morlet小波時(shí)頻原子中的尺度因子對(duì)信號(hào)的時(shí)域形態(tài)和頻域范圍的影響進(jìn)行了分析，分析結(jié)果表明：尺度因子取值較小時(shí)，Morlet小波具有相對(duì)較窄的時(shí)寬和較寬的頻寬，尺度因子取值較大時(shí)，Morlet小波具有較寬的時(shí)寬和較窄的頻寬。因此Morlet小波時(shí)頻原子中的尺度因子對(duì)時(shí)頻原子的時(shí)域形態(tài)具有重要的影響，進(jìn)而將影響時(shí)頻原子與信號(hào)局部特征的匹配特性(信號(hào)向時(shí)頻原子投影的投影值大小表征時(shí)頻原子與信號(hào)局部特征匹配特性的好壞)。在設(shè)定時(shí)頻原子中其它參數(shù)(頻率、相位和時(shí)延)相同，利用某一尺度因子生成信號(hào)，利用不同的尺度因子生成時(shí)頻原子的情況下，本文計(jì)算分析了信號(hào)與不同尺度因子生成的時(shí)頻原子的投影值大小，分析結(jié)果表明：在其它參數(shù)相同僅尺度因子不同的情況下，信號(hào)向時(shí)頻原子投影的投影值不同。在真實(shí)尺度因子處，投影值最大。越偏離真實(shí)尺度因子，投影值越小。進(jìn)而驗(yàn)證了尺度因子對(duì)時(shí)頻原子和信號(hào)局部特征匹配特性的控制作用。

鑒于以上結(jié)論，本文在以Morlet小波作為時(shí)頻原子進(jìn)行匹配追蹤算法實(shí)現(xiàn)時(shí)，首先利用希爾伯特變換獲得時(shí)頻原子頻率、相位、時(shí)延等參數(shù)的初始值，利用這些初始參數(shù)值對(duì)尺度因子進(jìn)行一維尋優(yōu)。在尋找到最優(yōu)的尺度因子后以其作為尺度因子的初值與其它參數(shù)一起進(jìn)行局部微調(diào)，提高了計(jì)算效率。通過(guò)模型試算形象說(shuō)明了匹配追蹤算法的稀疏信號(hào)分解的本質(zhì)，對(duì)算法在去噪和薄層厚度求取等方面的應(yīng)用進(jìn)行了模型測(cè)試，驗(yàn)證了算法的準(zhǔn)確性。

參考文獻(xiàn)：

[1]ALLEN J B，RABINER L R. A unified approach to short-time Fourier analysis and synthesis [J]. Proceedings of the IEEE, 1977, 65(11): 1558-1564.

[2]GABOR D. Theory of communication [J]. Journal of the IEEE, 1946, 93: 429-441.

[3]KOENIG W, DUNN H K, LACY L Y. The sound spectrograph [J]. The Journal of the Acoustical Society of America, 1946, 18(1): 19-49.

[4]邊海龍. 非平穩(wěn)信號(hào)聯(lián)合時(shí)頻分析方法的若干問(wèn)題研究與應(yīng)用[D]. 成都：電子科技大學(xué), 2008.

[5]董建華，顧漢明，張星. 幾種時(shí)頻分析方法的比較及應(yīng)用[J]. 工程地球物理學(xué)報(bào), 2007(4): 312-316.

[6]吳勇. 基于小波的信號(hào)去噪方法研究[D]. 武漢：武漢理工大學(xué), 2007.

[7]劉麗娟. 時(shí)頻分析技術(shù)及其應(yīng)用[D]. 成都：成都理工大學(xué), 2008.

[8]吳偉龍. 基于小波變換的地震信號(hào)瞬時(shí)參數(shù)提取方法研究[D]. 大慶：東北石油大學(xué), 2011.

[9]張賢達(dá)，保錚. 非平穩(wěn)信號(hào)分析與處理[M]. 北京: 國(guó)防工業(yè)出版社, 1998: 446.

[10]姜鐳. 基于時(shí)頻譜特征的薄互層分析[D]. 成都理工大學(xué), 2009.

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

[12]楊陽(yáng). 廣義S變換時(shí)頻分析的應(yīng)用研究[D]. 哈爾濱：哈爾濱工程大學(xué)，2011.

[13]STOCKWELL R G, MANSINHA L, LOWE R P. Localization of the complex spectrum: the S transform [J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.

[14]MALLAT S G, ZHANG Z F. Matching pursuits with time-frequency dictionaries [J]. IEEE Transactions on Signal Processing, 1993, 41(12): 3397-3415.

[15]CASTAGNA J P, SUN S, SIEGFRIED R W. Instantaneous spectral analysis: detection of low-frequency shadows associated with hydrocarbons [J]. The Leading Edge, 2003, 22(2): 120-127.

[16]孫萬(wàn)元，張會(huì)星，杜藝可. 匹配追蹤時(shí)頻分析及其在油氣檢測(cè)中的應(yīng)用[J]. 山東科技大學(xué)學(xué)報(bào)：自然科學(xué)版,2011,30(4): 51-57.

[17]PARTYKA G, GRIDLEY J, LOPEZ J. Interpretational applications of spectral decomposition in reservoir characterization [J]. The Leading Edge, 1999, 18(3): 353.

[18]LIU J. Time-frequency decomposition based on Ricker wavelet [J]. SEG Technical Program Expanded Abstracts, 2004, 23(1): 1937.

[19]LIU J. Matching pursuit decomposition using Morlet wavelets [J]. SEG Technical Program Expanded Abstracts, 2005, 24(1): 786.

[20]張繁昌，李傳輝. 非平穩(wěn)地震信號(hào)匹配追蹤時(shí)頻分析[J]. 物探與化探, 2011, 35(4): 546-552.

[21]張繁昌，李傳輝，印興耀. 基于動(dòng)態(tài)匹配小波庫(kù)的地震數(shù)據(jù)快速匹配追蹤[J]. 石油地球物理勘探, 2010, 45(5): 667-673.

[22]張繁昌，李傳輝. 基于正交時(shí)頻原子的地震信號(hào)快速匹配追蹤[J]. 地球物理學(xué)報(bào), 2012, 55(1): 277-283.

[23]邵君. 基于MP的信號(hào)稀疏分解算法研究[D]. 成都：西南交通大學(xué), 2006.

[24]WANG Y. Seismic time-frequency spectral decomposition by matching pursuit [J]. Geophysics, 2007, 72(1): V13-V20.