基于改進(jìn)廣義S變換的復(fù)賽譜時(shí)頻域反褶積方法

劉春成, 周 慰, 牛 聰, 龍 丹, 巫南克

(1.中海油研究總院,北京 100027; 2.成都理工大學(xué)地球物理學(xué)院,四川成都 610059)

隨著地震勘探中對復(fù)雜構(gòu)造儲層精確描述的難度逐漸增大,運(yùn)用反褶積方法來實(shí)現(xiàn)高分辨處理已經(jīng)成為一種最常用的手段。Robinson等[1]假設(shè)地層反射系數(shù)滿足白譜特性,提出了預(yù)測反褶積模型,后來又補(bǔ)充引入Wiener濾波,得到了廣義的最小平方反褶積。為了擺脫這種假設(shè),Robert等[2]將地震記錄振幅譜轉(zhuǎn)換到復(fù)賽譜中,假定子波是低頻的,反射系數(shù)是高頻的,以一種非線性濾波的方式來提取子波。Wiggins[3]提出最小熵反褶積,通過求解線性算子,使道集中的尖脈沖特性得到加強(qiáng)來進(jìn)行反褶積。Milton等[4]提出一種混合相位反褶積方法。以上反褶積方法都是基于地震波始終保持不變的理想條件下,但是實(shí)際傳播過程中其會受到地下介質(zhì)的影響,產(chǎn)生能量的耗散效應(yīng),主頻會隨著傳播時(shí)間的增長而向低頻偏移,因此傳統(tǒng)的褶積模型不吻合于這個(gè)動(dòng)態(tài)過程。為了更加符合實(shí)際地震波傳播規(guī)律,同時(shí)也考慮到大地濾波效應(yīng),Clarke[5]最早提出了時(shí)域非平穩(wěn)反褶積方法。對于時(shí)頻域地震模型,Rosa[6]結(jié)合Ricker研究結(jié)論認(rèn)為，假設(shè)地震子波振幅譜是光滑的,可以平滑地震道將子波提取出來。唐博文、趙波等[7]對地震記錄求二次譜,提出了一種新的子波提取方法,實(shí)現(xiàn)了譜模擬反褶積。Margrave等[8-9]將地震記錄轉(zhuǎn)換到Gabor域,假設(shè)地下反射系數(shù)滿足白譜特性,通過對每個(gè)時(shí)間點(diǎn)振幅譜平滑,直接估計(jì)出衰減子波振幅譜,實(shí)現(xiàn)了時(shí)頻域反褶積。周懷來等[10]在S域中實(shí)現(xiàn)動(dòng)態(tài)反褶積,彌補(bǔ)了Gabor變換時(shí)窗固定的缺陷,但是沒有從根本上擺脫反射系數(shù)滿足白譜特性的假設(shè)。筆者通過對比Gabor變換、S變換[11]、改進(jìn)廣義S變換這3種時(shí)頻分析方法,證明改進(jìn)廣義S變換在時(shí)頻分析過程中的優(yōu)越性;通過調(diào)節(jié)改進(jìn)廣義S變換得到參數(shù)獲得衰減地震道的時(shí)頻譜,將其轉(zhuǎn)換到復(fù)賽譜時(shí)頻域,通過多項(xiàng)式平滑的方式提取出子波振幅譜,從而推導(dǎo)出一種基于改進(jìn)廣義S變換的復(fù)賽譜時(shí)頻域反褶積方法。

1 方法原理

1.1 改進(jìn)廣義S變換

對給定平方可積的信號x(t)的S變換[11]為

S(τ,f)=

(1)

(2)

因此得到對給定信號x(t)的一種兩參數(shù)改進(jìn)廣義S變換為

Sg(τ,f)=

(3)

式中,Sg(τ,f)為x(t)在S域中時(shí)頻譜值;λ和r為高斯窗調(diào)節(jié)參數(shù)。

改進(jìn)廣義S變換與快速傅里葉變換直接聯(lián)系,因此可以利用快速算法實(shí)現(xiàn)其逆過程:

(4)

1.2 復(fù)賽譜時(shí)頻域反褶積模型

傳統(tǒng)的時(shí)間域地震道卷積算法可以表示為

(5)

式中,地震道s(t)可以由時(shí)間域的反射系數(shù)r(t)和震源子波ω(t)的卷積表示;n(t)為干擾噪聲。

在不考慮噪聲的影響下,引入衰減函數(shù)a(τ,f)和衰減因子Q來描述地震子波在頻率域中的衰減,衰減函數(shù)表示為

|a(τ,f)|=exp(-πfτ/Q).

(6)

綜合推導(dǎo)出頻率域的衰減褶積模型[11]:

(7)

將式(7)擴(kuò)展到S域中,可以將衰減地震道的改進(jìn)廣義S變換時(shí)頻譜Sg(τ,f)近似等同于子波頻譜ω(f)、衰減函數(shù)a(τ,f)和反射系數(shù)時(shí)頻譜值Rg(τ,f)之間的乘積[10,12],即

Sg(τ,f)≈ω(f)a(τ,f)Rg(τ,f).

(8)

只考慮地震道的振幅譜,式(8)可以表示[12]為

|Sg(τ,f)|≈|ωa(τ,f)||Rg(τ,f)|.

(9)

其中衰減子波的振幅譜值|ωa(τ,f)|可以表示為|ωa(τ,f)|=|ω(f)||a(τ,f)|。

時(shí)頻域地震模型可以看作是一個(gè)時(shí)變系統(tǒng),直接求出時(shí)變子波比較困難。Rosa[6]結(jié)合Ricker研究結(jié)論認(rèn)為:假設(shè)反射系數(shù)具有白譜特性,其只影響地震道振幅譜的細(xì)節(jié),而子波決定地震道振幅譜大體形狀,因此子波可以通過平滑地震道的方式提取出來。傳統(tǒng)的子波提取函數(shù)為

(10)

為了擺脫這種反射系數(shù)是白譜特性的假設(shè),本文中將每個(gè)時(shí)間點(diǎn)的非平穩(wěn)地震道時(shí)頻譜變換到復(fù)賽譜時(shí)頻域。記衰減地震道時(shí)頻譜Sg(τ,f)變換到復(fù)賽時(shí)頻域中為Slog(τ,f),即

Slog(τ,f)=log(|Sg(τ,f)|).

(11)

由于反射系數(shù)一般都是高頻的,子波是低頻的,在復(fù)賽時(shí)頻域中能夠很好地將兩者分離,減小反射系數(shù)對子波形狀的影響,因此對復(fù)賽時(shí)頻譜使用多項(xiàng)式擬合(ployfit),平滑地震道消除反射系數(shù)的影響:

Sfit(τ,f)=ployfit(Slog(τ,f)).

(12)

換算平滑后的復(fù)賽時(shí)頻譜,近似得到子波時(shí)頻譜|ωfit(τ,f)|為

|ωfit(τ,f)|≈exp(Sfit(τ,f)).

(13)

在S域中利用Sg(τ,f)和|ωfit(τ,f)|就可以得到反射系數(shù)時(shí)頻譜值R(τ,f),再對其進(jìn)行S反變換得到時(shí)間域的r(t):

(14)

為了防止分母出現(xiàn)零值引入?yún)?shù)μ,Amax表示動(dòng)態(tài)子波時(shí)頻譜|ωfit(τ,f)|的最大值。

由于以上模型,受μ值影響很大,使得算法缺乏穩(wěn)定性,本文中對以上模型進(jìn)行優(yōu)化,得

(15)

2 模型論證

以40 Hz雷克子波與隨機(jī)反射系數(shù)褶積合成一個(gè)非平穩(wěn)地震記錄[12]。利用改進(jìn)廣義S變換對單道地震記錄進(jìn)行時(shí)頻譜分析[13],取其中一個(gè)時(shí)間點(diǎn)振幅譜,如圖1(c)黑色實(shí)線,展示出原始的非平穩(wěn)地震道振幅譜明顯受反射系數(shù)影響波動(dòng)很大。利用式(11)對其取復(fù)賽譜,結(jié)果如圖1(a)紅色實(shí)線所示;在復(fù)賽譜中反射系數(shù)對低頻的子波影響更小,用式(12)對其平滑得到結(jié)果如圖1(a)黑色點(diǎn)線所示,將平滑后的復(fù)賽譜利用(13)式恢復(fù)到時(shí)間域得到子波振幅譜如圖1(c)中黑色點(diǎn)線所示。對原始雷克子波求振幅譜,結(jié)果如圖1(c)中紅色實(shí)線,將復(fù)賽譜域提取的子波與其對比,可以明顯地得到二者吻合,驗(yàn)證了復(fù)賽譜提取子波的可靠性;依次將非平穩(wěn)地震道時(shí)頻譜轉(zhuǎn)換到復(fù)賽譜時(shí)頻域,并對每個(gè)時(shí)間點(diǎn)平滑處理得到結(jié)果如圖1(b)。對平滑結(jié)果利用式(13)換算得到時(shí)變子波時(shí)頻譜,結(jié)果如圖1(d)所示,得到時(shí)頻域子波振幅譜。

圖1 復(fù)賽譜時(shí)頻域提取子波振幅譜Fig.1 Estimated wavelet amplitude spectrum in cepstral time-frequency domain

為了檢驗(yàn)各個(gè)時(shí)頻分析方法的效果,用20～120 Hz線性遞增變化的調(diào)諧信號和120～20 Hz線性遞減變化的調(diào)諧信號組合成一個(gè)復(fù)合信號如圖2(a)[14],分別使用Gabor變換、S變換[15]、和改進(jìn)廣義S變換對其進(jìn)行時(shí)頻分析得到結(jié)果分別如圖2(b)、圖2(c)、圖2(d)?？梢钥吹紾abor變換對調(diào)諧信號的頻率變化趨勢有很好的描述,但是整體頻率分辨率不夠高,且在中間位置不能很好地將兩個(gè)不同的調(diào)諧信號分開;當(dāng)復(fù)合信號處于低頻時(shí),S變換具有很高的時(shí)頻分辨率,但隨著復(fù)合信號由低頻到高頻變化,其時(shí)頻分辨能力越來越差;取改進(jìn)廣義S變換λ=0.8,r=1.4,增加窗函數(shù)對頻率變化的敏感程度,提高時(shí)頻聚焦性,對復(fù)合信號進(jìn)行分析,結(jié)果如圖2(d)所示,對比上述兩種方法,該方法能更好地檢測到復(fù)合信號的頻率變化趨勢,時(shí)頻分辨率更優(yōu)。

圖2 復(fù)合信號時(shí)頻分析Fig.2 Time-frequency analysis of composite signals

建立理論模型驗(yàn)證改進(jìn)后的時(shí)頻域反褶積算法的穩(wěn)定性。模擬地下反射系數(shù)如圖3(a)所示,使用主頻為40 Hz的雷克子波利用式(5)得到合成單道地震記錄結(jié)果如圖3(b)所示,利用式(14)對圖3(a)單道地震記錄進(jìn)行反褶積處理,其中μ取值為10-5,得到反褶積結(jié)果如圖3(c)所示;取同一μ值利用本文改進(jìn)之后的算法式(15),得到反褶積結(jié)果如圖3(d)所示;由于提取子波時(shí)算法近似帶來的誤差和μ值不確定性的影響,導(dǎo)致反褶積結(jié)果出現(xiàn)圖3(c)中額外假震蕩,經(jīng)過改進(jìn)后的算法,提高了系統(tǒng)穩(wěn)定性,抗μ值干擾能力增強(qiáng)。

圖3 改進(jìn)算法穩(wěn)定性研究Fig.3 Stability research of improved algorithm

圖4 反褶積結(jié)果對比Fig.4 Comparison of deconvolution results

將復(fù)賽譜時(shí)頻域反褶積與Gabor反褶積結(jié)果進(jìn)行對比。生成一個(gè)隨機(jī)反射系數(shù)如圖4(a)所示,利用式(5)使用主頻為40 Hz的雷克子波與其褶積產(chǎn)生一個(gè)平穩(wěn)地震記錄如圖4(b)所示,對其時(shí)頻分析結(jié)果如圖5(a)所示;引入Q值為50的衰減函數(shù)[16],結(jié)合式(7)得到一個(gè)非平穩(wěn)地震道[17]如圖4(c)所示,其時(shí)頻分析結(jié)果如圖5(b)所示,對非平穩(wěn)地震道進(jìn)行Gabor反褶積得到的結(jié)果如圖4(d),其時(shí)頻分析結(jié)果如圖5(c)所示,圖4(e)是使用本文復(fù)賽譜時(shí)頻域反褶積處理得到的結(jié)果,其時(shí)頻分析結(jié)果如圖5(d)所示;觀察發(fā)現(xiàn)復(fù)賽譜反褶積結(jié)果整體比Gabor反褶積結(jié)果更加收斂,且與圖4(a)中的反射系數(shù)對照更加吻合;對比圖4中紅色虛線標(biāo)記的位置,可以發(fā)現(xiàn)Gabor反褶積很大程度上只是對衰減地震道能量上的恢復(fù),而基于改進(jìn)廣義S變換的復(fù)賽譜時(shí)頻域反褶積,在恢復(fù)地震波衰減能量的基礎(chǔ)上,壓縮了地震子波,能清楚地區(qū)分出兩個(gè)薄層,提高了地震分辨率[18];從圖5時(shí)頻分析結(jié)果也可以看出復(fù)賽譜反褶積結(jié)果在恢復(fù)能量的同時(shí),拓寬了頻帶,這說明本文的方法更優(yōu)?；谇懊鎸蔚赖卣鹦盘柼幚硇Ч敿?xì)分析描述基礎(chǔ)上,建立一個(gè)二維速度模型[19](圖6),對比分析兩種反褶積算法分辨率。分別在200、300 m的位置假定兩個(gè)薄層,在400 m位置設(shè)計(jì)一個(gè)楔形。將速度模型轉(zhuǎn)換成反射系數(shù),用30 Hz的雷克子波[20]與其卷積得到二維地震剖面,結(jié)果如圖6(a)所示;結(jié)合式(7)得到衰減地震剖面如圖6(b)所示,分別用Gabor反褶積和復(fù)賽譜反褶積對衰減地震剖面處理,得到結(jié)果分別如圖6(c)、圖6(d)所示?？梢灾庇^看出,兩種反褶積方法對衰減剖面的能量都有恢復(fù)作用;在200 ms位置,復(fù)賽譜反褶積處理后兩個(gè)負(fù)反射層能夠很明顯地分開,但是Gabor反褶積處并沒有完全消除子波旁瓣的影響,對薄層識別能力不強(qiáng);同樣,在300 ms位置處,復(fù)賽譜反褶積對子波壓縮的能力更強(qiáng)。注意圖中黑色豎線位置,復(fù)賽譜反褶積在第10道就可以將楔形模型的兩個(gè)反射層完全分離,而Gabor反褶積在第14道才能分離這兩個(gè)子波,這說明復(fù)賽譜反褶積方法時(shí)間分辨率更高。

圖5 不同信號時(shí)頻譜Fig.5 Time spectrum of different signal

圖6 速度地質(zhì)模型Fig.6 Velocity geological model

圖7 二維模型反褶積結(jié)果對比Fig.7 Comparison of two-dimensional seismic model deconvolution result

3 實(shí)際資料處理

為了驗(yàn)證本文方法對海上數(shù)據(jù)的實(shí)用性,選取中國南海某海域地震數(shù)據(jù)進(jìn)行復(fù)賽譜時(shí)頻域反褶積處理。研究區(qū)域構(gòu)造背景上是深水沉積水道砂巖儲層氣藏,砂巖呈條帶分布狀,且分布局限,為了更好地識別儲層,使用本文的方法進(jìn)行高分辨處理。選取目的層區(qū)塊如圖8(a)所示,圖8(b)和圖8(c)是經(jīng)過Gabor反褶積和復(fù)賽譜時(shí)頻域反褶積方法處理之后的地震剖面。總體來看,兩種方法反褶積處理后對整個(gè)剖面的同相軸都有一定的壓縮作用。在圖8中紅色橢圓標(biāo)記的位置,復(fù)賽譜時(shí)頻域反褶積處理后成功識別出兩個(gè)薄氣層,而Gabor反褶積沒有區(qū)分出這兩個(gè)薄層;在圖中兩個(gè)紅色箭頭位置,復(fù)賽譜時(shí)頻域反褶積也成功地恢復(fù)了弱反射層的能量,且與井曲線對照一致,準(zhǔn)確地識別出下方的薄氣層和油層;同時(shí)對比發(fā)現(xiàn)復(fù)賽譜時(shí)頻域反褶積處理之后的剖面整體成層性更好,與井曲線對照更加準(zhǔn)確;這說明了本文方法在南海數(shù)據(jù)高分辨處理的有效性。將這兩個(gè)剖面進(jìn)行綜合振幅譜分析得到結(jié)果如圖9所示,藍(lán)色線是原始剖面的綜合振幅譜,綠色線是Gabor反褶積處理后的剖面綜合振幅譜,紅色線是復(fù)賽譜時(shí)頻域反褶積處理后的剖面綜合振幅譜,可以看到復(fù)賽譜時(shí)頻域反褶積處理后,在保留了原始剖面的低頻信息的基礎(chǔ)上,更好地拓寬了高頻部分頻帶。

圖8 反褶積前后剖面對比Fig.8 Comparison of profile before and after deconvolution

圖9 反褶積前后剖面振幅譜對比Fig.9 Comparison of amplitude spectra before and after deconvolution

4 結(jié) 論

(1)本文中引入改進(jìn)廣義S變換,克服了傅里葉變換的時(shí)窗大小問題,彌補(bǔ)了Gabor變換時(shí)窗不能自適應(yīng)頻率變化的缺陷。

(2)在復(fù)賽譜時(shí)頻域提取子波避免了反射系數(shù)白譜特性的假設(shè),更加符合實(shí)際情況。

(3)改進(jìn)后的算法更加穩(wěn)定,能更好地恢復(fù)地層對地震波的衰減,同時(shí)提高了對薄儲層的識別能力。

(4)地層對地震波的衰減是相對的過程,將反褶積擴(kuò)展到復(fù)賽譜時(shí)頻域后,每個(gè)時(shí)刻是相互獨(dú)立的,因此不用考慮Q值的影響,直接恢復(fù)地震波的能量。