韓海英，王志章，王宗俊,張雨晴

(1.中國石油大學(北京)地球科學學院，北京102249;2.中國石油天然氣股份有限公司勘探開發(fā)研究院，北京100083;3.中海油研究總院，北京100027)

近年來，信號時頻分解技術在地震資料處理和解釋中得到了廣泛的應用，如高分辨率處理[1]、烴類檢測[2-3]等。常規(guī)的時頻分解(如窗口傅里葉變換和小波變換等)，時頻窗受不確定準則約束，時頻分辨率難以同時達到最佳。為了克服常規(guī)方法的缺點，Mallat等[4]和Qian等[5]提出了匹配追蹤分解算法(Matching Pursuit，MP)。MP算法基于最佳匹配原則，通過不斷迭代，從原子庫中選取與殘余信號結構最佳匹配的原子，將原始信號分解為一系列原子。由于在迭代尋優(yōu)過程中，時窗長度均由信號局部特征決定，因此這些最優(yōu)原子反映了信號的局部特征，能較好地描述時變地震信號的時頻分布特性。

經典的MP算法基于Gabor函數構建超完備子波庫，為了適應信號的不同結構特征，Liu等[6-7]先后提出了基于Ricker子波和Morlet子波的匹配追蹤算法。對研究薄互層儲層特征而言，Morlet小波分析薄層的效果欠佳;宋新武等[8]認為Ricker子波波形簡單，延遲時間短，收斂較快，基于Ricker子波匹配追蹤算法的頻譜分解技術在薄互層分析中具有較高的時頻分辨率，能較好地刻畫地質體形態(tài)。原子庫的合理選取，對于信號的匹配追蹤分解效果具有重要作用。

Ricker類地震子波是由Ricker子波演變得到的零相位理論子波，具有3個控制參量，旁瓣能量的收斂速度可變，延遲時間可控，子波庫波形更豐富，與地震子波的匹配較為靈活[9]。我們在研究子波波形隨3個控制參量變化關系的基礎上，對該子波作適當的簡化，建立參量與波形特征之間的半定量關系。最后，以簡化的子波構建原子庫，通過三步法原理[10]實現信號的匹配追蹤分解。

1 基本原理

1.1 Ricker類地震子波

Ricker類地震子波頻域解析式為[9]

(1)

對于確定的k，m越大，高頻能量相對越強;對于確定的m，k越大，高頻能量相對越弱，同時頻譜的有效頻帶變窄，導致子波信號的旁瓣變化復雜。文獻[9]建議k不宜過大。事實上，參量k和m都是控制波譜高低頻相對強弱的參數，一般情況下，相同的增(減)量，m對頻譜形態(tài)的改變較為平緩。

(2)

這里記C(f)為二參量Ricker類地震子波，簡稱C子波。其中，fm和c為待定參量;A為歸一化系數。通過簡化和改寫，將子波譜表述成Ricker譜(主頻為fm)的c次方(c=1時對應Ricker子波)。fm為波譜主頻，控制波譜的等比例伸縮;c控制波譜主頻兩側能量的衰減速度。

下面考察參量c與子波波形的關系。首先，針對連續(xù)相位數小于1的對稱波形(圖1中藍線)，我們定義主、旁瓣極值比PR和主、旁瓣寬度比WR兩個參量，用于描述子波的波形特征。

主旁瓣極值比(PR)：旁瓣極值與主瓣極值之比的絕對值;

主旁瓣寬度比(WR)：旁瓣寬度(極小值點間距)與主瓣寬度(零值點間距)之比。

圖1 對稱C子波示意圖解

圖2為C子波波形特征值(紅色為極值比PR線，藍色為寬度比WR線)隨參量c的變化曲線。隨著c的增加，主旁瓣極值比增加，旁瓣能量越強;主旁瓣寬度比增加，說明在相同主瓣寬度的情況下，旁瓣能量達到極值需要的時間變長。因此，對C子波有如下結論：參數fm控制子波的橫向伸縮，fm越大，子波越短;參數c控制子波的波形形態(tài)，c越大，子波的主旁瓣極值比、寬度比越大，連續(xù)相位數增加。當c→0時，C子波趨近于脈沖信號，PR→0，WR→1；當c→∞時，C子波趨近于單頻信號，PR→1，WR→2。

圖2 主旁瓣寬度比(藍)和極值比(紅)曲線

1.2 基于C子波的匹配追蹤

匹配追蹤是基于原子庫掃描的信號自適應分解。記D為進行信號分解的超完備子波庫，可由Hilbert空間的窗口函數通過擴張、調制和平移生成，定義D={wγ:γ=1,2,3,…}，wγ為單個原子，滿足歸一化條件‖wγ‖=1。記待分解信號為s(t)，經N次迭代分解如下：

(3)

式中：wγ,n為第n次迭代得到的匹配子波;an為wγ,n對應的振幅；RNs為迭代N次后的殘差。

匹配追蹤算法中子波庫的波形與原信號結構越相似，分解效果越好，這里選擇C子波來構建匹配子波庫。匹配子波wγ,n由參量γn={tn,fm,n,cn,φn}來控制，其中，tn，φn為子波的中心時間和相位;fm,n，cn為波譜控制參量。匹配追蹤算法基于不斷迭代的貪婪算法，常規(guī)算法的每一次迭代都要從子波庫中尋找最佳子波。目前，通常采用“三步法”[10]原理來提高分解速度。

(4)

(5)

3) 估算最優(yōu)子波的幅值。按(6)式得到最優(yōu)子波的幅值an：

(6)

每一次迭代，都將生成匹配子波anwγ,n和殘差信號Rn+1s。分解完成后，即可用提取的匹配子波近似表示原始信號，即

(7)

1.3 時頻表征

常規(guī)的匹配追蹤時頻表征通過各匹配子波的Wigner-Ville分布疊加實現，但只能給出振幅分布，無法得到相位信息。這里，采用(8)式所示的時頻表征方法[12]，對每個匹配子波的復譜求和，獲取信號的時頻譜：

(8)

式中：Wγ,n(f)是匹配子波wγ,n(t)的頻譜;env[wγ,n(t)]是wγ,n(t)的瞬時包絡。

圖3是單個C原子(主頻fm取50Hz，c取0.7)的時頻譜。該時頻表征法具有較高的時頻分辨率，能同時獲取振幅和相位等信息的時頻分布，并且避免了窗口截斷效應。

圖3 C原子時頻譜

2 模型試算

設計兩個合成信號進行基于C子波庫的匹配追蹤算法測試，時間采樣率均為1ms。

圖4為信號1基于C子波的匹配追蹤重構結果，其中,圖4a為原始信號，圖4b為重構信號，圖4c 為重構相對誤差，重構誤差保持在2%以內，對信號的損傷較小(重構精度與算法終止條件有關，精度越高，運行時間越長)。信號組成方式如圖5a 所示，其中，第1道為原始信號(圖4a)，為其余4道的疊加;圖5b為信號1的匹配追蹤時頻譜。圖6a為信號1的短時Fourier變換(STFT)譜，圖6b 為S變換譜。對比發(fā)現，匹配追蹤與常規(guī)時頻分解方法在分辨率上存在較大的差異。在300ms和400ms處，3個時頻譜上都有單個能量團出現，無論是時間尺度還是頻率尺度，匹配追蹤的分辨率都是最高的;在500ms處，3種方法都沒能較好地把兩個頻率分辨出來，但匹配追蹤的時間分辨率還是最高的;在600ms處，S變換和STFT譜的時間分辨率不夠，匹配追蹤則較好地識別出2個反射;在700ms和800ms處，匹配追蹤譜的時頻分辨能力明顯高于STFT譜和S變換譜，能反映更多的地層薄層信息。因此，匹配追蹤時頻譜具有較高的時頻聚焦性。

圖4 信號1重構圖解a 原始信號; b 重構信號; c相對誤差

信號2(圖7)由Ricker類理論子波與隨機反射系數褶積生成。為符合實際情況，由淺到深理論子波的主頻和帶寬逐漸減小，且主頻與帶寬的比值隨機改變。圖7a為隨機反射系數，圖7b為合成記錄，從上往下信號的分辨率降低，圖7c為匹配追蹤結果，圖7d為相對誤差百分比。為作對比，對信號2同時進行基于Ricker子波的匹配追蹤(算法的終止條件不變)，時頻譜如圖8所示。其中，圖8a為圖7b基于C子波的匹配分解結果，圖8b為基于Ricker子波的匹配分解結果。對比發(fā)現，兩個分解結果存在一定的差異：C子波匹配追蹤譜主頻隨時間呈近似下降的趨勢，符合模型設計，而Ricker子波匹配譜在600ms后主頻呈遞增的趨勢，與實際不符;C子波匹配譜的強能量團與信號的峰值分布符合度較高，如信號2在680～740ms有兩個強峰值，C子波匹配譜上有2個強能量團，而Ricker子波匹配譜上則有3個，在800ms附近存在強峰值，C子波匹配譜上存在強能量顯示，Ricker子波匹配譜上能量不明顯。

圖5 信號1的分解結果a 原始信號組成示意; b 基于C子波的匹配追蹤時頻譜

圖6 信號1的短時Fourier變換(a)和S變換(b)時頻譜

圖7 信號2a 反射系數; b 合成記錄; c 匹配結果; d 相對誤差

圖8 信號2匹配追蹤時頻譜a 基于C子波; b 基于Ricker子波

對比結果表明，匹配追蹤的分解效果依賴于原子庫與信號基本結構的匹配度。C原子庫的參量調控靈活，波形較Ricker子波庫豐富，對實際地震信號基本結構(時變地震子波)的適應性更強。因此，對于子波特征值(主旁瓣極值比和寬度比)時變的地震信號，基于C子波的匹配效果優(yōu)于基于Ricker子波的匹配效果。

3 實際資料處理

圖9為陸上某探區(qū)的過井測線地震剖面。該區(qū)地下構造簡單，地層平緩，斷裂較少，淺表氣遍布全區(qū)，第四系沉積以砂泥巖互層為主，成巖性較差，特殊的地質條件造就了該區(qū)地震資料的特殊性。經勘探發(fā)現該井存在工業(yè)氣流，綜合研究后認為該地區(qū)氣藏為鼻狀構造背景上的巖性氣藏。根據以往經驗，含油氣地區(qū)一般會有明顯地震異常。從圖9 可以看出，在井附近同相軸的連續(xù)性變差，有下拉趨勢，主頻變低，在主測線和聯絡線上地震相外形呈塔狀。這種地震異常雖然存在，但并非特別明顯。針對該地區(qū)的氣藏特征，應用基于C子波的匹配追蹤算法，并對資料進行分頻處理，為該地區(qū)開展巖性氣藏研究提供依據。

圖9 某工區(qū)過井測線地震剖面

提取該過井剖面井旁道目的層段(1000～2000ms)地震記錄，分別用Ricker原子和C原子對該記錄進行匹配追蹤，結果如圖10所示。對比發(fā)現，兩種匹配結果都能實現地震信號的高精度匹配分解，較好地反映信號的時頻特征;在局部細節(jié)上，C原子匹配的時間分辨率優(yōu)于Ricker原子匹配，如：在1900ms附近(黑框標注部分)，C原子匹配較好地將20Hz主頻的兩個能量團區(qū)分開。因此，利用C原子對地震記錄進行匹配追蹤，更有助于薄層的精細刻畫。

圖11是對過井地震剖面進行分頻處理得到的單頻調諧體剖面，不同頻率剖面顯示的振幅相對強弱各不相同。通過頻率切片，可以清晰地觀察到各層位能量不同的展布和變換趨勢。圖11中黑色橫線為氣層所處深度(1650ms附近)，觀察井(紅線)附近目的層的能量，在12Hz左右的分頻剖面上達到最大，之后隨著頻率的增加而降低，這是因為薄層的調諧頻率在12Hz左右。在紅線與黑線相交位置附近，右側的能量在12Hz之后逐漸減小，而左側在25Hz剖面上仍有較強的能量，這是由于右側含氣層的影響。匹配追蹤分解結果為目的層油氣藏的預測和描述提供了可靠的依據。

圖10 過井地震剖面目的層段井旁單道匹配結果a 地震記錄; b Ricker原子匹配時頻譜; c C原子匹配時頻譜

圖11 過井地震剖面分頻處理得到的單頻調諧體剖面

4 結束語

Ricker類子波是一類新型的理論子波，波譜主頻以及高低頻的衰減速度由三參量聯合控制。根據實際需要，將該類子波簡化成二參量C子波，波譜主頻、帶寬由二參量獨立控制。針對子波波形，提出以主旁瓣極值比PR和主旁瓣寬度比WR來定量表征，并闡述了C子波波形比Ricker子波、寬帶Ricker子波更豐富，更能匹配時變地震子波。

根據C子波建立原子庫，用于信號的匹配追蹤算法。模型試算表明，與短時Fourier變換和S變換相比，匹配追蹤得到的時頻譜圖具有較好的時頻聚焦性;對于時變地震記錄，C子波與信號基本結構的匹配度較Ricker子波高，匹配追蹤效果更好。最后，將基于C子波的匹配追蹤應用于實際資料的分頻處理，取得了較好的解釋效果。

但同時需要指出的是，基于C子波的匹配追蹤算法的運算量較大，運行速度的提升是需要進一步改進優(yōu)化的方向。

參 考 文 獻

[1] 熊曉軍，賀振華，黃德濟，等．廣義S變換在地震高分辨處理中的應用[J]．勘探地球物理進展，2006，29(6)：415-418

Xiong X J,He Z H,Huang D J,et al.Application of generalized S transform in seismic high resolution processing[J].Progress in Exploration Geophysics,2006,29(6):415-418

[2] 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

[3] Sinha S,Routh P S,Anno P D,et al.Spectral decomposition of seismic data with continuous-wavelet transforms[J].Geophysics,2005,70(6):19-25

[4] Mallat S,Zhang Z.Matching pursuits with time-frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415

[5] Qian S,Chen D.Signal representation via adaptive normalized Gaussian function[J].Signal Processing,1994,36(1):1-11

[6] Liu J,Wu Y,Han D,et al.Time-frequency decomposition based on Ricker wavelet[J].Expanded Abstracts of 74thSEG Annual Internat Mtg,2004,1937-1940

[7] Liu J,Marfurt K J.Matching pursuit decomposition using Morlet wavelets[J].Expanded Abstracts of 75thSEG Annual Internat Mtg,2005,786-789

[8] 宋新武，鄭浚茂，范興燕，等．基于Ricker子波匹配追蹤算法在薄互層砂體儲層預測中的應用[J]．吉林大學學報(地球科學版)，2011，41(增刊1)：387-392

Song X W,Zheng J M,Fan X Y,et al.Application of the thin-interbedded reservoir prediction based on Ricker wavelet match tracing algorithm[J].Journal of Jilin University(Earth Science Edition),2011,41(S1):387-392

[9] 梅金順，王潤秋．Ricker類地震子波[J]．石油地球物理勘探，2012，47(增刊1)：8-14

Mei J S,Wang R Q.Ricker-like seismic wavelets[J].Oil Geophysical Prospecting,2012,47(S1):8-14

[10] Wang Y.Seismic time-frequency spectral decomposition by matching pursuit[J].Geophysics,2007,72(1):13-20

[11] Barnes A E.Instantaneous frequency and amplitude at the envelope peak of a constant-phase wavelet[J].Geophysics,1991,56(7):1058-1060

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

Zhang F C,Li C H.Matching pursuit time-frequency analysis of non-stationary seismic signals[J].Geophysical and Geochemical Exploration,2011,35(4):546-552

[13] 張繁昌,李傳輝.匹配追蹤保幅地震AVF剖面及同頻率剖面的構建[J].石油物探,2013,52(3):234-239

Zhang F C,Li C H.Amplitude-preserved AVF and isofrequency section construction based on seismic matching pursuit[J].Geophysical Prospecting for Petroleum,2013,52(3):234-239

[14] 李傳輝,張繁昌.地震信號可變分辨率匹配追蹤頻譜成像方法[J].石油物探,2012,51(3):213-218

Li C H,Zhang F C.Variable resolution matching pursuit spectrum imaging of seismic signals[J].Geophysical Prospecting for Petroleum,2012,51(3):213-218

[15] 張顯文,韓立國,王宇,等.地震信號譜分解匹配追蹤快速算法及其應用[J].石油物探,2010,49(1):1-6

Zhang X W,Han L G,Wang Y,et al.Seismic spectral decomposition fast matching persuit algorithm and its application[J].Geophysical Prospecting for Petroleum,2010,49(1):1-6

[16] 陳發(fā)宇,楊長春.對依據頻率的匹配追蹤快速算法的改進[J].石油物探,2009,48(1):80-83

Chen F Y,Yang C C.An improvement for frequency-domainated fast matching pursuits algorithm[J].Geophysical Prospecting for Petroleum,2009,48(1):80-83

[17] 趙天姿,宋煒,王尚旭.基于匹配追蹤算法的時頻濾波去噪方法[J].石油物探,2008,47(4):367-371

Zhao T Z,Song W,Wang S X.Time-frequency filtering de-noise method based on matching pursuit algorithm[J].Geophysical Prospecting for Petroleum,2008,47(4):367-371