稀疏反射系數(shù)頻率域正余弦分量協(xié)同反演方法

張繁昌 何晉越 桑凱恒* 秦廣勝 張佳佳

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②中國石化中原油田分公司，河南濮陽 457099)

1 引言

隨著油氣資源勘探程度的不斷加深，對勘探精度的要求也越來越高，即對地震資料的分辨率提出了越來越高的要求。由于受調(diào)諧效應(yīng)的影響，應(yīng)用常規(guī)地震資料及傳統(tǒng)處理方法無法識別厚度小于調(diào)諧厚度的薄層和小尺度地質(zhì)構(gòu)造，難以提供高精度的解釋成果。識別調(diào)諧厚度以下的信息，突破常規(guī)地震分辨率的極限是近年來地震領(lǐng)域資料處理研究的重點之一[1]。Portniaguine等[2，3]提出了譜反演方法，它不依賴測井資料等先驗信息，只需利用地震數(shù)據(jù)本身的頻譜信息就能反演地層稀疏反射系數(shù)； Chopra等[4，5]對譜反演方法進(jìn)行了詳細(xì)的論證、推導(dǎo)； Castagna等[6]提出奇偶分解理論； Puryear等[7]對奇偶分解理論進(jìn)行詳細(xì)闡述，把稀疏脈沖反射系數(shù)分解成奇、偶脈沖分量，進(jìn)一步提高了譜反演對小于調(diào)諧厚度的薄層的分辨能力； Yuan等[8]給出了一種奇偶分解譜反演的等價矩陣形式； 劉萬金等[9]、陳祖慶等[10]、張華等[11]及劉潔等[12]應(yīng)用譜反演理論并結(jié)合不同反演算法對實際地震資料進(jìn)行處理，均取得不錯的效果。

在時間域?qū)Ψ瓷湎禂?shù)進(jìn)行奇偶分解，需給定奇偶脈沖分量的時間厚度，該參數(shù)對反演結(jié)果的影響很大且不易確定。研究發(fā)現(xiàn)，在時間域?qū)⒎瓷湎禂?shù)序列分解成奇、偶脈沖分量，與在頻率域?qū)⒎瓷湎禂?shù)序列的頻譜分解成正、余弦分量是等價的。基于此等價關(guān)系，本文提出一種直接在頻率域利用地震數(shù)據(jù)譜的正、余弦分量進(jìn)行協(xié)同反演的方法，推導(dǎo)出反演目標(biāo)函數(shù)。對反演目標(biāo)函數(shù)的求解有多種方法[13-16]，為了獲得稀疏反射系數(shù)序列，本文根據(jù)壓縮感知理論，采用梯度投影稀疏重構(gòu)(Gradient Projection for Sparse Reconstruction，GPSR)基追蹤算法[17，18]對目標(biāo)函數(shù)進(jìn)行求解。通過一維稀疏脈沖反射系數(shù)模型和實際地震資料反演，對頻率域正、余弦分量協(xié)同反演結(jié)果與稀疏脈沖反演結(jié)果進(jìn)行對比，證明頻率域正、余弦分量協(xié)同反演方法能取得高分辨率反射系數(shù)序列，有助于識別原始地震剖面上無法分辨的薄層，且反演剖面的信噪比得到提高。

2 方法原理

2.1 正、余弦分量協(xié)同反演基本原理

正弦和余弦信號的傅里葉變換為

F[sin(ω0t)]=jπ[δ(ω+ω0)-δ(ω-ω0)]

(1)

F[cos(ω0t)]=π[δ(ω+ω0)+δ(ω-ω0)]

(2)

由傅里葉變換的對稱性質(zhì)，即已知F(ω)=F[f(t)]，則F[F(t)]=2πf(-ω)，可得到

F[δ(t+t0)-δ(t-t0)]=2jsin(t0ω)

(3)

F[δ(t+t0)+δ(t-t0)]=2cos(t0ω)

(4)

其對應(yīng)波形及頻譜如圖1所示。

由此可知，時間間隔為2t0的奇脈沖分量在頻率域表現(xiàn)為以t0為周期的正弦分量，時間間隔為2t0的偶脈沖分量在頻率域表現(xiàn)為以t0為周期的余弦分量。根據(jù)Puryear的奇偶分解理論，假設(shè)反射系數(shù)是稀疏脈沖的，任意反射系數(shù)序列都可以分解成奇、偶脈沖分量，因此，對應(yīng)于頻率域，任意反射系數(shù)的頻譜可以分解成正、余弦分量。 因此將反射系數(shù)頻譜分解成正、余弦分量，構(gòu)建用正、余弦分量表示的目標(biāo)函數(shù)，在頻率域進(jìn)行正、余弦分量協(xié)同反演。

圖1 奇脈沖分量(a)、偶脈沖分量(c)及其對應(yīng)頻譜(b，d)

反射系數(shù)序列r(t)可表示為

(5)

式中：rn和tn分別為第n個采樣點上的反射系數(shù)值和對應(yīng)的時間位置；N為采樣點個數(shù)。

根據(jù)褶積模型[19]，地震信號s(t)可表示為

(6)

式中：w(t)為地震子波； “*”表示褶積運算。對式(6)做傅里葉變換，即得褶積模型在頻率域表達(dá)式

(7)

式中：f為頻率；S(f)、W(f)和R(f)分別為地震信號s(t)、地震子波w(t)和反射系數(shù)序列r(t)的頻譜。因此反射系數(shù)序列r(t)的頻譜可表示為

(8)

對式(8)分別取實部和虛部得到

(9)

(10)

式中aWmax為白噪因子，可防止分母出現(xiàn)零值。其中：Wmax為地震子波頻譜W(f)的最大振幅值；a是一個很小的系數(shù)(一般取0.01)。將式(4)、式(3)分別代入式(9)、式(10)，得到

(11)

(12)

由此可見，反射系數(shù)頻譜實部是不同系數(shù)、不同頻率的余弦分量的疊加，這里的系數(shù)就是各采樣點的反射系數(shù)值rn，頻率即為該反射系數(shù)值對應(yīng)的時間位置tn，同時，這一系列余弦分量對應(yīng)于時間域一系列偶脈沖分量。同理，反射系數(shù)頻譜虛部是正弦分量的疊加，在時間域?qū)?yīng)于一系列奇脈沖分量。

地震信號頻帶寬度都是有限的，即地震信號都是帶限信號。以M區(qū)一道地震記錄為例，其振幅譜如圖2所示，可知其只在5～85Hz范圍內(nèi)有明顯的振幅值(紅框內(nèi)部分)。設(shè)選取的地震信號頻率范圍是[fmin，fmax]，頻率采樣間隔為Δf，有

fm=(m-1)×Δf+f1m=1,2,…,M

(13)

圖2 M區(qū)一道地震記錄的振幅譜

在[fmin，fmax]內(nèi)對式(9)、式(10)進(jìn)行離散，建立反射系數(shù)離散譜方程

(14)

將式(14)改寫成矩陣形式

(15)

其中

Cmn=cos(2πtnfm)

Smn=sin(2πtnfm)

r=[r1,r2,…,rn]T

式中：C即余弦分量矩陣；S即正弦分量矩陣；r為反射系數(shù)序列組成的向量。式(15)即正、余弦分量協(xié)同反演目標(biāo)函數(shù)表達(dá)式。

2.2 GPSR算法

為了求解稀疏的反射系數(shù)序列，應(yīng)使用稀疏約束的基追蹤算法對目標(biāo)函數(shù)進(jìn)行求解。目前有大量文獻(xiàn)提供了各種基追蹤求解算法[20]，本文選用梯度投影稀疏重構(gòu)(GPSR)算法對式(15)進(jìn)行求解。

將目標(biāo)函數(shù)式(式(15))改寫為基追蹤的標(biāo)準(zhǔn)形式

(16)

將待反演的反射系數(shù)序列r分解為正數(shù)序列u和負(fù)數(shù)序列-v，則r可寫成r=u-v(u≥0,v>0)，因此r的L1范數(shù)可寫成

(17)

式中l(wèi)n=[1,1,…,1]1×n。將式(17)代入式(16)，則求解式(16)轉(zhuǎn)為求解目標(biāo)函數(shù)

(18)

進(jìn)一步進(jìn)行下列變換

(19)

即得到新的目標(biāo)函數(shù)表達(dá)形式

(20)

通過以下步驟的迭代，最終求得式(20)的解。

(1)初始化， 給定z(0)，λ=0.01， 選擇參數(shù)β∈(0,1)，μ∈(0,0.5)，令k=0；

(2)利用式(21)和式(22)計算α0， 并用mid(αmin,α0,αmax) 代替α0；

(21)

(22)

(23)

(4)進(jìn)行收斂判定，并在結(jié)果滿足精度要求時終止迭代;如果不滿足精度要求，則令k=k+1，返回步驟(2)。

3 模型和實際數(shù)據(jù)測試

3.1 理論模型試驗

為了體現(xiàn)正余弦分量協(xié)同反演的優(yōu)勢，將時域褶積模型直接作為目標(biāo)函數(shù)，并采用業(yè)內(nèi)常用的稀疏脈沖反演方法[24，25]進(jìn)行反演，將其結(jié)果作為對比。稀疏脈沖反演目標(biāo)函數(shù)為

(24)

式中：W為子波矩陣；s為觀測得到的地震信號組成的向量。

本模型的稀疏脈沖反演結(jié)果如圖3e所示，它與25Hz的Ricker子波合成的地震記錄如圖3f所示。觀察圖3d和圖3f可發(fā)現(xiàn)，兩種方法得到的反射系數(shù)恢復(fù)出的合成地震記錄與模型正演地震記錄匹配較好，但對比圖3c與圖3e可看到，基于GPSR的正余弦分量協(xié)同反演結(jié)果與模型匹配得很好， 反射系數(shù)的位置、幅值都能較為精確地匹配，而稀疏脈沖反演結(jié)果中，大部分反射系數(shù)值恢復(fù)不到位，第2、6組不能將兩個脈沖信號完全分開，而第4、8組反演結(jié)果時間間隔偏大，且在相反極性對兩側(cè)出現(xiàn)模型不存在的脈沖。試驗說明，基于GPSR的正余弦分量協(xié)同反演能夠識別時間間隔小于分辨率極限的反射界面，且反演結(jié)果穩(wěn)定而精準(zhǔn)。

圖3 一維理論模型試驗結(jié)果

3.2 實際資料應(yīng)用

為了進(jìn)一步分析正余弦分量協(xié)同反演的實際應(yīng)用效果，將其應(yīng)用于中國東部E區(qū)一過井測線的實際地震資料(圖4a)。基于GPSR的正余弦分量協(xié)同反演過程中采用的地震子波由譜模擬法[26]估計得到，同時根據(jù)地震頻譜選取5～70Hz為地震有效頻帶，反演獲得的反射系數(shù)剖面如圖4b所示。作為對比，稀疏脈沖反演獲得的反射系數(shù)剖面如圖4c所示。 對比圖4a與圖4b可見，相比原始地震剖面基于GPSR的正余弦分量協(xié)同反演獲得的反射系數(shù)剖面明顯具有更高的分辨率，能將由于相干形成的復(fù)波同相軸分離開來，識別原始地震剖面上無法分辨的、厚度小于調(diào)諧厚度的薄層和小尺度地質(zhì)體，反映地層真實的接觸關(guān)系，刻畫地層真實的尖滅位置，提供高分辨率的、更加豐富的地質(zhì)信息。對比圖4b與圖4c可見，基于GPSR的正余弦分量協(xié)同反演結(jié)果相較于稀疏脈沖反演結(jié)果在保持良好信噪比的同時分辨率有明顯提升，能識別出更小厚度的薄層，且橫向連續(xù)性更好，有助于薄互層、小尺度地質(zhì)體的精細(xì)識別。觀察三張圖黑色虛線框內(nèi)的部分，圖4b的基于GPSR的正余弦分量協(xié)同反演結(jié)果能精細(xì)刻畫小尺度地質(zhì)體和薄層以及它們之間的接觸關(guān)系，而在圖4c的稀疏脈沖反演結(jié)果中薄層反射軸無法分開，有效信息被淹沒，不能反映地質(zhì)構(gòu)造的真實樣式。

圖4 實際資料應(yīng)用結(jié)果

4 結(jié)束語

研究發(fā)現(xiàn)，時間域的奇偶脈沖分量與頻率域的正余弦分量是彼此等價的。在時間域?qū)Ψ瓷湎禂?shù)進(jìn)行奇偶分解，需要給定奇偶脈沖分量的時間厚度，該參數(shù)對反演結(jié)果影響很大而且不易確定。而以此等價關(guān)系為基礎(chǔ)，直接在頻率域利用地震數(shù)據(jù)譜的正余弦分量進(jìn)行協(xié)同反演可避免此問題。從反演算法看出，本文反演方法不依賴于測井資料等先驗信息，只需要地震數(shù)據(jù)本身的頻譜信息就能反演得到地層稀疏反射系數(shù)。

從模型測試和實際數(shù)據(jù)應(yīng)用結(jié)果看，相比于時域褶積模型反演結(jié)果，正余弦分量協(xié)同反演在保持良好信噪比的同時，分辨率明顯提升，能夠識別原始地震剖面上分辨不出的、小于調(diào)諧厚度的薄層和小尺度地質(zhì)體，更有利于地震沉積學(xué)研究。