基于局部峰值約束的快速正交匹配追蹤地震數(shù)據(jù)分解方法

杜澤源, 楊 森, 陶永慧, 蔡杰雄, 何兵紅

(1. 中國石油化工股份有限公司 石油物探技術(shù)研究院，南京 211103；2. 中國石油大學(xué)(華東) 地球科學(xué)與技術(shù)學(xué)院，青島 266580)

0 引言

地震數(shù)據(jù)包含了豐富的地下介質(zhì)信息，通過對地震數(shù)據(jù)分解可以得到多個含有不同信息的地震資料，有利于提高對地下不同介質(zhì)情況及數(shù)據(jù)特征的認(rèn)識與辨別。地震數(shù)據(jù)分解是分析地震數(shù)據(jù)結(jié)構(gòu)、屬性特征，提取地質(zhì)信息的有效途徑。在地震信號處理領(lǐng)域，地震數(shù)據(jù)稀疏分解可以用于數(shù)據(jù)去噪、數(shù)據(jù)重建[1]以及時頻分析[2]和地震反演等問題[3]。數(shù)據(jù)分解可以通過多種途徑實現(xiàn)，傅里葉變換將時間域信號變換到頻率域，使人們可以更好地分析信號的特點，在頻率域?qū)π盘栠M(jìn)行處理在很多情況下可以帶來很多便利[4]。但是傅里葉變換只能對信號整體進(jìn)行分析，不能獲得局部的信息[5]，這對于非平穩(wěn)地震信號的處理是不利的。小波變換以其多分辨率的特點能夠得到更好的信號分解，小波分析的思想被廣泛應(yīng)用到地震數(shù)據(jù)分析中[6-7]。Mallat[8-9]提出了基于塔式多分辨率分析的小波分解與信號重構(gòu)的快速算法(Mallat算法)；許康生[10]研究了基于小波變換的地震信號多尺度分解，證明了重建地震信號無畸變和虛假信號產(chǎn)生，而且多尺度分解可用于地震記錄的弱震相時間信息的提取；陳飛等[11]利用小波分析進(jìn)行地震信號多尺度分解及去噪處理，通過數(shù)據(jù)測試證明能夠得到較好的去噪效果；杜麗英[12]提出小波包最佳基地震數(shù)據(jù)分解與重建，過最佳小波基的選擇，使在空間域局部性很好的情況下，頻域的局部性也達(dá)到最佳狀態(tài)；周靜毅等[13]研究了基于小波變換的地震奇異性分析，利用小波變換良好的時間-頻率局部化特征，通過地震數(shù)據(jù)奇異性屬性分析來劃分地層旋回。但小波變換母小波的選擇以及分解的層數(shù)，對信號多尺度分解的結(jié)果影響較大，有時候難以獲得理想的尺度分解數(shù)據(jù)。

稀疏表示理論為信號分析處理帶來了一次革命性的推動，它的核心是通過稀疏分解，借助字典中少量原子的線性組合形式對原信號所涵蓋的信息進(jìn)行表達(dá)。Mallat等[14]提出使用過完備原子庫用于信號稀疏分解的匹配追蹤算法。匹配追蹤算法具有自適應(yīng)選擇匹配原子的優(yōu)點[15-16]，在地震信號處理領(lǐng)域的應(yīng)用發(fā)展迅速。何胡軍等[17]研究了基于匹配追蹤算法子波分解技術(shù)識別薄互層儲層，該技術(shù)能較好地識別薄互層儲層，分辨出儲層在空間的展布規(guī)律；張繁昌等[18]利用基于匹配追蹤的地震瞬時譜分解方法實現(xiàn)了三角洲砂巖尖滅先的識別。匹配追蹤是一種貪婪算法，在信號分解過程中容易出現(xiàn)重復(fù)選擇原子的情況，影響計算效率，許多學(xué)者也提出了改進(jìn)算法。Wang[19]提出了多道匹配追蹤算法，并應(yīng)用于地震數(shù)據(jù)的分解；汪瑞良等[20]研究了快速動態(tài)匹配追蹤算法，實現(xiàn)了匹配追蹤高分辨率時頻譜計算。正交匹配追蹤(Orthogonal Matching Pursuit，OMP)算法在每次迭代計算過程中保證選擇的原子相互正交，計算效率更高[21]?？焖儆行У剡M(jìn)行匹配追蹤計算是一個挑戰(zhàn)性問題，對于最初的正交匹配追蹤算法，由于原子字典規(guī)模很大，對原子遍歷搜索，計算成本依然很高，對于正交匹配追蹤算法的優(yōu)化依然是研究的熱點。

針對正交匹配追蹤在通過內(nèi)積實現(xiàn)子波與信號匹配時不能保證子波與信號具有相同的位置及振幅的問題，避免在子波匹配時對字典內(nèi)所有子波都進(jìn)行掃描而大大降低計算效率。筆者采用局部峰值約束的正交匹配追蹤方法來解決這一問題，構(gòu)建信號與子波峰值距離的目標(biāo)泛函，通過快速掃描算法進(jìn)行高效的求解計算。模型測試結(jié)果表明，該方法可以高效實現(xiàn)地震數(shù)據(jù)的高精度分解與重構(gòu)，驗證了本文算法的有效性和高效性。

1 理論方法

1.1 正交匹配追蹤求解方法原理

匹配追蹤是基于稀疏表示理論發(fā)展起來的，稀疏表示為地震信號處理及地震反演等方面帶來了很多便利。地震信號的稀疏表示可以通過幾個原子(基本波形)的線性加權(quán)組合表示：

(1)

其中：s(t)為地震信號；ω(t)為信號原子(基本子波)；ωk(t)是從一個原子庫(子波字典)中選取的；原子庫Φ=[ω1(t),ω2(t),…,ωK(t)]是由信號原子為列向量構(gòu)成的維度為N×K的矩陣；N為子波長度；αk為信號原子對應(yīng)的系數(shù)；α=[α1,α2,…,αK]為信號在原子庫Φ中的表示系數(shù)。

匹配追蹤算法的主要思想是，首先建立一個過完備原子庫，在每次迭代求解過程中，從原子庫中選出與信號最為匹配的原子，來稀疏逼近信號，進(jìn)而求出信號與所選原子的殘差；然后從原子庫中選出與殘差信號最匹配的原子，同時計算新的殘差；重復(fù)多次迭代計算以上過程，直到殘差小于設(shè)定的閾值或迭代次數(shù)達(dá)到設(shè)定的上限迭代次數(shù)時，將選出的最佳原子進(jìn)行線性組合就可以表示原始信號。但是匹配追蹤存在一個問題，在匹配分解過程中很可能會重復(fù)選擇原子，這嚴(yán)重影響計算效率。正交匹配追蹤算法要求每次求得的信號殘差與之前選出的所有原子相正交，這樣在匹配時就會避免重復(fù)地選擇原子，能夠更快地收斂到全局極值[21-22]。

所述正交匹配追蹤算法的目標(biāo)泛函可表示為：

(2)

其中，σ為迭代終止閾值。

在稀疏求解問題中，如果原子ωk(t)正交，系數(shù)αk可以用信號s(t)與原子ωk(t)的內(nèi)積來計算得到，即：

(3)

匹配迭代的過程為：首先從匹配原子庫Φ中選出與給定地震信號s(t)匹配最佳的原子ω1(t)，匹配最佳的要求則是它們的內(nèi)積[s(t),ω1(t)]為匹配原子庫中所有原子與該地震信號s(t)內(nèi)積最大的一個：

|[s(t),ω1(t)]|>|[s(t),ωk(t)]|,k≠1

(4)

即，第一次迭代后，信號s(t)被分解為最佳原子ω1(t)上的分量和信號殘差R1(t)，

[s(t),ω1(t)]ω1(t)+R1(t)

(5)

然后對殘差信號R1(t)重復(fù)上述相同的匹配分解過程：從Φ中篩選出與R1(t)最為匹配的第二個原子，即第二次迭代可表示為：

[R1(t),ω2(t)]ω2(t)+R2(t)

(6)

不斷重復(fù)上述過程，第n次分解后，獲得殘差信號Rn(t)：

[Rn-1(t),ωn(t)]ωn(t)+Rn(t)

(7)

所以經(jīng)過n次分解后，信號可表示為式(8)。

(8)

其中R0(t)為原始地震信號s(t)。式(8)表達(dá)的意思即為經(jīng)過n次迭代分解計算后，原始信號s(t)可以用n個匹配原子與殘差Rn(t)的合成來表示。

由于殘差‖Rn(t)‖的衰減特性，對有限長度的信號，‖Rn(t)‖隨著迭代次數(shù)n的增大而衰減至小于閾值σ，即信號的稀疏分解跟信號的近似程度足夠好或殘差信號滿足設(shè)定的閾值時，迭代停止。一般，用少數(shù)M個原子就可以表示信號主要成分，滿足數(shù)據(jù)處理的需求，且閾值σ是很小的數(shù)，則地震信號(可稀疏表示)為最佳原子的線性組合式(9)。

(9)

1.2 子波字典構(gòu)建

匹配追蹤算法能夠快速有效計算的關(guān)鍵之一是，選取合適的信號原子構(gòu)建原子庫。經(jīng)典的匹配追蹤算法是基于Garbor原子來構(gòu)建過完備匹配原子庫。為適應(yīng)地震信號不同的結(jié)構(gòu)特征，對于原子庫可以由一系列能夠表達(dá)地震信號特征的子波波形構(gòu)成，即過完備匹配子波庫是一系列匹配子波構(gòu)成的子波矩陣。

對于地震數(shù)據(jù)處理，子波字典實質(zhì)上就是一系列基本子波ωδ(t)進(jìn)行適當(dāng)?shù)臅r移、調(diào)制和相位變化得到的子波矩陣。若ωδ(t)為子波庫中滿足匹配條件的一個子波，其控制參數(shù)用δ(τ,f,φ,c)來表示，子波函數(shù)記為：

ωδ(t)=ω(t-τ)exp{i[πf(t-τ)/c+φ]}

(10)

其中：τ表示時移；f表示中心頻率；φ表示相位；c表示尺度因子；i表示虛數(shù)單位。

子波字典的構(gòu)建有很多子波可選擇，可以根據(jù)所分析信號的特征及處理目標(biāo)進(jìn)行選擇。對于地震數(shù)據(jù)來說，可以選擇Ricker子波、高斯函數(shù)或Morlet小波等基本子波構(gòu)建子波字典。Ricker子波是比較常用的子波，由于其波形簡單，延遲時間短，時頻分辨率較高，控制參數(shù)改變靈活，易于構(gòu)建波形豐富的子波字典，充分表達(dá)地震數(shù)據(jù)[23]。另外，Morlet小波有較好的尺度變化能力，在時間和頻率的局部化之間有著很好的平衡，切與地震信號有較好的相似性，也可以選取Morlet小波創(chuàng)建子波字典[24]。

1.3 局部峰值約束正交匹配追蹤

正交匹配追蹤的目標(biāo)是最大化所選匹配子波與信號殘差之間的內(nèi)積，然而在內(nèi)積實現(xiàn)時往往只是在全局時間尺度上度量原子與殘差之間的信號相似性，并不能保證原子波峰與殘差波峰具有相同的位置。

為了解決這一問題，采用局部峰值約束的方法來確定匹配子波參數(shù)，將最小化信號殘差波峰與原子峰值距離為目標(biāo)泛函。為了更有效地對信號及原子波峰進(jìn)行掃描，首先掃描信號包絡(luò)，獲取其局部極大值來限定信號的掃描范圍，然后再尋找信號峰值并掃描原子進(jìn)行匹配。因此定義了以下目標(biāo)函數(shù)：

(11)

具體地，對信號進(jìn)行Hilbert變換，構(gòu)造復(fù)數(shù)信號，通過時間序列確定信號采樣率和信號長度，利用復(fù)數(shù)信號求取信號的瞬時包絡(luò)與瞬時相位，則τi可粗略估計為瞬時包絡(luò)局部極大值峰值處所對應(yīng)的時間ti。在復(fù)數(shù)域內(nèi)對子波庫中的原子進(jìn)行頻率參數(shù)的掃描，對于某次迭代確定出的一組中心時間τb，得到信號在這些時間處的頻率先驗值fb，假設(shè)基本子波為零相位子波，則可利用τb和fb確定若干的先驗信息點(τb,fb,0)。以這些先驗信息點為中心，選取頻率掃描范圍，同時搜索多個原子，來提高匹配追蹤效率，確定匹配子波的波峰頻率fi。

在確定了最佳匹配子波ω(t)的中心時間和主頻后，需要確定這些匹配子波的振幅和相位。根據(jù)殘差能量盡可能小的原則，可以利用殘差能量阻尼最小二乘方法[25]計算這些匹配子波的振幅和相位信息。

1.4 地震數(shù)據(jù)多尺度分解

地震信號稀疏處理一般可以表示為線性形式[1]：

dobs(t)=Fs(t)+o(t)

(12)

其中:dobs(t)為觀測地震數(shù)據(jù)；o(t)為噪聲；F為有界線性算子，對于不同的數(shù)據(jù)處理問題有不同的表達(dá)形式。對于本研究中，F(xiàn)表示多尺度分解算子；s(t)為稀疏求解數(shù)據(jù)。

地震資料可以分解成不同的分量(或尺度)。重建的地震信號可以表示為地震信號分量的線性疊加：

s(t)=sδ(f1)(t)+sδ(f2)(t)+…+sδ(fC)(t)

(13)

其中，δ(fi)代表不同的頻帶尺度。子波字典中子波由δ(τ,f,φ,c)參數(shù)控制包含頻率因子，將頻帶尺度的子波歸為一類，則子波字典可以分為相應(yīng)頻率尺度下Φδ(fi)。

地震信號稀疏表示的解為s(t)=Φα，則稀疏約束反演問題在相應(yīng)頻帶尺度下可表示為求解以下泛函：

(14)

對地震數(shù)據(jù)稀疏分解則可以直接提取不同的頻率分量，利用正交匹配追蹤算法求解多尺度分解目標(biāo)泛函(14)，即可將地震數(shù)據(jù)分解為不同頻率尺度下的分量，基于本方法的數(shù)據(jù)分解有較好的時頻分辨率，能夠較好地服務(wù)于多尺度數(shù)據(jù)分析或多尺度反演。

2 模型測試

為驗證文中正交匹配追蹤算法的效果，設(shè)計幾個基本波形疊加得到單道信號，如圖1(a)所示。然后利用正交匹配追蹤信號分解算法對合成地震信號進(jìn)行分解，將分解后的信號進(jìn)行重構(gòu)，得到重構(gòu)地震信號，見圖1(b)。

對合成信號速度分解求解迭代次數(shù)為23次，從圖1可以看出,重構(gòu)信號與原始合成信號波形基本一致，圖1(c)為原始信號與重構(gòu)信號的殘差，信號殘差小于2‰，滿足計算精度要求。圖2為將原始信號與重構(gòu)信號匹配對比，也可以看出原始信號與重構(gòu)信號基本吻合，由此說明本文分解方法能夠有效地匹配出最佳子波，并精確計算子波的振幅及相位。通過模型測試驗證了本文稀疏分解重構(gòu)算法能夠?qū)⑿盘柧_分解重構(gòu)，效果較好。

圖1 合成信號分解重構(gòu)的結(jié)果Fig.1 The result of decomposition and reconstruction of synthetic signal(a)合成信號;(b)重構(gòu)信號;(c)信號殘差

圖2 原始信號與重構(gòu)信號對比Fig.2 Comparison of original signal and reconstructed signal

匹配追蹤有較好的信號分解與重構(gòu)精度，同時具有較高的時頻分辨率?；谄ヅ渥粉櫟亩虝r傅里葉變換可以對每一個匹配子波進(jìn)行變換，時窗調(diào)節(jié)比較靈活，相對于短時傅里葉變換時窗固定的特點，匹配追蹤算法的時頻分辨率更高。Wigner-Ville分布雖然相對于短時傅里葉變換時頻分辨率提高，但是易產(chǎn)生交叉干擾項。分別用短時傅里葉變換、Wigner-Ville分布及基于匹配追蹤的方法分析信號的時頻譜，結(jié)果如圖3所示。

圖3 不同算法下的信號時頻譜Fig.3 Time spectrum of signals under different algorithms(a)短時傅里葉時頻譜;(b)Wigner-Ville時頻譜; (c)基于OMP的短時傅里葉時頻譜;(d)基于OMP的Wigner-Ville時頻譜

從圖3可以看出,基于OMP的短時傅里葉變換時頻譜(圖3(c))相對于常規(guī)短時傅里葉變換時頻譜(圖3(a))分辨率更高。常規(guī)Wigner-Ville分布時頻譜(圖3(b))比常規(guī)短時傅里葉變換時頻譜分辨率高，但是產(chǎn)生了交叉項，而通過正交匹配追蹤算法改進(jìn)的Wigner-Ville分布時頻譜(圖3(d))消除了干擾項，同時保留了較高的時頻譜分辨率，能量更加聚焦。通過對比說明基于匹配追蹤的時頻譜具有較高的時頻分辨率和能量聚焦特性，能夠準(zhǔn)確表示地震信號的局部頻譜特征。

3 實際資料分析

采用某實際資料分析本文方法應(yīng)用效果，選取一條二維剖面(圖4)。地震數(shù)據(jù)共131道，每道數(shù)據(jù)為0.5 s，采樣間隔為2 ms。使用本文正交匹配追蹤數(shù)據(jù)分解方法對圖4中地震數(shù)據(jù)進(jìn)行分解與重構(gòu)，圖5展示了原始地震剖面(圖5(a))與重構(gòu)后地震剖面(圖5(b))，通過對比圖5(a)和圖5(b)可以看出，重構(gòu)后的剖面與原始剖面基本完全一致，基本恢復(fù)了原始數(shù)據(jù)中包含的細(xì)節(jié)部分。

圖4 原始地震剖面Fig.4 Original seismic section

圖5 原始地震剖面與OMP分解重構(gòu)剖面Fig.5 Original seismic section and OMP decomposition reconstruction section(a)原始地震剖面；(b)OMP分解重構(gòu)剖面

抽取地震剖面中第75道數(shù)據(jù)，如圖6(a)所示，在0.33 s處有一含氣儲層。對地震道進(jìn)行時頻分析，圖6(b)為其時頻譜，可以看出在含氣儲層位置處譜能量較為聚焦，時頻分辨率較高，可以清楚地區(qū)分上下地層。通過時頻譜可以非常直觀的了解含氣儲層地震響應(yīng)的頻帶特征。

圖6 第75道地震數(shù)據(jù)及其時頻譜Fig.6 The 75th trace seismic data and its time-frequency spectrum(a)第75道地震數(shù)據(jù)；(b)地震道時頻譜

使用本文多尺度分解方法，對地震剖面進(jìn)行尺度分解。由圖7可以看出，分解后的地震剖面符合原始剖面中的地質(zhì)特征，同相軸與原始剖面中形態(tài)一致，沒有假頻產(chǎn)生，頻譜的混疊也能有效分離，也可以看出不同的地層有不同的頻帶相應(yīng)特征。在第75道位置處有一口井，聲波測井?dāng)?shù)據(jù)見圖7，分解后的地震同相軸指示地層與測井速度分層相吻合。因此，可以不同地層地震響應(yīng)頻帶特征為依據(jù)，使用分解后的地震數(shù)據(jù)在目標(biāo)優(yōu)勢尺度下進(jìn)行數(shù)據(jù)分析或地震反演，會得到更好的效果。通過地震數(shù)據(jù)分解，可以為面向各種目標(biāo)的數(shù)據(jù)處理或地震反演問題選擇提供可靠的數(shù)據(jù)分量。

圖7 地震數(shù)據(jù)多尺度分解剖面Fig.7 Multi-scale decomposition sections of seismic data(a)主頻為10 Hz的單頻剖面；(b)主頻為16 Hz的單頻剖面； (c)主頻為22 Hz的單頻剖面；(d)主頻為28 Hz的單頻剖面

4 結(jié)論

匹配追蹤算法是一種貪婪算法，難以避免原子的重復(fù)選擇，正交匹配追蹤算法要求每次迭代后的殘差信號和之前選出的所有原子相正交，這樣在匹配時就會避免重復(fù)地選擇原子。但是正交匹配追蹤在內(nèi)積的實現(xiàn)時，往往只是在全局時間尺度上度量原子與殘差之間的信號相似性，并不能保證原子波峰與殘差波峰具有相同的位置，對信號原子的掃描也存在冗余。針對這一問題，筆者研究了基于局部峰值約束的正交匹配追蹤方法，并在原子掃描時進(jìn)行優(yōu)化，避免過多的選取原子，影響效率。在算法實現(xiàn)的基礎(chǔ)上，通過模型數(shù)據(jù)與實際資料測試證明了本文數(shù)據(jù)分解與重構(gòu)方法的有效性。

地震數(shù)據(jù)稀疏表示分解方法能夠高效精確地將地震數(shù)據(jù)分解為不同頻帶尺度，在數(shù)據(jù)稀疏分解時夠克服頻譜混疊和波的干涉效應(yīng)，提高時頻分辨率，并且不會引入假頻。得到的多尺度(分量)地震數(shù)據(jù)，可以為數(shù)據(jù)分析或多尺度地震反演提供可靠的數(shù)據(jù)支持。