周東勇, 文曉濤, 賀振華, 蒲 勇, 黃德濟, 胡軍輝

(1.油氣藏地質(zhì)及開發(fā)工程國家重點實驗室(成都理工大學(xué)),成都 610059;2.中國石化股份有限公司 勘探南方分公司,成都 610041)

波阻抗反演開始出現(xiàn)于20世紀70年代初，是一種以研究巖性油氣藏為主的地震勘探技術(shù)。目前，波阻抗反演的方法有很多，如模擬退火法、遺傳算法等等。本文研究采用匹配追蹤算法對地震數(shù)據(jù)進行分解，得相對反射系數(shù)，進而依據(jù)波阻抗公式進行波阻抗反演。匹配追蹤算法(Matching Pursuits，本文統(tǒng)一簡寫為MP算法)由Mallat和 Zhang于1993年首先提出[1]，之后很快被應(yīng)用于地震信號處理領(lǐng)域。2003年Castagna等將MP算法用于地震譜分解[2]; 2007年Y.H.Wang將該算法應(yīng)用于時頻分析[3]；2012年張繁昌等將該算法用于瞬時譜識別來確定三角洲砂體尖滅線[4]等等。雖然應(yīng)用MP算法對信號稀疏分解時，計算量十分大，但該算法有很多優(yōu)良的性質(zhì),如對信號自適應(yīng)的靈活表達，較高的時間、頻率分辨率,暫態(tài)結(jié)構(gòu)的局部自適應(yīng)性,信號結(jié)構(gòu)的參數(shù)表示更加靈活等[5]。

1 地震信號的稀疏分解

信號表示實際上就是把給定的信號在已知的函數(shù)集(函數(shù)集中的任意一個函數(shù)稱為基函數(shù))上進行分解, 然后在變換域上表達原始信號或原始信號的主要成分。若在變換域上用盡量少的基函數(shù)來表示原始信號或原始信號的主要成分, 就是信號的稀疏表示, 得到信號稀疏表示的過程就是信號的稀疏分解[6]。若用這種分解過程分解地震信號，則是地震信號的稀疏分解。MP算法就是實現(xiàn)地震信號稀疏分解的算法之一。

設(shè)研究的地震信號為f,信號長度為N,且f∈H，H為希爾伯特空間。MP算法將原始地震信號f稀疏表示為

(1)

式中:gγi為函數(shù)集中的基函數(shù)(在MP算法中函數(shù)集被稱為過完備原子庫，基函數(shù)被稱為原子);ci為地震信號f在基函數(shù)gγi上的投影分量[7]。式(1)和條件n?N集中體現(xiàn)了地震信號稀疏表示的思想。地震信號在過完備原子庫上的分解結(jié)果一定是稀疏的[1]。

2 過完備原子庫的建立

如何建立合理的過完備原子庫對能否精確匹配地震信號至關(guān)重要[8]。雷克子波與地震子波的波形相似，適合分析地震信號的特性，在地震反演中得到了廣泛的應(yīng)用。此處，筆者提出利用非零相位雷克子波來構(gòu)建過完備原子庫的思想，實現(xiàn)對地震信號的稀疏分解。

2.1 零相位雷克子波

Ricker子波表達式如下[9]

r(t-μ)=[1-2(πfp(t-μ))2]e-(πfp(t-μ))2

(2)

式中:fp為峰值頻率[10]；μ為時間延遲；e為自然對數(shù)。其波形是由1個主瓣2個旁瓣組成，圖1為30 Hz、50 Hz雷克子波波形圖。

圖1 30 Hz、50 Hz雷克子波波形圖Fig.1 The shape of Ricker wavelets for 30 Hz and 50 Hz

2.2 非零相位雷克子波

本文采用的非零相位雷克子波，由解析信號推導(dǎo)而來[11]。以下是Ricker子波經(jīng)ω相位旋轉(zhuǎn)后的表達式

r′(t)=r(t)cosω-r*(t)sinω

(3)

其中:r(t)=[1-2(πfpt)2]e-(πfpt)2;r*(t)為r(t)的Hilbert變換。

對Ricker子波分別進行30°、90°、150°和180°相位旋轉(zhuǎn)后如圖2所示。

2.3 過完備原子庫

過完備原子庫D={gγ}γ∈Γ由非零相位雷克子波公式經(jīng)離散化、歸一化形成，為敘述方便現(xiàn)將非零相位雷克子波公式寫成

gγ∈Γ(t-u)=r(t-u)cosω-r*(t-u)sinω

(4)

式中:r(t)為雷克子波；r*(t)為雷克子波的希爾伯特變換；γ=(fp,ω,u)是時頻參數(shù)；Γ是所有時頻參數(shù)的集合。fp為主頻參數(shù)；ω為相位參數(shù)；u為位移參數(shù)。時頻參數(shù)可以按下式離散化

γ=(i·dfp,j·dω,n·du)

(5)

然后通過式‖gγ(t)‖=1將每個原子歸一化。

在過完備原子庫中，1個原子gγ(t)由參數(shù)(fp,ω,u)決定。參數(shù)(fp,ω)決定原子gγ(t)的波形，而參數(shù)u只決定原子gγ(t)出現(xiàn)的位置。

3 基于MP算法的波阻抗反演流程

3.1 MP算法原理

設(shè)f為待分解的地震信號，D為用于分解地震信號的過完備原子庫，gγ為過完備原子庫D中的原子，gγ∈D, ‖gγ‖=1。

MP算法分解地震信號f過程如下。

(1)從過完備原子庫D中選擇與地震信號f最為匹配的原子(即為最佳原子)gγ0，使其滿足以下條件

|〈f,gγ0〉|=sup|〈f,gγ〉|

(6)

其中，〈f,gγ〉錯誤！未找到引用源。f與原子庫中任意原子gγ的內(nèi)積，sup在泛函分析中稱為上確界，此處即為最大值的意思。

此時可將地震信號f分解為

f=〈f,gγ0〉gγ0+R1f

(7)

圖2 30°、90°、150°和180°相位雷克子波波形圖Fig.2 The shape of Ricker wavelets for 30°, 90°, 150° and 180°

其中R1f為地震信號f與最佳原子gγ0匹配后的剩余地震信號(即為殘差)。

由于gγ0與R1f正交，所以其能量表示為

(8)

(2)將剩余地震信號(或殘差)看做新地震信號，重復(fù)以上過程，不斷將新地震信號匹配到過完備原子庫上，經(jīng)過k+1次匹配可以得到

Rkf=〈Rkf,gγk〉gγk+Rk+1f

(9)

其中g(shù)γk滿足

|〈Rkf,gγk〉|=sup|〈Rkf,gγ〉|

(10)

由于Rk+1f與grk正交，則有

(11)

(3)經(jīng)過n次分解后最終把地震信號f分解為

(12)

能量分解為

(13)

此處，令f=R0f,Rnf為原地震信號分解為n個最佳匹配原子的線性組合后所產(chǎn)生的誤差。

(4)分解完成的判定。MP算法是一個不斷迭代的過程，需要一定的標(biāo)準(zhǔn)來判定分解是否完成。筆者認為，通過設(shè)定殘差Rnf的能量閾值，來作為分解結(jié)束條件的方法是較好的判定標(biāo)準(zhǔn)。即當(dāng)能量值滿足‖Rnf‖2<ε(N)時，分解停止；否則繼續(xù)運行直到滿足條件。其中ε(N)為與采樣點N相關(guān)的閾值，可根據(jù)經(jīng)驗自行設(shè)定。

3.2 波阻抗反演流程

為更好地描述波阻抗反演過程，設(shè)計如圖3所示流程圖。整個流程大體可分為兩部分：首先通過MP算法不斷匹配地震信號或剩余地震信號，獲取相對反射系數(shù)及其對應(yīng)位置；第二依據(jù)波阻抗遞推公式(14)計算波阻抗。

(14)

圖3 基于MP算法地震波阻抗反演流程流程圖Fig.3 The workflow of seismic wave impedance inversion based on MP algorithm

圖中，因?qū)嶋H地震數(shù)據(jù)中的反射系數(shù)一般在±0.3之間，故將相對反射系數(shù)sk約束到±0.3之間，得近似反射系數(shù)Sk；初始波阻抗Z0由模型給出。

4 基于雷克子波原子庫的MP算法用于地震波阻抗反演的可行性

為了驗證算法的有效性，筆者用合成單道地震信號和楔形地震正演模型2種方法分別進行了MP算法反演測試。其中，前者主要是驗證MP算法用于地震信號分解的可行性，測試其對地震信號分解的精度；后者主要是簡單模擬實際地震數(shù)據(jù)，驗證其對地震波阻抗反演的宏觀效果。

4.1 MP算法用于地震信號稀疏分解的可行性及其分解精度分析

為驗證其對地震信號稀疏分解的可行性及分解精度，設(shè)計了3種合成地震信號。3種合成地震信號除第二層反射系數(shù)的位置不同外，其他均采用相同參數(shù)的雷克子波(此處雷克子波長度λ為40個采樣點，2 ms采樣，其參數(shù)如表1所示)合成長度為N(此處N為128個采樣點，2 ms采樣)的地震信號。合成地震信號第二層反射系數(shù)的位置分別取在第90、60、50個采樣點上，分別依次為兩層間無調(diào)諧、λ/2調(diào)諧、λ/4調(diào)諧3種情況。

表1 合成地震信號所需雷克子波參數(shù)Table 1 The Ricker wavelet parameter for the synthetic seismic signal

4.1.1 兩層間無調(diào)諧(兩層間厚度大于λ)時

對比圖(圖4-C和圖4-A對比、圖4-D和圖4-B對比)和列表(表2和表1相應(yīng)參數(shù)對比)可得出：兩層間無調(diào)諧時，層位、主頻、相位及系數(shù)的相對大小都能被準(zhǔn)確反演出；重構(gòu)地震信號和合成地震信號幾乎完全相同；殘差基本為零。

4.1.2 兩層間λ/2調(diào)諧(層間厚度為λ/2)時

表2 MP算法分解兩層間無調(diào)解的地震信號所得雷克子波參數(shù)Table 2 The Ricker wavelet parameters from the synthetic seismic signal

對比圖(圖5-C和圖5-A對比、圖5-D和圖5-B對比)和列表(表3和表1相應(yīng)參數(shù)對比)可得出：兩層間λ/2調(diào)諧時，反演出的層位、主頻及系數(shù)的相對大小與真實值相比有偏差，但偏差不大；反演出的相位與真實值相比偏差較大；重構(gòu)地震信號和合成地震信號基本相同；殘差趨近于零。

4.1.3 兩層間λ/4調(diào)諧(層間厚度為λ/4)時

表3 MP算法分解兩層間λ/2調(diào)解的地震信號所得雷克子波參數(shù)Table 3 Ricker wavelet parameters from that MP decomposes seismic signal when λ/2 tunes

圖4 層間無調(diào)解時合成地震信號及反演結(jié)果Fig.4 Synthetic seismic signal & inversion results without the tuning between layers(A)反射系數(shù); (B)合成地震信號; (C)相對反射系數(shù); (D)重構(gòu)地震信號; (E)殘差

圖5 層間λ/2調(diào)解時合成地震信號及反演結(jié)果Fig.5 Synthetic seismic signal and the inversion result when λ/2 tunes(A)反射系數(shù); (B)合成地震信號; (C)相對反射系數(shù); (D)重構(gòu)地震信號; (E)殘差

圖6 層間λ/4調(diào)解時合成地震信號及反演結(jié)果Fig.6 Synthetic seismic signal and the inversion result when λ/4 tunes(A)反射系數(shù); (B)合成地震信號; (C)相對反射系數(shù); (D)重構(gòu)地震信號; (E)殘差

對比圖(圖6-C和圖6-A對比、圖6-D和圖6-B對比)和列表(表4和表1相應(yīng)參數(shù)對比)可得出：兩層間λ/4調(diào)諧時，反演出的主頻、相位與真實值相比偏差較大，層位和系數(shù)相對大小分別與真實值相比都有偏差，但偏差不大，仍可反映層的大體位置和系數(shù)的相對大小關(guān)系。重構(gòu)地震信號和合成地震信號也基本相同，殘差趨近于零。

表4 MP算法分解兩層間λ/4調(diào)解的地震信號所得雷克子波參數(shù)Table 4 Ricker wavelet parameters from that MP decomposes seismic signal when λ/4 tunes

通過以上對比分析可知，MP算法可用于地震信號的稀疏分解；在精度方面，兩層間λ/4調(diào)諧時，仍可大體反映層位和系數(shù)的相對大小關(guān)系。

4.2 MP算法用于楔形地震正演模型波阻抗反演分析

本節(jié)主要通過反演楔形正演模型來驗證該算法對實際地震數(shù)據(jù)的適用性，同時可觀察其對模型的波阻抗反演宏觀效果。為不失一般性，本節(jié)所用楔形正演模型(圖7-C)是由任意頻率、任意相位的雷克子波合成。

通過反演結(jié)果與模型對比(圖7中，(E)與(B)對比、(D)與(A)對比)可以得出：兩層間>λ/4時，可基本完全反演出模型中的信息；兩層間≤λ/4時，與模型相比雖略有些偏差，但仍可反映楔形模型的楔形趨勢。

5 實際地震數(shù)據(jù)波阻抗反演結(jié)果分析

該研究區(qū)域為塔河YY區(qū)，該區(qū)三疊系構(gòu)造較平緩，主要受NE向斷裂帶和鹽體的塑性活動控制，沿NE方向延伸形成大型圈閉群。三疊系發(fā)育的儲層主要有上、中、下3個油組的砂體。本文針對中油組砂體進行研究。從圖8-B中可看出，橢圓區(qū)域可能為砂體儲層，與實鉆情況吻合，驗證了該方法的有效性。

6 結(jié) 論

本文將Mallat和 Zhang提出的MP算法應(yīng)用于地震波阻抗反演中，模型試算和實際地震數(shù)據(jù)都對其有效性進行了驗證，并得出以下結(jié)論。

a.MP算法能夠?qū)崿F(xiàn)對地震信號的自適應(yīng)分解, 模型試算表明，其對地震信號分解的適用性，且分解精度相對較高。

b.對實際地震數(shù)據(jù)的運算與實鉆情況吻合，驗證了該方法的有效性，可為儲層預(yù)測提供有利依據(jù)。

c.地層厚度<λ/4時，MP算法計算出的相對反射系數(shù)誤差可能會偏大。因此，進一步提高該算法應(yīng)用于地震信號分解的精度，更精確地反演地下地質(zhì)體的波阻抗信息是筆者下一步的研究方向。

在本文成文過程中，作者與蔡涵鵬博士、高剛博士、楊小江同學(xué)及賈玉婷同學(xué)進行了有益的探討，在此向他們表示衷心的感謝。

圖7 地震楔形模型及MP算法反演成果圖Fig.7 Seismic model and the inversion result of MP algorithm(A)波阻抗模型; (B)反射系數(shù)模型; (C)合成地震記錄; (D)反演出的相對波阻抗; (E)反演出的相對反射系數(shù); (F)重構(gòu)地震信號

圖8 實際地震數(shù)據(jù)及其反演成果圖Fig.8 Actual seismic data and the result of inversion(A)塔河YY地區(qū)實際地震剖面圖; (B)實際地震數(shù)據(jù)波阻抗反演成果圖

[參考文獻]

[1] Mallat S, Zhang Z. Matching pursuits with time-frequency dictionaries [J]. IEEE Transactions on Signal Processing, 1993, 41: 2297-3415.

[2] Castagna J P, Sun S J, Roberb W S. Instantaneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons[J]. The Leading Edge, 2003, 22(2): 120-127.

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

[4] 張繁昌，李傳輝，印興耀．三角洲砂巖尖滅線的地震匹配追蹤瞬時譜識別方法[J]．石油地球物理勘探，2012，47(1)：82-88．

Zhang F C, Li C H, Yin X Y. Delta fringe line recognition based on seismic matching pursuit instantaneous spectral characteristics[J]. Oil Geophysical Prospecting, 2012, 47(1): 82-88. (In Chinese)

[5] 屈念念，劉江平，李家斌．Riker子波匹配追蹤算法及其改進[J]．工程地球物理學(xué)報，2009,6(6)：740-745．

Qu N N, Liu J P, Li J B. Revised matching pursuit algorithm using Ricker wavelets[J]. Chinese Journal of Engineering Geophysics, 2009, 6(6): 740-745. (In Chinese)

[6] 程文波，王華軍．信號稀疏表示的研究及應(yīng)用[J]．西南石油大學(xué)學(xué)報:自然科學(xué)版，2008，30(5)：148-151．

Cheng W B, Wang H J. The research and application of sparse signal representation[J]. Journal of Southwest Petroleum University (Science& Technology Edition), 2008, 30(5): 148-151. (In Chinese)

[7] 邵君，尹忠科，王建英，等．信號稀疏分解中過完備原子庫的集合劃分[J]．鐵道學(xué)報， 2006，28(1)：68-71．

Shao J, Yin Z K, Wang J Y,etal. Set partitioning of the over-complete dictionary in signal sparse decomposition[J]. Journal of the China Rall Way Society, 2006, 28(1): 68-71. (In Chinese)

[8] 黃捍東，郭飛，汪佳蓓，等．高精度地震時頻譜分解方法及應(yīng)用[J]．石油地球物理勘探，2012，47(5)：774-780．

Huang H D, Guo F, Wang J B,etal. The method and application of high precision seismic spectrum decomposition[J]. Oil Geophysical Prospecting, 2012, 47(5): 774-780. (In Chinese)

[9] Ricker N. The form of laws of propagation of seismic wavelets[J]. Geophysics, 1953, 18(1): 10-40．

[10] 云美厚，丁偉．地震子波頻率淺析[J]．石油物探，2005，44(6)：578-581．

Yun M H, Ding W. Analysis of seismic wavelet frequency[J]. Geophysical Prospecting for Petroleum, 2005, 44(6): 578-581. (In Chinese)

[11] 王純偉．MP算法在地震信號去噪中的應(yīng)用研究[D]．成都：西南交通大學(xué)檔案館，2010．

Wang C W. The Application Research of Matching Pursuit in Seismic Signal Denoising[D]. Chengdu: The Archive of Southwest Jiaotong University, 2010. (In Chinese)