基于Ricker類(lèi)地震子波的匹配追蹤

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

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

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

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

Ricker類(lèi)地震子波是由Ricker子波演變得到的零相位理論子波，具有3個(gè)控制參量，旁瓣能量的收斂速度可變，延遲時(shí)間可控，子波庫(kù)波形更豐富，與地震子波的匹配較為靈活[9]。我們?cè)谘芯孔硬úㄐ坞S3個(gè)控制參量變化關(guān)系的基礎(chǔ)上，對(duì)該子波作適當(dāng)?shù)暮?jiǎn)化，建立參量與波形特征之間的半定量關(guān)系。最后，以簡(jiǎn)化的子波構(gòu)建原子庫(kù)，通過(guò)三步法原理[10]實(shí)現(xiàn)信號(hào)的匹配追蹤分解。

1 基本原理

1.1 Ricker類(lèi)地震子波

Ricker類(lèi)地震子波頻域解析式為[9]

(1)

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

(2)

這里記C(f)為二參量Ricker類(lèi)地震子波，簡(jiǎn)稱C子波。其中，fm和c為待定參量;A為歸一化系數(shù)。通過(guò)簡(jiǎn)化和改寫(xiě)，將子波譜表述成Ricker譜(主頻為fm)的c次方(c=1時(shí)對(duì)應(yīng)Ricker子波)。fm為波譜主頻，控制波譜的等比例伸縮;c控制波譜主頻兩側(cè)能量的衰減速度。

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

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

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

圖1 對(duì)稱C子波示意圖解

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

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

1.2 基于C子波的匹配追蹤

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

(3)

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

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

(4)

(5)

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

(6)

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

(7)

1.3 時(shí)頻表征

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

(8)

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

圖3是單個(gè)C原子(主頻fm取50Hz，c取0.7)的時(shí)頻譜。該時(shí)頻表征法具有較高的時(shí)頻分辨率，能同時(shí)獲取振幅和相位等信息的時(shí)頻分布，并且避免了窗口截?cái)嘈?yīng)。

圖3 C原子時(shí)頻譜

2 模型試算

設(shè)計(jì)兩個(gè)合成信號(hào)進(jìn)行基于C子波庫(kù)的匹配追蹤算法測(cè)試，時(shí)間采樣率均為1ms。

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

圖4 信號(hào)1重構(gòu)圖解a 原始信號(hào); b 重構(gòu)信號(hào); c相對(duì)誤差

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

圖5 信號(hào)1的分解結(jié)果a 原始信號(hào)組成示意; b 基于C子波的匹配追蹤時(shí)頻譜

圖6 信號(hào)1的短時(shí)Fourier變換(a)和S變換(b)時(shí)頻譜

圖7 信號(hào)2a 反射系數(shù); b 合成記錄; c 匹配結(jié)果; d 相對(duì)誤差

圖8 信號(hào)2匹配追蹤時(shí)頻譜a 基于C子波; b 基于Ricker子波

對(duì)比結(jié)果表明，匹配追蹤的分解效果依賴于原子庫(kù)與信號(hào)基本結(jié)構(gòu)的匹配度。C原子庫(kù)的參量調(diào)控靈活，波形較Ricker子波庫(kù)豐富，對(duì)實(shí)際地震信號(hào)基本結(jié)構(gòu)(時(shí)變地震子波)的適應(yīng)性更強(qiáng)。因此，對(duì)于子波特征值(主旁瓣極值比和寬度比)時(shí)變的地震信號(hào)，基于C子波的匹配效果優(yōu)于基于Ricker子波的匹配效果。

3 實(shí)際資料處理

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

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

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

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

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

圖11 過(guò)井地震剖面分頻處理得到的單頻調(diào)諧體剖面

4 結(jié)束語(yǔ)

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

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

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

參 考 文 獻(xiàn)

[1] 熊曉軍，賀振華，黃德濟(jì)，等．廣義S變換在地震高分辨處理中的應(yīng)用[J]．勘探地球物理進(jìn)展，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子波匹配追蹤算法在薄互層砂體儲(chǔ)層預(yù)測(cè)中的應(yīng)用[J]．吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，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] 梅金順，王潤(rùn)秋．Ricker類(lèi)地震子波[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)地震信號(hào)匹配追蹤時(shí)頻分析[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剖面及同頻率剖面的構(gòu)建[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] 李傳輝,張繁昌.地震信號(hào)可變分辨率匹配追蹤頻譜成像方法[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] 張顯文,韓立國(guó),王宇,等.地震信號(hào)譜分解匹配追蹤快速算法及其應(yīng)用[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ā)宇,楊長(zhǎng)春.對(duì)依據(jù)頻率的匹配追蹤快速算法的改進(jìn)[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] 趙天姿,宋煒,王尚旭.基于匹配追蹤算法的時(shí)頻濾波去噪方法[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