周 慰， 黃 饒， 周懷來,2， 牛 聰， 巫南克

(1.成都理工大學 地球物理學院，成都 610059；2.“油氣地質(zhì)及開發(fā)工程”國家重點實驗室，成都 610059；3.中海油研究總院，北京 100028)

0 引言

傳統(tǒng)的反褶積基于地震波在地下傳播的過程中能量沒有發(fā)生變化的前提,但地震波在實際的傳播過程中主頻會隨著傳播時間的增長向低頻偏移，傳統(tǒng)的褶積模型不能有效模擬這個動態(tài)過程[1-2]。

考慮到地下介質(zhì)的衰減作用，彭才[3]假設子波是最小相位，利用非平穩(wěn)地震記錄提取時變子波；Margrave[4]將衰減地震記錄轉(zhuǎn)換到Gabor域，結(jié)合Rosa提出的平滑函數(shù)，對每個時窗振幅譜平滑，直接估計出衰減子波振幅譜，提出了一種時頻域褶積模型；王蓉蓉[5]對比各種時頻分析方法，對時變子波提取方法進行研究并總結(jié)。反褶積的關鍵是能否更好的分離子波和反射系數(shù)，在頻率域中分析關鍵在于能否更準確地提取子波振幅譜。Rosa[6]提出了直接對地震記錄的振幅譜進行譜模擬，提取子波振幅譜，實現(xiàn)了譜模擬反褶積；周懷來[7-9]在S域反褶積，結(jié)合了Rosa譜模擬提取子波是思想，對子波時頻譜進行譜模擬，提取動態(tài)子波時頻譜，實現(xiàn)S域反褶積。

地震信號的時頻譜分析方法有Gabor變換、小波變換、S變換等。其中Gabor變換受時窗的影響，不能根據(jù)實際信號的頻率變化進行靈活處理；小波變換需要對小波基進行合理選擇，高頻區(qū)分辨率差；S變換基本函數(shù)固定，高頻信號識別能力低。筆者對S變換窗函數(shù)進行擴展，提高其對高頻信號的分辨率，推導出一種改進廣義S變換，并在此基礎上使用傅里葉變換得到非平穩(wěn)地震記錄二次時頻譜，通過低通濾波提取子波時頻譜，實現(xiàn)二次時頻譜反褶積，彌補了傳統(tǒng)平滑函數(shù)擬合階數(shù)不確定性的缺陷，擺脫了反射系數(shù)滿足白譜特性的假設，提高了算法的實用性，從而達到在不影響地震資料信噪比的前提下，恢復地震波被地層吸收所引起的衰減，并提高對薄儲層地識別能力。

1 方法原理

1.1 改進廣義S變換

由于S變換基本窗函數(shù)基本形態(tài)固定，不能根據(jù)實際情況調(diào)節(jié)窗口的大小，近年來廣大專家學者提出很多改進廣義S變換，以達到提高時頻聚焦性的目的。為了彌補S變換在信號高頻部分分辨率低的缺陷，這里將其高斯窗函數(shù)進行擴展，令σ=λ/|f|r，使得當信號處于高頻時，加快窗口的收縮速度[10]。

(1)

式中：t為時間;f為頻率;λ、r用來調(diào)節(jié)窗函數(shù)的寬度和高度。

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

exp(-i2πt)dt

(2)

1.2 二次譜反褶積模型

由于地震波在地下傳播的過程中會受到地下介質(zhì)的影響，主頻向低頻偏移。為了更加符合這一實際情況，筆者引入一個時間延遲因子τ和衰減函數(shù)來描述這個衰減的動態(tài)過程，推導出一個頻率域的衰減褶積模型[11]：

(3)

其中：ω(f)是震源子波的振幅譜；r(t)是時間域的反射系數(shù)；α(τ，f)用來描述地震子波在頻率域中的衰減，表示為式(4)。

a(τ，f)=exp(πfτ/Q)

(4)

將式(4)擴展到S域中，可以將衰減地震道的時頻譜G(τ，f)近似等同子波頻譜ω(f)、衰減函數(shù)a(τ，f)和反射系數(shù)改進廣義S變換后時頻譜R(τ，f)之間的乘積，如果只從振幅譜出發(fā)則可以表示成式(5)。

|G(τ，f)|≈|ω(f)‖a(τ，f)‖R(τ，f)|

(5)

非平穩(wěn)地震記錄的動態(tài)子波振幅譜W(τ，f)可以用其靜態(tài)子波振幅譜ω(f)與衰減函數(shù)振幅譜a(τ，f)的乘積表示，即：

|G(τ，f)|≈|W(τ，f)‖R(τ，f)|

(6)

基于Rosa提出的傳統(tǒng)的子波提取函數(shù)為式(7)[12]。

(7)

由于式(7)譜模擬時受數(shù)據(jù)體干擾大，平滑函數(shù)受擬合階數(shù)k、n的不確定性影響，當多項式擬合階數(shù)低時，擬合誤差大，階數(shù)高時計算效率低且容易溢出。為了彌補這一缺陷，我們將非平穩(wěn)地震記錄時頻譜轉(zhuǎn)換到二次時頻譜譜中，通過濾波的方式提取子波時頻譜[13]。提取非平穩(wěn)地震記錄時頻譜|G(τ，f)|每個時間點振幅譜，對其傅里葉變換得到二次時頻譜，其結(jié)果記為|G(τ，f)|(2)，則式(6)可以表示成為式(8)。

(|R(τ，f)|)(2)

(8)

其中:|W(τ，f)|(2)表示動態(tài)子波的二次時頻譜；|R(τ，f)|(2)表示反射系數(shù)的二次時頻譜。

基于Rosa理論，地震子波振幅譜是光滑的，子波和反射系數(shù)是可以分開的，在非平穩(wěn)地震記錄二次時頻譜中，低頻穩(wěn)定的部分是子波的能量，而高頻擾動的部分是由于反射系數(shù)引起的。因此可以通過濾波的方式來保留低頻的子波能量，再通過相位恢復，反變換得到地震子波時頻譜。這里記低通濾波后的地震記錄二次時頻譜為|G(τ，f)|fit(2)，通過對其傅里葉反變換即可以提取動態(tài)子波振幅譜：

|W(τ，f)|fit≈ifft(|G(τ，f)|(2))

(9)

則反射系數(shù)時頻譜值R(τ，f)可以表示成式(10)。

(10)

為了防止分母出現(xiàn)零值引入?yún)?shù)μ，Amax表示動態(tài)子波時頻譜|ωa(τ，f)|fit的最大值。再對R(τ，f)進行S反變換就可以獲得時間域的反射系數(shù)r(t)為式(11)。

(11)

2 模型論證

使用30 Hz到120 Hz線性遞增變化的chirp信號與60 Hz到250 Hz二次函數(shù)遞增的chirp信號組合成一個復合信號如圖1(a)所示。利用圖1(a)來檢測Gabor變換、S變換和改進廣義S變換(GST)的時頻分析效果。我們可以通過調(diào)節(jié)GST參數(shù)，調(diào)節(jié)時頻分析結(jié)果,圖1(b)、圖1(c)、圖1(d)相應的是Gabor變換、S變換(λ=1、r=1)、改進廣義S變換(λ=2，r=0.8)對圖1(a)的時頻分析結(jié)果。與圖1(b)中的Gabor變換相比，圖1c中的S變換分析結(jié)果在低頻時分辨率高，隨著信號頻率的升高，時頻分辨率降低，圖1(d)中的GST方法能精確的檢測到兩個chirp信號的時間-頻率變化趨勢，從低頻到高頻，時頻分析結(jié)果優(yōu)于Gabor變換和S變換。

圖1 復合信號時頻分析對比圖Fig.1 The comparison of time-frequency analysis of composite signal(a)Chirp信號;(b)對a短時傅里葉變換結(jié)果;(c)對aS變換結(jié)果;(d)對a改進廣義S變換結(jié)果

圖2 模擬地震記錄二次時頻譜反褶積結(jié)果圖Fig.2 The quadratic time-spectrum deconvolution results of simulation seismic records(a)模擬反射系數(shù)；(b)平穩(wěn)地震記錄；(c)非平穩(wěn)地震記錄；(d)反褶積之后的地震記錄

為了檢驗二次時頻譜反褶積實用性，建立一維模型對其分析。如圖2(a)是一個加微弱噪聲的隨機反射系數(shù)，將其與主頻為40 Hz的雷克子波進行卷積得到平穩(wěn)地震記錄(圖2(b))；為了滿足實際衰減地震模型，引入Q取值50的衰減函數(shù)，結(jié)合式(3)得到衰減地震道如圖2(c)所示。將衰減地震記錄進行改進廣義S變換分析得到其時頻譜如圖4(b)，取其中一個時刻的時頻譜值進行單獨子波提取(圖3)，圖3(a)中藍色線表示原始的非平穩(wěn)地震道振幅譜，其明顯受反射系數(shù)影響很大，將其傅里葉變換得到二次譜如圖3(a)紅色線表示，將圖3(a)中的二次譜進行低通濾波并反傅里葉變換得到提取后的子波振幅譜(圖3(b)中紅色線表示)，與原始子波振幅譜(圖3(b)綠色線表示)進行對比能很好地得到吻合，這說明二次譜提取子波振幅譜的準確性相比能量有明顯的衰減；再用二次譜模擬(式(8)～式(9))依次對每個時間點進行二次時頻譜提取，得到衰減子波時頻譜如圖4(c)所示，再利用(式(10)～式(11))并結(jié)合地震道時頻譜實現(xiàn)二次時頻譜反褶積得到結(jié)果如圖2(d)所示。可以發(fā)現(xiàn)，經(jīng)過本文的反褶積方法處理后，原始衰減剖面的能量不僅得到了恢復，同時與原始的平穩(wěn)地震記錄相比結(jié)果更收斂；對比圖2中兩個紅色框位置可以看出，原始平穩(wěn)地震記錄的兩個疊加在一起的地震波被有效的分開,且能和圖2(a)相對應，驗證了本文的方法能在恢復地震波的能量的基礎上，提高對薄層地識別能力。

圖3 子波振幅譜提取結(jié)果圖Fig.3 The result of wavelet amplitude spectrum extraction(a)子波二次譜；(b)提取子波振幅譜對比

圖4 不同信號時頻譜Fig.4 The time spectrum of different signal(a)平穩(wěn)信號時頻譜；(b)非平穩(wěn)信號時頻譜；(c)對(b)平滑后信號時頻譜；(d)反褶積后時頻譜

圖5 二維地震模型反褶積結(jié)果對比圖Fig.5 The comparison of two-dimensional seismic model deconvolution results(a)模擬反射系數(shù);(b)平穩(wěn)地震記錄;(c)衰減地震記錄;(d)反褶積地震記錄

設計一個二維楔形模型驗證算法的魯棒性，如圖5(a)設定一個二維反射系數(shù)，在200 ms和300 ms的位置分別是是一負一正的兩個薄層，400 ms位置建立一個楔形模型；用一個主頻為35 Hz雷克子波和圖5(a)中的反射系數(shù)進行褶積得到地震記錄如圖5(b)；利用式(3)得到一個衰減二維地震剖面(如圖5(c))，其中衰減函數(shù)Q值取50；圖5(d)是用本文方法對圖5(c)中的衰減地震記錄進行反褶積得到的結(jié)果。對比分析可以得出結(jié)論：經(jīng)過反褶積處理后整個剖面被衰減的能量得到了恢復，子波被有效壓縮，提高了地震記錄的縱向分辨率。在楔形模型位置，觀察到原始地震剖面(圖5(b))在第19道之后兩個子波才能被分離(紅色實線位置)，經(jīng)過反褶積處理后的剖面(圖5(d))在11道(紅色實線位置)之后就可以完全的將兩個子波識別出來，這驗證了本文的方法可以有效地提高地震資料的時間分辨率。綜上模型論證說明，基于改進廣義S變換二次時頻譜反褶積能夠恢復地震波能量，同時壓縮子波，提高地震記錄分辨率。

3 實際資料處理

圖6 反褶積前后剖面對比圖Fig.6 The comparison of the deconvolution result(a)原始地震剖面；(b)Gabor反褶積結(jié)果；(c)二次時頻譜反褶積結(jié)果

圖7 井資料解釋對比圖Fig.7 The comparison of well data to profile(a)原始地震剖面；(b)Gabor反褶積結(jié)果；(c)二次時頻譜反褶積結(jié)果

圖8 反褶積前后剖面振幅譜對比Fig.8 The comparison of the amplitude spectra before and after deconvolution

選取中國南海流花工區(qū)數(shù)據(jù)進行二次時頻譜反褶積處理。流花構(gòu)造位于珠江口盆地白云凹陷東部，北部為白云東凹，處于白云東區(qū)深水扇砂體成藏帶與東沙隆起復式成藏帶交界處，為了達到更好地識別儲層的目的，使用本文的方法進行高分辨處理。選取目的層區(qū)域如圖6(a)所示，對其進行使用Gabor反褶積和二次時頻譜反褶積處理得到結(jié)果圖6(b)、圖6(c)，可以發(fā)現(xiàn)，二次時頻譜反褶積處理之后整個剖面的同相軸被壓縮，能量聚焦性更好，原始不明顯的地層經(jīng)過處理后更加地清楚，出現(xiàn)了更多的層位。為了證實處理結(jié)果的正確性，將圖6中紅色長方形標記的過井目的層，結(jié)合密度測井曲線進行井震標定，并放大如圖7(a)、圖7(b)和圖7(c)所示。對比原始剖面和Gabor反褶積處理的剖面(圖7中紅色橢圓位置)，二次時頻譜反褶積處理后層位與井曲線對照更加準確且清晰，在2.5 s位置，原始剖面兩氣層之間的薄蓋層也被準確分辨出來。結(jié)合測井資料我們對地震剖面進行解釋，從二次時頻譜反褶積處理后的剖面可以準確的區(qū)分出四個氣層和兩個油層，其中Gas2和Gas4是薄氣層，而原始剖面只能區(qū)分出兩個氣層和兩個油層，這為后期儲層預測提供借鑒。

通過對剖面進行綜合振幅譜分析，定量評估了兩種反卷積方法的應用效果，對圖7的三個剖面求平均振幅譜如圖8所示。其中二次時頻譜反褶積后的頻譜(紅色線)比Gabor反褶積數(shù)據(jù)(綠曲線)的頻譜更寬，且都優(yōu)于原始數(shù)據(jù)的頻譜(藍色線)。另外，對比三個頻譜的主頻，二次時頻譜反褶積的主頻明顯高于Gabor反褶積高于原始剖面，說明本文的方法應用效果更優(yōu)。

4 結(jié)論與認識

1)引入改進廣義S變換，將時間域反褶積擴展到時頻域，更好適應地震波傳播時的時變性。

2)基于改進廣義S變換的二次時頻域反褶積后，能有效地恢復地層對地震波的衰減。同時能有效地壓縮子波，提高對薄儲層的識別能力。

3)利用非平穩(wěn)地震記錄二次譜提取子波擺脫了反射系數(shù)要滿足白譜特性的假設，同時沒有平滑函數(shù)階數(shù)不確定性的問題。

4)由于不同單道記錄的時頻譜最大值不同，可以求出整個剖面時頻譜最大值的平均值，以此來代替Amax進行反褶積效果更好。