基于Morlet小波尺度參數尋優的匹配追蹤時頻分析*

范興利,成 谷,2

(1. 中山大學地球科學與地質工程學院，廣東 廣州 510275；2. 廣東省地質過程與礦產資源探查重點實驗室，廣東 廣州 510275)

對于地震信號這種典型的非平穩信號來說，傳統的傅里葉分析方法不能完整地刻畫其時變特征，需采用時間-頻率的聯合分析(時頻分析)方式，將一維的時間信號映射到二維的時頻平面，得到全面反映被觀測信號的時間-頻率聯合特征。常用的時頻分析方法有短時傅里葉變換(Short-Time Fourier Transform, 簡稱STFT)、連續小波變換(Continuous Wavelet Transform，簡稱CWT)、S變換(S Transform，簡稱ST)、廣義S變換(Generalized S Transform，簡稱GST)、匹配追蹤(Matching Pursuit，簡稱MP)方法等。

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

連續小波變換繼承和發展了短時傅里葉變換的局部化思想，引入了尺度因子，克服了短時傅里葉變換單一分辨率的缺點，窗口大小隨頻率變化，具有更高的時頻分辨率，對信號具有自適應性[6]。但小波變換中尺度因子與頻率沒有直接關聯起來，實際上是利用時間-尺度函數來分析非平穩信號，而不是嚴格的時間-頻率的方式，結果不是真正意義上的時頻譜；對于中等或者高頻信號要想獲得較高的頻率分辨率，小波變換的方式并不高效。此外利用不同母函數進行小波變換分析信號時時常會得到不同的結論和效果，甚至出現無解的情況[6-9]。

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

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

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

1.1 匹配追蹤算法思想

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

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

從數學上表示，匹配追蹤算法思想如下：設H表示希爾伯特空間，定義H中的一個時頻原子庫D，gγ∈D為時頻原子庫中的一個時頻原子，是一單位長度的矢量，即數學上滿足

‖gγ‖=1

(1)

對信號f(f∈H)，匹配追蹤時頻分析的思想是找到時頻原子逐步最佳地表征信號的某一特征，數學上表述為將信號進行如下分解

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

(2)

〈f，gγ0〉表示f與gγ0的內積，R1(f)表示用時頻原子gγ0表示信號f所產生的誤差。顯然，gγ0與R1(f)是正交的，所以有

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

(3)

為了使殘差能量‖R1(f)‖最小，須選擇gγ0∈D使得‖〈f,gγ0〉‖最大。匹配追蹤采用迭代算法，開始時設R0(f)=f，進行到第k次迭代時，上一步迭代的殘差信號為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)

重復此分解過程，直到殘差能量小于所設閾值，信號最終被分解為

(5)

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

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

1.2 匹配追蹤算法流程

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

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

1)將實地震道信號通過希爾伯特變換轉換成復地震道信號；

2)尋找復地震道信號包絡最大值處，此處的時間作為時間延遲u，瞬時頻率作為頻率參數ζ，瞬時相位作為相位參數φ；

3)以利用希爾伯特變換求得的時間延遲、瞬時頻率、瞬時相位等參數作為待尋找的時頻原子參數的初值，在其周圍一定鄰域內進行擾動構建多個時頻原子，并將待分解信號向時頻原子進行投影，對應投影值最大的參數即為尋找的最佳時頻原子的參數。

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

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

上述流程的第3步驟中，在計算機實現時通常需采用多重(循環的重數取決于時頻原子的參數個數)循環的貪婪迭代尋找算法，因此計算效率相對較低。

2 基于Morlet小波尺度參數尋優的匹配追蹤

在對地震信號進行匹配追蹤時頻分析時，主要采用Ricker和Morlet兩種小波作為時頻原子。Morlet小波在時域中的表達式為

(6)

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

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

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

(7)

其中，mn為在分解過程中與信號局部特征匹配最佳的Morlet小波時頻原子，對應圖1中綠線所示矢量，an為分解過程中信號向最佳時頻原子mn投影值(圖1中虛線向細實線方向投影值)大小，RN(f)為剩余量(圖1中粗實線所示)，當達到閾值后可認為是噪聲 。

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

2.1 Morlet小波時頻原子中尺度參數的作用

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

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

圖2 尺度對Morlet小波在時頻域展布的影響

2.2 基于Morlet小波尺度參數尋優的匹配追蹤算法流程

基于Morlet小波尺度參數尋優的匹配追蹤算法流程如下：

1)將實地震道信號通過希爾伯特變換轉換成復地震道信號；

2)尋找復地震道信號包絡最大值處，此處的時間作為時間延遲u，瞬時頻率作為頻率參數ζ，瞬時相位作為相位參數φ；

3)在固定其它參數情況下，對尺度參數進行-維尋優選代求取尺度參數σ，σ通過下式求得

(8)

4)求取最佳振幅an

(9)

其對應的尺度參數即為最佳的尺度參數。

5)在求得的瞬時頻率、瞬時相位、時間延遲、最佳尺度參數的基礎上進行局部微調得到當前步驟最佳的時頻原子。從地震信號中減去上面生成的時頻原子的實部，并將剩余值作為新的地震信號；

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

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

(10)

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

圖3 時頻原子的尺度參數與投影能量的關系

3 模型測試與分析

3.1 算法準確性測試

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

圖4 合成信號

圖5 匹配追蹤得到的時頻原子

3.2 算法應用測試

采用基于Morlet小波尺度參數尋優的匹配追蹤時頻分析方法，本文針對模型對算法在去噪和薄層厚度求取等方面的應用進行了測試。

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

圖6 四種時頻分析方法比較

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

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

(11)

式中，v是薄層砂巖體速度，f是單頻剖面對應的頻率，h是根據調諧理論計算出的對應該頻率的調諧厚度。觀察40 Hz單頻剖面圖，可以看出最強能量振幅出現在第7道，根據模型設置參數，30道對應厚度100 m，則第7道對應的厚度為

(12)

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

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

(13)

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

(14)

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

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

4 結 語

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

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

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

參考文獻：

[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]邊海龍. 非平穩信號聯合時頻分析方法的若干問題研究與應用[D]. 成都：電子科技大學, 2008.

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

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

[7]劉麗娟. 時頻分析技術及其應用[D]. 成都：成都理工大學, 2008.

[8]吳偉龍. 基于小波變換的地震信號瞬時參數提取方法研究[D]. 大慶：東北石油大學, 2011.

[9]張賢達，保錚. 非平穩信號分析與處理[M]. 北京: 國防工業出版社, 1998: 446.

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

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

[12]楊陽. 廣義S變換時頻分析的應用研究[D]. 哈爾濱：哈爾濱工程大學，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]孫萬元，張會星，杜藝可. 匹配追蹤時頻分析及其在油氣檢測中的應用[J]. 山東科技大學學報：自然科學版,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]張繁昌，李傳輝. 非平穩地震信號匹配追蹤時頻分析[J]. 物探與化探, 2011, 35(4): 546-552.

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

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

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

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