姚振岸 孫成禹 李紅星 楊安根

(①東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點學(xué)科實驗室，江西南昌 330013；②中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島，266580；③青島海洋科學(xué)與技術(shù)試點國家實驗室，山東青島，266237)

0 引言

傳統(tǒng)地震反演是基于褶積模型，也就意味著地震道是平穩(wěn)的，實際上由于吸收衰減和噪聲干擾等因素的影響，地震數(shù)據(jù)往往是非平穩(wěn)的。因而，常子波假設(shè)不適用于實際地震道。為了解決這個問題，目前已經(jīng)發(fā)展多種方法建立時變子波或非平穩(wěn)褶積模型，以滿足地震數(shù)據(jù)的非平穩(wěn)性質(zhì)。Van der Baan等[1-2]利用最大峰度提出了一種時變子波估計方法；Margrave等[3]考慮到地層的衰減提出非平穩(wěn)褶積模型，并基于此模型實現(xiàn)了高分辨率地震反褶積；Dai等[4]利用廣義S變換對地震道進行時頻分解，并基于局部相似度估計時變子波，通過反褶積驗證了時變子波的適用性，繼而又提出時頻分析與自適應(yīng)分段結(jié)合的時變子波提取方法[5]。為了獲得符合實際的混合相位子波，孔德輝等[6]提出一種基于在線字典學(xué)習(xí)的時變子波估計方法；Zhang等[7]基于局部譜時頻分析提取了時變子波，并將其用于地震反演，得到了更準(zhǔn)確的波阻抗反演結(jié)果；Wang[8]提出了廣義地震子波解析式，能更好地表征地震數(shù)據(jù)頻譜的不對稱性，用于更精確的地震子波構(gòu)建。

時頻分解技術(shù)能刻畫時間與瞬時頻率之間的非平穩(wěn)關(guān)系，已逐漸成為地震信號處理和解釋的有力工具，廣泛用于地震薄層厚度分析、氣藏低頻陰影檢測等方面。經(jīng)典的短時傅里葉變換(STFT)用窗函數(shù)截取信號，展示信號的局部頻率特征[9]，但存在固有的窗口效應(yīng)，易引起頻譜模糊。將STFT的正弦基函數(shù)改為小波函數(shù)，以克服窗口效應(yīng)，就形成了目前廣泛使用的基于小波的頻譜分解方法[10]。該類方法在描述微幅構(gòu)造和儲層特征方面都體現(xiàn)出良好的性能，常見的有連續(xù)小波變換[11]和匹配追蹤技術(shù)[12-13]。Liu等[14]和Liu等[15]發(fā)展了一種基于局部屬性的S變換方法[16]，通過正則化的非平穩(wěn)回歸[17]實現(xiàn)頻譜分解。這種基于局部屬性的譜分解技術(shù)已成功應(yīng)用于地滾波壓制、多分量數(shù)據(jù)匹配和地層解釋[15,18]。以上譜分解技術(shù)都各有其優(yōu)點，但都不可避免地在時間分辨率與頻率分辨率之間進行權(quán)衡。經(jīng)驗?zāi)Ｊ椒纸?EMD)[19]算法能將信號分離成局部常頻率分量，同時獲得較高的時間和頻率分辨率。隨后很多種改進算法被提出，如EEMD和CEEMD[20]。近年來，基于算法結(jié)合(如同步壓縮變換SST[21-22])和反演方案(如基追蹤[23])的時頻分解技術(shù)在改善時頻分辨率方面呈現(xiàn)越來越大的潛能。

本文基于基追蹤譜分解提取時變廣義地震子波，并將其應(yīng)用于非平穩(wěn)地震道的基追蹤反演，提高了反演結(jié)果的分辨率，增強了反演結(jié)果的地質(zhì)連續(xù)性。首先對非平穩(wěn)地震道做基追蹤譜分解，將其映射到時頻空間，于各時間采樣點處獲得局部功率譜；然后采用廣義地震子波構(gòu)建方法，提取時變子波并構(gòu)建時變子波核矩陣；繼而基于時變奇偶反射波形字典，最終實現(xiàn)非平穩(wěn)地震道的基追蹤反演。對比常子波基追蹤地震道反演結(jié)果，說明了時變子波基追蹤反演對反演結(jié)果分辨率和地質(zhì)連續(xù)性的改善，也證明了提取時變子波的必要性。

1 基本理論與技術(shù)流程

1.1 基追蹤譜分解

基追蹤(BP)與匹配追蹤(MP)譜分解[12]的基本原理類似，都是試圖將地震信號分解到預(yù)先定義的波形字典的各個基函數(shù)上，求得基函數(shù)表示系數(shù)，經(jīng)尺度到頻率的映射，最終得到高分辨率譜分解結(jié)果?；粉欁V分解有兩個典型特點：一是引入了一個最小化項減少重構(gòu)基函數(shù)的數(shù)量和幅值，以期得到地震信號的稀疏表示[23]；二是在匹配追蹤譜分解過程中，采用逐步壓縮的方法，也就是依據(jù)內(nèi)積最大原則逐次識別并移除基函數(shù)，而基追蹤譜分解是同時鑒別所有基函數(shù)，通過將識別和移除兩個步驟融入到一個反演問題中予以實現(xiàn)[24-25]。

地震道s(t)可被表示為一族基函數(shù)ψ(t,n)與其對應(yīng)的系數(shù)序列a(t,n)的卷積

(1)

式中：N是基函數(shù)的個數(shù)；n是控制基函數(shù)頻率特征的尺度參數(shù)。采用矩陣記法，式(1)寫作

(2)

式中：s是地震道s(t)的向量記法；Ψn表示基函數(shù)ψ(t,n)的卷積矩陣；an表示對應(yīng)Ψn的系數(shù)；D為波形字典；a是由an組成的列向量，也即由全部系數(shù)a(t,n)按列排布形成的列向量；η是隨機噪聲。這樣譜分解結(jié)果就是地震道映射在波形字典D中的權(quán)系數(shù)a(t,n)于時間—頻率空間的分布。

式(2)中權(quán)系數(shù)a的求解問題即為典型的基追蹤去噪問題。設(shè)定目標(biāo)函數(shù)為

(3)

式中λ是權(quán)衡因子。J的第一項代表基于L2范數(shù)的數(shù)據(jù)誤差項，也就是地震道與重構(gòu)數(shù)據(jù)之間的最小二乘誤差；J的第二項即為L1范數(shù)約束正則化項，λ控制數(shù)據(jù)誤差和解稀疏度的相對強弱?？捎枚喾N方法求解式(2)，使目標(biāo)泛函(式(3))最小化[26-28]。

與匹配追蹤不同，基追蹤是一種非貪婪算法，它從一個初始模型出發(fā)，通過不斷迭代，調(diào)整波形字典基函數(shù)，最終收斂到一個局部最優(yōu)解。與匹配追蹤類似，方程解的精度在一定程度上依賴于波形字典的選擇，波形字典越完備，解越精確，但同時會增加計算時間，降低計算效率[23,29]。

1.2 廣義地震子波的提取

在地震反演中，子波提取非常關(guān)鍵，子波估計的精度會直接影響地震反演的效果。當(dāng)有井信息時，往往通過井震標(biāo)定實現(xiàn)子波提??；當(dāng)沒有井信息時，往往假設(shè)反射系數(shù)序列是白色的或者近似是白色的，也就是有一個平穩(wěn)的功率譜，進而可以采用譜比法或者相關(guān)法進行子波估計。為了更好地表示實際地震道特征，Wang[8]構(gòu)建了廣義地震子波解析式，該式廣泛適用于多種非對稱振幅譜。廣義地震波解析表達(dá)式被定義為高斯函數(shù)的導(dǎo)數(shù)。廣義地震波的頻譜具有如下的解析形式

(4)

式中：τ0為對稱中心的時間位置；u為分?jǐn)?shù)階，取值范圍為0.4～2.2，它控制了頻譜的對稱性質(zhì)；ω0是參考角頻率，控制子波的帶寬和主頻。當(dāng)u=2.0時，廣義地震子波就變成經(jīng)典的雷克子波。廣義地震子波頻譜W(ω)的平均頻率及其標(biāo)準(zhǔn)差分別為

(5)

(6)

式中|A(ω)|2是實值地震信號的功率譜。分?jǐn)?shù)階數(shù)u可以通過平均頻率ωm與其標(biāo)準(zhǔn)差ωσ的比值唯一地給出

(7)

式中Γ(u)是Gamma函數(shù)。

對于這個單變量的非線性方程，式中的兩個Gamma函數(shù)的比值可被表示為漸進級數(shù)序列[30]

(8)

(9)

當(dāng)確定分?jǐn)?shù)階數(shù)u和參考角頻率ω0之后，即可求得廣義地震子波的頻譜W(ω)，然后再做反傅里葉變換，并最終得到時間域子波。

1.3 時變子波核矩陣的構(gòu)建

在地震反演中，褶積模型[31]被廣泛使用，即地震道是地震子波與反射率的褶積

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

(10)

式中：w(t)表示地震子波；r(t)表示反射系數(shù)序列；η(t)表示隨機噪聲。

地震子波與反射率的卷積運算等價于矩陣與向量相乘，即

s=Wr+η

(11)

式中：s、r和η分別代表地震記錄向量、反射系數(shù)向量和隨機噪聲向量。假設(shè)W代表子波核矩陣，該矩陣通過子波w(t)沿著方陣主對角線方向平移構(gòu)建而成，W在反演問題中即為映射矩陣。

基于褶積模型假設(shè)，采用統(tǒng)計方法[32-36]或者井震標(biāo)定的方法[37-38]，可以從地震數(shù)據(jù)中估計得到一個常子波w(t)。

考慮到吸收衰減以及噪聲干擾造成的地震道的非平穩(wěn)性，首先對地震道做基追蹤譜分解，然后逐時間點提取廣義地震子波，再用時變的廣義地震子波逐列取代子波核矩陣W的常子波w(t)，形成時變子波核矩陣，用于后續(xù)的地震反演。

2.4 時變子波基追蹤反演

基追蹤是一種L1范數(shù)約束的優(yōu)化方法，考慮到地下地層的稀疏性，近年來該方法逐漸被應(yīng)用于地震反演[39-43]，提高了反演結(jié)果的分辨率，其中目標(biāo)函數(shù)被定義為

min[‖s-Gr‖2+‖r‖1]

(12)

式中G是映射矩陣，可選為常子波核矩陣或時變子波核矩陣。

為了表征地層的層稀疏，而不是反射率稀疏，Zhang等[39]采用奇偶反射脈沖對代替地層反射率

bn,mro(t,m,n,Δt]

(13a)

或?qū)憺榫仃囆问?/p>

r=Drx

(13b)

式中：ro(t,m,n,Δt)和re(t,m,n,Δt)分別為奇、偶反射對分量；an,m和bn,m分別為對應(yīng)奇、偶反射對的系數(shù)；Δt為采樣間隔；m、n分別為層頂、底界面所對應(yīng)的采樣點位置；M為單道地震道的采樣點數(shù)；N為最大地層數(shù)，根據(jù)實際地震資料選定；Dr為反射系數(shù)字典；x為反射系數(shù)字典基函數(shù)的稀疏表示系數(shù)向量。

求反射系數(shù)序列r(t)與地震子波w(t)的卷積，得到地震道

bn,m[w(t)*ro(t,m,n,Δt)]}

(14a)

或?qū)憺榫仃囆问?/p>

s=Dsx

(14b)

即地震數(shù)據(jù)是奇、偶反射對波形的疊加。式中Ds為奇、偶反射對波形字典?？紤]到地震數(shù)據(jù)的非平穩(wěn)性，將常子波w(t)替換為時變子波w(t,m,n)，可得

bn,m[w(t,m,n)*ro(t,m,n,Δt)]}

(15a)

或?qū)憺榫仃囆问?/p>

(15b)

反演過程中，目標(biāo)函數(shù)變?yōu)?/p>

(16)

采用基追蹤優(yōu)化算法求解上述方程，即可得到優(yōu)化稀疏表示系數(shù)向量x，然后利用式(13)重構(gòu)得到反射系數(shù)的最終反演結(jié)果。

2 非平穩(wěn)基追蹤反演模型測試

為了說明非平穩(wěn)基追蹤反演的重要性，建立如圖1a所示的聲波阻抗模型，計算對應(yīng)的反射系數(shù)如圖1b所示。為了獲得非平穩(wěn)地震道，對反射系數(shù)分段采用主頻為60Hz和100Hz兩種雷克子波進行褶積(圖1c)。

分別采用60Hz、100Hz雷克子波波形字典以及由60Hz、100Hz雷克子波形成的時變子波波形字典，進行基追蹤反演，得到反射系數(shù)反演結(jié)果(圖2)。

圖1 模型及合成地震記錄

圖2 基追蹤反射系數(shù)反演結(jié)果

可以看出，采用60Hz雷克子波波形字典時，地震道中采用100Hz雷克子波褶積形成的部分不能被有效地反演出來，反之亦然。總的來說，低頻成份高頻子波反演結(jié)果使反射率變“胖”，分辨率降低，高頻成份低頻子波反演結(jié)果的反射率個數(shù)增多，出現(xiàn)虛假界面信息。采用時變子波字典時，基追蹤地震反演能獲得準(zhǔn)確的反射率序列。

因此，基追蹤地震反演中，子波字典的正確性和完備性至關(guān)重要，故對于實際的非平穩(wěn)地震數(shù)據(jù)，時變子波的有效估計也是地震反演精度的有力保障。

3 實際資料測試

3.1 基追蹤時頻分解

本文采用來自荷蘭北海F3區(qū)塊的實際地震數(shù)據(jù)體和測井?dāng)?shù)據(jù)，闡述基追蹤譜分解時變子波提取與地震反演的技術(shù)流程。圖3為過F06井(位于Crossline387)的Inline244地震剖面。對該地震剖面逐道進行基追蹤譜分解(圖4)，可以看出時頻譜分解能量團具有較高的分辨率，與地震波形吻合較好，隨著時間增加，時頻譜主頻整體變小，但局部有明顯錯動。圖5為整個地震剖面基追蹤譜分解結(jié)果的三維展示，可以看出，地震記錄的譜分解結(jié)果不僅隨時間軸變化，還沿著Crossline方向變化，在整個地震剖面每一道的每個時間點處都提取一個廣義地震子波，則對整個地震剖面而言，提取到的地震子波不僅是時變的，還是空變的。

圖3 過F06井的Inline244地震剖面

3.2 時變子波矩陣構(gòu)建

傳統(tǒng)地震子波的估計是基于一些假設(shè)，包括卷積模型、隨機反射率以及最小相位等[32-34]。如Inline244地震剖面，采用測井?dāng)?shù)據(jù)和井旁地震道做井震標(biāo)定，利用譜除法獲得常地震子波(圖6)。該地震子波在傳統(tǒng)的地震反演中用于整條地震剖面，距井越遠(yuǎn)的地震道反演結(jié)果誤差越大，且在單地震道內(nèi)，未考慮地震道的非平穩(wěn)特征，因此常子波地震反演不能精確地捕獲地層的局部特征。

從圖4和圖5可以清晰地觀測到地震數(shù)據(jù)時頻譜沿著時間軸和空間水平軸的變化特征。基追蹤譜分解結(jié)果代表地震數(shù)據(jù)沿時間和頻率的能量分布，可以看作是每個地震道在每個采樣點的局部功率譜。基于單點功率譜，采用廣義地震子波構(gòu)建方法，提取整個地震剖面的時變子波。沿圖4a中的四條白色虛線提取功率譜(圖7上)，從而得到時間點0.52、0.80、1.00、1.30s處的廣義地震子波(圖7下)。

圖4 Inline244單道地震記錄(左)及其基追蹤譜分解結(jié)果(右)

圖5 Inline244地震剖面基追蹤譜分解結(jié)果

圖6 經(jīng)井震標(biāo)定提取的常地震子波

圖7 四個時間點處的局部譜(上)及提取的時變子波(下)

可以看出，地震子波從淺到深逐漸變化，對0.12s以上的中新世—更新世地層，廣義地震子波的主瓣基本一致，旁瓣有差異；而對于0.12s以下的大型河流相基底，廣義地震子波主頻降低，子波長度有所增大。

傳統(tǒng)地震反演將常地震子波沿方陣主對角線方向設(shè)定獲得的子波核矩陣，作為聯(lián)系地震記錄與反射系數(shù)序列的映射算子。用逐點提取的時變子波替換常子波，形成時變子波核矩陣(圖8)，形成非平穩(wěn)地震道反演的基本映射關(guān)系。從圖8a可見，通過對Inline244地震剖面的每一道提取時變廣義地震子波，可以得到空變的子波核矩陣；圖8b和圖8c所示的子波核矩陣切片分別對應(yīng)圖8a中水平和垂直白色虛線。從圖8可以清晰地觀察廣義地震子波沿時間和Crossline方向的變化特征，這也驗證了基于基追蹤譜分解進行子波估計的可行性。獲得時變子波矩陣之后，對應(yīng)每個時間點上，在下一個波長范圍內(nèi)構(gòu)建楔形體奇、偶波形字典，即將傳統(tǒng)的反射系數(shù)稀疏變成層稀疏，繼而求解目標(biāo)函數(shù)(式(16))，最終反演得到反射系數(shù)。

從圖9可見，采用常地震子波合成地震記錄時，深層部分與原始地震記錄不能很好地匹配，若采用時變地震子波，整道地震記錄與原始地震記錄匹配較好。采用常地震子波時，合成地震記錄與原始地震記錄的相關(guān)系數(shù)為0.31，而采用時變子波時，其相關(guān)系數(shù)為0.52。因此，采用時變地震子波能有效地改善井震標(biāo)定效果，為后續(xù)地震反演提供保障。

圖8 時變子波核矩陣

(a)Crossline400地震記錄時變子波核矩陣; (b)沿Crossline方向水平時間0.80s處的時變子波核矩陣切片; (c)沿Crossline方向垂直時間0.80s處的時變子波核矩陣切片; (d)圖a中紅色實線矩形區(qū)域的細(xì)節(jié)展示

圖9 測井?dāng)?shù)據(jù)及井震標(biāo)定

3.3 時變子波基追蹤反演

為了比較基追蹤譜分解時變地震子波與井震標(biāo)定常地震子波的差別，用實際地震數(shù)據(jù)進行測試。反演采用相同的初始模型、收斂條件及最大迭代次數(shù)。對比圖10a與圖10b可以看出，采用時變子波可提高反演結(jié)果的分辨率，尤其對深部地層的構(gòu)造細(xì)節(jié)(圖中方框區(qū)域)展現(xiàn)得更清晰，連續(xù)性(圖中橢圓形區(qū)域)得到明顯增強。

圖10 采用不同子波Inline244地震剖面反射系數(shù)反演剖面

4 結(jié)論

(1)基于基追蹤的譜分解法將時頻分解問題融入反演框架，并結(jié)合L1范數(shù)稀疏約束，能同時獲得較高的時間和頻率分辨率，克服了傳統(tǒng)時頻分解方法中時間與頻率高分辨率之間的矛盾。

(2)廣義地震子波被定義為高斯函數(shù)的分?jǐn)?shù)階導(dǎo)數(shù)，基于地震道譜分解的時間點譜可唯一地估計廣義地震子波，它能更好地表征地震數(shù)據(jù)的局部屬性和非平穩(wěn)特征?；跁r變廣義地震子波構(gòu)建的非平穩(wěn)波形字典更精確和完備，為后續(xù)的基追蹤地震反演提供了保障。

(3)實際資料測試表明，采用本文提出的基追蹤譜分解廣義時變子波提取方法，能實現(xiàn)高效井震標(biāo)定，基于時變子波的非平穩(wěn)基追蹤地震反演結(jié)果有更高的分辨率、更強的連續(xù)性和更清晰的地質(zhì)細(xì)節(jié)展示。