印興耀, 許 璐, 宗兆云, 李 坤

(1.中國石油大學(華東)地球科學與技術學院,山東青島 266580; 2.海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島 266071)

匹配追蹤算法可自適應地將信號在過完備原子庫中分解為最佳匹配時頻原子的線性組合。Mallat等[1]于1993年首次提出了基于Gabor型時頻原子庫的匹配追蹤算法(Matching Pursuit),Chakraborty和Okaya[2]于1995年最早將匹配追蹤算法應用于地震信號的分析中。與短時傅里葉變換[3]、連續(xù)小波變換[4]和S變換[5]等常規(guī)時頻分析方法相比,匹配追蹤算法的時頻分辨率相對較高,是一種高精度的地震信號分解與重構算法。但傳統(tǒng)的匹配追蹤算法計算速度非常慢,不適用于大型數(shù)據(jù)的處理過程。針對匹配追蹤算法的計算量龐大和運算效率極低的問題, 國內外專家學者提出了許多改進算法。Liu等[6-7]對Mallat的匹配追蹤算法作了改進,率先提出了動態(tài)匹配追蹤算法。張繁昌等[8]提出雙參數(shù)動態(tài)子波庫匹配追蹤快速算法,既減少掃描參數(shù)的個數(shù),又采用多原子分解策略,提高了匹配分解的效率。劉漢卿等[9]針對瞬時頻率存在負值的情況提出了利用連續(xù)相位求取瞬時頻率,基本解決了由主值相位和解纏相位求導獲得的瞬時頻率存在負值這一問題,但該方法也存在一定的問題,利用連續(xù)相位直接求取的瞬時頻率容易受到噪聲影響產生波動,出現(xiàn)一些頻率異常值。地震屬性大都是瞬時屬性[10],像瞬時頻率、瞬時振幅、瞬時傾角等,這些描述地震頻率特征的屬性依附于每個信號點的瞬時值。瞬時頻率屬性受相移、噪聲及其他因素的影響,其可解釋性差,很難反映微小的地質變化特征。局部屬性描述地震特征并不依賴每一個數(shù)據(jù)點的瞬時值,而是參考這個數(shù)據(jù)點的局部鄰近值,采用最小二乘估計策略,引進整形正則化對數(shù)據(jù)進行平滑處理,保證了局部頻率的穩(wěn)定性和可靠性。局部頻率在各個領域都得到廣泛的應用,陳常樂[11]應用局部屬性進行火山巖識別技術的研究,陳雙全等[12]應用局部頻率屬性檢測天然氣藏?；谝陨涎芯勘尘?筆者針對匹配追蹤的瞬時頻率存在負異常的情況進行深入的研究,引入局部頻率計算模式[13-14],將局部頻率作為約束動態(tài)子波庫的頻率搜索參數(shù)。

1 基于局部頻率約束的動態(tài)快速匹配追蹤方法

1.1 動態(tài)快速匹配追蹤算法

匹配追蹤算法是一種基于子波庫掃描的地震信號分解方法。假定被表示的信號為f,過完備原子庫為D,H表示Hilbert空間,D由一組向量{gγ1,gγ2,…,gγk}構成,其中每個向量可以稱為原子(atom),其長度與被表示信號f的長度相同,并且對這些向量已進行歸一化處理,即‖gγ‖=1。

令信號f∈H,第1次迭代分解從冗余原子庫D中選取一個與信號f最佳匹配的原子,即滿足下式:

〈f,gγ1〉=supi∈(1,2,…,k)〈f,gγi〉.

(1)

式中,i為字典矩陣D的列索引;gγ1為第1次迭代搜索得到的最佳匹配原子;〈f,gyi〉代表信號與最佳匹配子波的內積運算;sup為〈f,gγi〉的上確界。

這樣,經過一次迭代,信號f就被分解為在最佳匹配原子gγ1的垂直投影分量和信號殘差兩部分,即

f=〈f,gγ1〉gγ1+R1f.

(2)

式中,R1f為第一次迭代后得到的殘差信號。匹配追蹤算法是一個不斷迭代過程[15],將殘差信號R1f繼續(xù)投影到原子庫中,找到一個最佳匹配的原子,如此迭代下去。第n+1次迭代分解過程可以用下式表示:

Rnf=〈Rnf,gγn+1〉gγn+1+Rn+1f.

(3)

式中,Rnf和Rn+1f分別為第n次和第n+1次迭代完成后得到的殘差信號;gγn+1為第n+1次迭代搜索得到的最佳匹配原子,其中gγn+1滿足下式:

〈Rnf,gγn+1〉=supi∈(1,2,…,k)〈Rnf,gγi〉.

(4)

不斷進行上述的迭代過程,直到匹配次數(shù)達到預設的迭代次數(shù)或信號殘差能量小于能量閾值,則迭代停止。經過m步分解后,信號f被分解為m個最佳匹配的時頻原子之和:

(5)

其中

R0f=f.

盡管傳統(tǒng)的匹配追蹤算法能夠較好地表示信號,但因迭代分解的每一步都要將信號殘差與冗余原子庫中每一個原子進行內積運算,導致計算量過于龐大。因此專家學者們致力于在不損害匹配追蹤算法的高分辨率的前提下,減少算法的計算量,使得算法的實際應用范圍更廣泛。Liu等[6]提出動態(tài)匹配追蹤算法,利用復地震道分析技術得到信號的主頻、中心時間和相位等先驗信息,只需要確定這3個控制參數(shù)即可快速創(chuàng)建原子庫以及實現(xiàn)對最優(yōu)時頻原子的掃描,由于每次迭代的時頻原子的搜索范圍會隨著最大包絡振幅位置發(fā)生變化,因此稱該搜索策略為動態(tài)最優(yōu)搜索方式,這使得貪婪的匹配追蹤算法運算效率得到顯著提高。

在對信號進行匹配分解過程中,每次迭代只在信號包絡振幅最大值處進行,這樣的分解策略使得每次迭代歷經搜索范圍內所有的原子,僅得到一個最佳匹配的時頻原子,影響了信號的收斂速度。張繁昌等[8]提出了雙參數(shù)動態(tài)子波庫匹配追蹤快速算法,不僅減少匹配子波的控制參數(shù),且采用多原子分解的策略,即每次迭代并不是只在包絡振幅最大值處進行匹配搜索,而是在滿足一定條件的極大值位置處都進行搜索,這樣一次迭代就能得到若干個匹配子波。各匹配子波的振幅可以通過阻尼最小平方的方法來確定,大大提高了分解的效率。

在迭代分解前,動態(tài)快速匹配追蹤算法首先考慮原始信號中所包含的瞬時頻率、瞬時相位、時間等先驗信息,進一步約束其搜索半徑,因此每一次的搜索方式可以簡單表述為

gγi={u0,ω∈U[ω(u0),δω],φ∈U[φ(u0),δφ]},i=1,2,…,M.

(6)

式中,M為每次迭代得到的最優(yōu)時頻原子的個數(shù);u0為時間中心,ω和φ表示地震信號瞬時頻率及瞬時相位的先驗信息[16];U[ω(u0),δω]、U[φ(u0),δφ]為頻率和相位的搜索鄰域;δω、δφ為相應參數(shù)搜索半徑。

經過k次迭代后,得到信號殘差Rkf,在第k+1次迭代中,假設共得到M個最優(yōu)時頻原子gγi(i=1,2,…,M),則第k次迭代得到殘差信號Rkf可以表示為

(7)

式中,ai(i=1,2,…,M)為第k+1次迭代得到的各時頻原子的振幅。

1.2 局部頻率約束方法

Taner等[17]提出復地震道分析技術,首次將瞬時頻率引入到地震領域,之后瞬時頻率成為一種重要的地震屬性。后來許多學者提出瞬時頻率的近似算法,Barnes等[18]提出了3種瞬時頻率的近似計算方法;高靜懷等[19]提出小波變換域瞬時頻率分析方法;尹繼堯等[20]提出基于TK能量的最大峰值瞬時頻率計算方法等。

考慮到傳統(tǒng)瞬時頻率的不穩(wěn)定性,本文中提出用局部頻率替代瞬時頻率并應用于動態(tài)快速匹配追蹤算法中。局部頻率算法不僅能處理噪聲強的數(shù)據(jù),即使是部分缺失的數(shù)據(jù)也能高效準確地處理。局部頻率算法應用的是整形理論(也稱數(shù)據(jù)規(guī)則化),整形理論是一種通過對矩陣進行操作包括逆矩陣計算的數(shù)據(jù)平滑算法。對于地震數(shù)據(jù)部分微弱或缺失的情況,局部頻率算法可以從鄰近的時間窗口內提取有效信息進行處理并得到合理的頻率結果。

假設x(t)代表地震道,是時間t的函數(shù),則相應復數(shù)道表示為

c(t)=x(t)+ih(t).

(8)

其中h(t)為x(t)的希爾伯特變換,也可以用振幅A(t)和瞬時相位φ(t)來表示復數(shù)道,如下式:

c(t)=A(t)eiφ(t).

(9)

其中

φ(t)=arctan[h(t)x-1(t)].

式中，A(t)為地震道的包絡振幅;φ(t)為瞬時相位。

瞬時頻率是瞬時相位的導數(shù),也就是瞬時相位的變化率,因此地震道的瞬時頻率屬性ω(t)的數(shù)學表達式為

(10)

將瞬時頻率的數(shù)學表達式寫成向量表達式的形式,把時間記錄看成是一個向量,通過定義一個對角矩陣和兩個向量,即

(11)

式中,k為地震道的采樣點數(shù)。利用式(11)將方程(10)寫成矩陣形式,即

W=D-1N.

(12)

式中,W為瞬時頻率向量ω(t);N為由式(10)中的分子n(t)構成的向量;D為由(10)式中的分母d(t)構成的對角矩陣算子。注意到,矩陣D中很多元素很小甚至是零,為了避免式(12)出現(xiàn)“除以零”導致計算錯誤的現(xiàn)象,在分母上加一阻尼項ε,因此方程(12)變?yōu)?/p>

Winst=(D+εI)-1N.

(13)

式中,I為單位運算符。這里阻尼項ε并不能讓瞬時頻率變得穩(wěn)定,但可以減弱瞬時頻率受噪音和不穩(wěn)定性因素的影響。

局部頻率屬性將式(13)作為線性反演的正則化形式,用一個一般化的正則化算子R代替單位矩陣I來改變正則化形式,因此局部頻率的定義如下:

Wloc=(D+εR)-1N.

(14)

式中,R為正則化算子,以確保局部頻率的連續(xù)性和光滑性,使得計算結果具有更高的信噪比。也可以選取整形正則算子S[21],對式(14)做進一步的約束,求解得到下式:

Wloc=(λ2I+S(D-λ2I))-1SN.

(15)

式中,λ為尺度因子,通常選自矩陣D的最小二乘范數(shù)。利用迭代反演方法求解時,通過λ進行縮放,可以保留原始物理維度,并能夠加快收斂速度。

為了驗證局部頻率的可靠性,設計了3個用于比較信號頻率屬性的測試信號。圖1(a)為一個頻率由10 Hz增大到60 Hz的線性調諧信號,圖1(b)為由頻率10 Hz的Morlet子波合成的理論信號,圖1(c)為一道實際地震信號。圖1(d)、(e)、(f)展示了依次對應于測試信號的瞬時頻率,由于瞬時相位的突變不可避免地出現(xiàn)負異?，F(xiàn)象,不能直接用于確定動態(tài)匹配追蹤控制參數(shù)的搜索范圍,圖1(g)、(h)、(i)展示了3個測試信號的局部頻率,可以明顯看到局部頻率的結果比較精確地貼近于真實值,克服了瞬時頻率出現(xiàn)負值的問題。

圖1 不同信號的瞬時頻率和局部頻率對比Fig.1 Comparison of instantaneous frequency and local frequency of different signals

驗證局部頻率的抗噪性。同樣采用圖2中的兩個理論信號,加入10%的噪聲,如圖2(a)和(d)所示。從兩個理論信號的瞬時頻率圖2(b)和(e)可以看出瞬時頻率受噪聲影響波動大,極不穩(wěn)定,而局部頻率圖2(c)和(f)的結果與理論值相符,基本不受噪聲的影響,比較穩(wěn)定。

對比以上的瞬時頻率與局部頻率的運算結果可以明顯地看出,局部頻率的結果比較精確地貼近于真實值。將局部頻率引入動態(tài)快速匹配追蹤方法中,用局部頻率替代傳統(tǒng)的以信號的瞬時頻率作為動態(tài)參數(shù)的掃描范圍。不僅避免了常規(guī)方法的負頻率異常現(xiàn)象,并且抗噪性也得到改善,計算結果穩(wěn)定。圖3為基于局部頻率的動態(tài)快速匹配追蹤算法的流程。

圖2 加入10%噪聲的理論信號的瞬時頻率和局部頻率對比Fig.2 Comparison of instantaneous frequency and local frequency of different signals with 10% noise

圖3 基于局部頻率約束的動態(tài)快速匹配追蹤方法流程Fig.3 Flowchart of dynamic and fast matching pursuit method based on local frequency constrains

2 理論信號測試

2.1 一維理論驗證

為了驗證基于局部頻率約束的動態(tài)快速匹配追蹤方法的可行性和穩(wěn)定性,本文中分別設計了一維和二維地震模型進行驗證。首先用11個不同頻率、相位和振幅的Morelet子波合成一道理論信號,構建一維地震模型,對其用基于局部頻率約束的動態(tài)快速匹配追蹤方法進行分解,迭代分解結果如圖4所示。可以看出理論信號可以完全由匹配子波重構得到,殘差曲線為零值,對其增加5%的噪聲,得到如圖5的地震模型。從分解結果可以看出,在加噪聲的情況下依舊能夠得到與無噪聲情況一致的匹配分解結果,迭代分解結果比較穩(wěn)定。

圖4 一維合成信號的迭代分解結果Fig.4 Iterative decomposition results of one-dimensional synthesized signal

圖5 加5%噪聲的一維合成信號的迭代分解結果Fig.5 Iterative decomposition results of one-dimensional synthesized signal with 5% noise

2.2 二維理論驗證

為了進一步驗證基于局部頻率約束的動態(tài)快速匹配追蹤方法的可靠性,建立了如圖6(a)的一個二維縱波阻抗模型。模型的復雜性有助于驗證本方法基于局部頻率的準確性。圖6(b)為利用28 Hz零相位雷克子波合成的地震數(shù)據(jù),對其用基于局部頻率約束的快速匹配追蹤方法進行匹配分解,圖6(c)為重構的地震記錄,并得到殘差剖面(圖6(d))。對圖6(b)合成地震記錄加入10%的噪聲得到圖7所示的結果。

圖6 無噪聲的合成信號迭代分解結果Fig.6 Iterative decomposition results of synthesized signal without noise

圖7 加入10%噪聲的合成信號迭代分解結果Fig.7 Iterative decomposition results of synthesized signal with 10% noise

從圖6的迭代分解結果可以看出,基于局部頻率的動態(tài)快速匹配追蹤算法能夠很好地對理論數(shù)據(jù)進行迭代分解,重構的地震剖面與原始地震剖面一致,殘差基本為零。對地震記錄加噪聲后(圖7),殘差剖面所示基本為噪聲信號，不含有效的反射信息，證明該算法具有一定的抗噪性。

3 實際資料處理

實際資料來自中國東部某勘探研究區(qū),圖8為某一實際地震剖面,為便于了解砂體的尖滅特征,將砂體的形態(tài)特征展布在地震剖面上,紅色的橢圓框為砂體真正尖滅點位置,圖9為用基于局部頻率的MP時頻譜分析的尖滅點識別結果。檢測結果與實際尖滅點位置一致,很好地驗證了本文中提出的改進匹配追蹤快速算法在砂體尖滅線識別方法方面的可行性。

圖8 實際地震剖面Fig.8 Real seismic profile

圖9 尖滅點識別結果Fig.9 Tip point recognition result

為了對比本文中提出的改進的動態(tài)快速匹配追蹤算法與傳統(tǒng)的三參數(shù)動態(tài)快速匹配追蹤算法,使用相同的迭代終止條件,對圖8所在研究區(qū)的一道地震記錄進行匹配追蹤分解,從圖10中可以清楚地看到,本文中提出的基于局部頻率約束的動態(tài)快速匹配追蹤算法提高了匹配分解的計算速度。

圖10 兩種匹配追蹤算法計算時間對比Fig.10 Comparison of computing time of two methods

4 結束語

針對瞬時頻率存在負異常現(xiàn)象且不穩(wěn)定的問題,本文中發(fā)展了基于局部頻率約束的動態(tài)快速匹配追蹤算法,局部頻率引進整形光滑算子,其計算過程并不是依附于每一個數(shù)據(jù)點的瞬時值,而是參考這個數(shù)據(jù)點的局部鄰近值求取的,因此即使是在部分信號微弱或缺失的情況下,局部頻率也能從鄰近的時窗內提取有效信息并得到合理的頻率值。相比于瞬時頻率由于瞬時相位求導存在負頻率現(xiàn)象,且受噪聲影響大,計算結果極不穩(wěn)定,本文中通過設計理論信號與實際信號驗證了局部頻率的可行性,其抗噪性得到增強,計算結果穩(wěn)定,可以直接將其應用于匹配追蹤的運算過程。在確定動態(tài)快速匹配追蹤先驗信息的搜索范圍時,以局部頻率替換瞬時頻率,加快了運算效率。最后將基于局部頻率約束的動態(tài)快速匹配追蹤方法用于識別砂體尖滅點,在單頻瞬時譜剖面上可以清楚看到砂體尖滅位置,且與實際位置一致,且運算速率明顯提高,從而驗證了本文中提出的改進匹配追蹤算法的可行性。