基于時變子波的基追蹤地震反演

周玉毅, 周懷來,2, 廖璐瑤, 王元君,2, 張 舜

(1.成都理工大學 地球物理學院，成都 610059；2.“油氣藏地質及開發工程”國家重點實驗室，成都 610059)

0 引言

地震子波提取的研究由來已久，許多研究者相繼提出和發展了多種提取地震子波方法[1-3]。地震子波估計準確程度直接影響著反演的精度與可靠性，傳統的地震子波估計方法應用前提是,地震波在傳播的過程中其波形以及能量保持不變，然而由于地下介質對地震波能量的衰減吸收，導致地震子波在地層傳播過程中將產生高頻成分缺失、相位畸變以及頻散現象，即實際地震記錄中的地震子波具有時變性，因此基于Robinson平穩褶積模型估計子波的理論是不完善的。

針對地震子波的非平穩特性，Clarke[4]首次提出非平穩褶積模型;Margrave等[5]在非平穩褶積模型基礎上研究了地震子波具有時變性的原因;彭才等[6]建立符合地震子波衰減特性的動態褶積模型，在隨機白噪反射系數和地震子波最小相位的假設前提下，通過對非平穩地震道進行加高斯窗傅里葉變換后，可以獲得地震道振幅譜，再經方形窗平滑處理得到動態子波；Margrave等[7]在Gabor變換的基礎上，模擬各個時窗內地震記錄振幅譜提取動態子波。上述方法應用的前提是反射系數序列為白噪聲的假設，后續研究者采用譜模擬技術對其進行規避。Rosa等[8]結合Ricker的研究結果，認為子波在頻率域內的曲線趨勢是光滑連續的，與之相對應的是反射系數在頻率域內是振蕩的，地震記錄的振幅譜主要受地震子波影響，而反射系數只是對細節部分有所影響，因此可以通過模擬地震記錄振幅譜來提取子波。由于譜模擬方法對于反射系數序列有很好的適應性，因此，諸多學者將該方法拓展到高精度時頻分析方法中來提高時頻聚焦性。李振春等[9]利用S變換將地震記錄轉換到頻率域，在頻率域內用數學方法的多項式擬合出時變子波的振幅譜，在地震子波相位為零的假設下，得到頻率域中的時變子波;張漫漫等[10]利用廣義S變換結合時頻濾波求取時頻譜，對時頻譜進行多項式譜模擬來提取時變子波;Zhou等[11]對譜模擬反褶積方法進行推廣，提出基于廣義S變換的動態反褶積算法，估算出的子波準確性有所提高。在考慮實際地層吸收衰減基礎上，時變子波譜可近似表示為震源子波譜與衰減函數的乘積形式，因此可以從實際地震記錄中提取時變子波進一步提高反射系數反演精度。姚振岸等[12]對非平穩地震道進行基追蹤譜分解獲得各時刻的局部功率譜，采用廣義地震子波構建方法提取時變子波，最終實現非平穩地震道的基追蹤反演;李婧等[13]基于改進廣義S變換提取符合地震數據非平穩特性的時變子波，建立時變子波矩陣，通過求解對應的稀疏約束正則化目標函數獲得與實際接近的反射系數。

筆者以前人的研究為基礎，利用廣義S變換譜模擬方法從地震資料中提取時變子波，用估算的時變子波逐點替換時不變子波構建時變子波庫，對地震記錄采用時變子波進行基追蹤反演。理論模型和實際工區地震資料測試結果表明，基于本文方法提取的時變子波基追蹤地震反演的結果,與傳統時不變子波基追蹤反演的結果相比，井震標定吻合度更高、橫向連續性更好以及細節展示更為豐富，說明對非平穩地震信號進行基追蹤反演時采用時變子波庫的有效性。

1 方法原理

1.1 基于廣義S變換的時變子波提取

Rosa結合Ricker的研究結果分析認為：子波振幅譜的曲線趨勢是光滑連續的，與之相對應，反射系數振幅譜的曲線趨勢是波動的。因此，只要設置適當的擬合參數，就能采用數學手段對地震記錄的振幅譜利用多項式擬合進行處理，最后得到子波振幅譜[8]。現在假設一個地震子波的振幅譜W(t,f)曲線趨勢與雷克子波單峰光滑連續的振幅譜相似，則可以將其擴展為基于時頻變換的時變子波，從而在頻率范圍內構成一種可用的子波數學分析模式：

(1)

其中:f為頻率；k、n為常數；an為關于f的多項式的系數。通常情況下，0≤k≤3，4≤n≤7[9]。

描述信號時間和頻率聯合分布的函數有很多種，如STFT、CWT、Hilbert-Huang Transform等，為了實現局部時頻分析，Stockwell等[14]繼承和擴展了STFT和CWT的優勢方面，提出了一種與其他時頻分析方法相同的無損可逆S變換，該方法對信號進行處理后，可以觀察到信號在頻域和時域的能量分布。設x(t)為平方可積的函數，則x(t)的S變換可以定義為：

(2)

式中：τ表示時間;f表示頻率;ST(x(τ,f))表示x(t)的S變換。

其中S變換的基本小波函數定義為：

gf(t)exp(-i2πft)

(3)

式中gf(t)為高斯窗函數，可定義為式(4)。

(4)

由于基本小波窗函數形態固定，S變換不能夠依據實際情況來調整窗口的大小，在高頻部分分辨率低，導致在實際應用中不能得到較好的效果。針對S變換的問題，陳學華等[15]在Stockwell的S變換基礎上提出了廣義S變換(Generalized S Transform)，通過加入兩個參數rgs和ρ來調節高斯窗函數的大小，其數學表達式為式(5)。

exp(-i2πft)dt

(5)

通過調整控制rgs和ρ的大小，該方法能夠對所有的時頻譜具有較高的時頻聚焦性，可以用來提取不同時間點的時變子波，利用廣義S變換時頻譜提取子波過程(圖1)如下：

圖1 時變子波提取流程圖

1)對沿t方向的非平穩地震記錄進行廣義S變換，得到時頻譜|GST(x(τ,f))|。

2)固定某一個時間T，通過設置合理的參數k和n，采用數學手段對該時間點地震記錄的振幅譜|GST(x(T,f))|利用多項式擬合進行處理，擬合得到相應的子波振幅譜|GST(W(T,f))|。

3)改變τ的取值，逐點計算子波振幅譜，從而實現“時變子波”的時頻譜提取。

4)在子波相位為零的假設前提下[9]，對各時刻的子波振幅譜進行傅里葉逆變換后，可求得各時刻零相位時變子波時域表達式。

1.2 時變子波基追蹤反演

地震記錄通常是由一個地震子波和地層內部反射界面的反射系數進行褶積而得到的，用公式表示如下：

s(t)=w(t)*r(t)

(6)

其中：w(t)為一個地震子波；r(t)為一個反射系數序列；“*”為褶積運算的符號。

將式(6)寫為矩陣形式:

s=Wr

(7)

其中：假設地震數據的采樣點為M；s是由地震數據所構成的M×1矩陣;W是由地震子波w(t)經過逐點移動構成的M×M對角矩陣；r是由地下反射系數構成的M×1矩陣。雖然式(6)和式(7)一個是褶積運算，一個是矩陣運算，但是它們的運算結果是完全相等的。

基追蹤是一種基于范數約束的最優化準則，基追蹤反演結果與實際地層相接近，分辨率較高，近年來逐漸被應用于地震反演，其中目標函數被定義為式(8)[16]。

J=argmin||s-Wr||2+λ||r||1

(8)

其中：J為目標函數；下標1和2為目標向量，分別用于表示L1和L2范數；λ是影響調整值的因子。滿足J→min即為一個解函數向量r的最終函數值。

為了表征地下介質的稀疏性，而不是反射系數的稀疏，根據Helholm定理，將原始反射系數對分解成如圖2所示的偶脈沖對與奇脈沖對的線性組合，則任意反射系數序列可表示為式(9)。

bm,nro(t,m,n,Δt)}

(9)

其中：ro(t,m,n, Δt)為奇楔反射系數對分量；re(t,m,n, Δt)為偶楔反射系數對分量；am,n和bm,n為對應奇、偶反射系數對對應的系數；Δt為采樣間隔；m、n分別為地下介質層頂、底界面所對應的采樣點位置；由于不同反射對出現的位置不固定，頂底反射系數以平移量mΔt沿著圖2中t的方向進行逐點移動，其中m的取值范圍為“0”到M，M、N由實際地震資料給定。

根據褶積定理，將反射系數對與地震子波進行褶積可得到地震信號:

re(t,m,n,Δt)}+bm,n[w(t,m,n)*

ro(t,m,n,Δt)]}

(10)

用矩陣的形式表示為式(11)。

s=Dwr

(11)

也就是說，地震響應是通過奇、偶楔形反射系數模型和地震子波褶積獲得的奇、偶地震響應的線性組合。式(11)中，Dw是奇、偶楔形模型地震響應。傳統的基追蹤反演是在褶積模型假設的前提下，用統計方法或井震標定方法從實際地震數據提取出地震子波w(t)。由于隨機噪聲干擾、地下介質吸收衰減作用的影響，從實際地層采集的地震數據呈現非平穩特性，運用本文方法逐道提取各時刻的子波，逐列替代時不變子波構建時變子波核矩陣，可得：

re(t,m,n,Δt)]+bm,n[w(t,m,n)*

ro(t,m,n,Δt)]}

(12)

反演過程中，目標函數變為：

(13)

采用基追蹤最優化準則從式(12)中求得參數am,n和bm,n，因為式(9)和式(12)采用相同的系數，因此將am,n、bm,n代入式(9)通過求和重構可得到有效的反射系數。在準確提取地震時變子波的前提下，通過基追蹤最優化準則可求得最優的反射系數值。

2 模型測試

2.1 不含噪聲模型

為了檢驗時變子波庫的建立對于基追蹤反演的重要性，設計如圖3(a)所示以1 ms為采樣間隔且時間長度為1 000 ms的反射系數序列，其中各反射系數脈沖的值分別設定為(120 ms，-0.15)、(165 ms，0.09)、(380 ms，-0.06)、(540 ms，-0.06)、(760 ms，0.13)、(900 ms，0.18)，非平穩地震記錄由主頻從60 Hz單調下降到40 Hz的Ricker子波與設計反射系數褶積合成，如圖3(b)所示。分別采用60 Hz子波字典、50 Hz子波字典、40 Hz子波字典以及由本文方法提取的隨時間變化子波構建的子波字典，對合成的非平穩地震記錄進行基追蹤反演，其反演結果如圖4所示。由圖4(a)～圖4(c)可知，反演所用子波的主頻越接近合成地震記錄所用的子波主頻時，反演出的反射系數才更接近實際，否則反演出的反射系數的值與實際相差較遠，甚至出現虛假的反射系數。由此可知，針對非平穩地震記錄而言，若采用固定主導頻率的時不變子波字典用于基追蹤反演，則難以確保反演結果的準確度和可靠性，而采用時變子波字典進行基追蹤反演時能獲得準確的反射系數系列(圖4(d))，說明時變子波的估計準確程度對于基追蹤反演精度有重要影響。

圖3 反射系數序列及合成地震記錄

圖4 基追蹤反演結果

2.2 含噪聲模型

為檢驗本文方法的抗噪性，對非平穩地震道分別添加SNR=15 dB和SNR=25 dB的隨機噪聲后,得到如圖5(a)與圖5(d)所示的含噪非平穩地震記錄，對加噪后的非平穩地震道采用上述方法處理,可得到含噪非平穩地震道的時頻譜和“時變子波”時頻譜。對圖5(b)和圖5(e)進行分析可知，廣義S變換在含噪的情況下依舊具有較好的時頻聚焦性。利用本文方法分別對含噪地震數據提取時變子波，進而基追蹤反演，反演結果如圖6所示，將反演出的反射系數與原始反射系數作對比可知，時變子波基追蹤反演結果與原始反射系數結果近似，由此可以得出本文子波提取及反演方法具有較好的抗噪性。

圖5 含噪地震記錄及時頻譜

圖6 含噪地震記錄基追蹤反演結果

3 實際應用

這里所用數據為南海某實際工區的疊后過井地震數據剖面如圖7所示，該地震剖面共有201道，每道采樣點數為256，截取的時窗范圍為2 401 ms～2 912 ms，其中測井曲線為密度曲線。

圖7 實際地震數據剖面

選擇圖7中紫色虛線標注的第86道地震記錄(圖8(a))進行廣義S變換處理后，可以獲得如圖8(b)所示的時頻譜圖，對時頻譜進行譜模擬后得到如圖8(c)所示的 “時變子波”時頻譜圖。對比圖8(b)和圖8(c)可知，時變子波時頻譜消除了地層反射系數的影響，與地震波形吻合度較高。

圖8 第86道地震記錄及時頻譜圖

利用廣義S變換譜模擬方法對該地震剖面逐道提取時變子波，其不僅沿時間軸方向變化，而且還沿著聯絡測線的方向變化。在整個地震剖面中，對每道地震記錄的所有時間點提取隨時間變化的地震子波，而對于地震剖面而言，提取的地震子波隨時間變化的同時，還隨空間變化。考慮到地震子波在不同道之間的差異，以提取的時變子波為基礎，對所估計的地震子波采用多地震道加權空間變化處理[17]。為了研究時不變子波和時變子波對反演結果精度的影響和差異，在此基礎上采用本文方法提取的時變子波和時不變子波,用于實際工區地震資料的基追蹤反演，反演結果如圖9所示。對比圖9(a)與圖9(b)可以看出，與時變子波基追蹤結果圖中的綠色方框部分相比，時變子波基追蹤地震反演結果的分辨率要高于時不變子波基追蹤反演結果的分辨率；時不變子波基追蹤反演出的反射系數出現虛假反射，時變子波基追蹤反演出的反射系數與原始剖面更為接近；同時，與時不變子波基追蹤結果圖中的綠色方框部分相比，時變子波基追蹤反演出的反射系數橫向連續性更好，究其原因，是因為本文方法在求取時變子波時還考慮了地震數據縱橫向的變化特性。

圖9 不同子波字典基追蹤反演剖面對比圖

4 結論

1)本文方法提取的時變子波符合地震子波在地下傳播時的頻率衰減特性，利用時變子波建立的子波庫具有較高的精度，為后續的基追蹤反演提供了保證。

2)筆者提取時變子波過程中，沒有考慮子波的相位影響，后續可以結合子波相位的估計來提取時變子波。

3)實際地震剖面的反演結果表明，基于本文方法提取的時變子波基追蹤地震反演的結果，與傳統時不變子波基追蹤反演得到的反演結果相比，能進一步提高反演精度，可用于后期儲層預測與精細刻畫。