應用時變子波的盲反射系數反演

江雨濛 曹思遠* 陳思遠 蔡明俊 張家良

(①中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249；②中國石油大學(北京)地球物理學院，北京 102249；③中國石油大港油田公司，天津 300280)

0 引言

反射系數反演是提高地震資料分辨率的重要手段之一。傳統反演方法大多基于穩態褶積模型，即假設地震子波已知，且其振幅譜和相位譜在地震波傳播過程中是時不變的。由于實際地層介質的非完全彈性，地震子波在傳播過程中具有動態衰減性且是未知的。因此，根據地震數據估計反射系數不僅是非穩態的，更是一個全盲的過程[1-2]。

針對這一問題，人們提出了多種方法估計時變子波[3]，進而實現從非穩態地震數據中反演反射系數。馮晅等[4]提出了分時窗提取子波并將其用于合成地震記錄。Van der Baan[5]利用旋轉相位和峰值最大化準則成功提取了非最小相位時變子波。

上述方法大多采用分段處理，并假設每一段內地震子波是時不變的，進而從每一段內提取一個時不變子波[6-7]。這類方法的精度易受時窗長度選擇的影響，且所提取的具有平均意義的子波并不能充分反映子波的時變特征[8-10]。近年來，為了消除時窗長度對提取子波的影響，戴永壽等[11]和王蓉蓉等[12]相繼提出了基于時頻譜模擬估計時變混合相位子波的方法。在此基礎上，Zhang等[13]利用局部譜提取技術求取時變子波并應用于地震反演中。李婧等[14]和姚振岸等[15]分別利用廣義S變換和基追蹤譜分解改進了該方法，實現從非穩態地震數據中反演反射系數。這類方法的核心理論是通過對非穩態地震數據進行時頻分析處理，利用每一采樣點的頻譜提取時變子波，這是現今逐點提取子波的常用有效方法。

時頻分析技術是刻畫非穩態信號頻譜變化特征的重要工具[16]，它通過時間頻率的聯合函數準確描述信號在不同時間和頻率的能量密度和強度。短時傅里葉變換是最常用的時頻分析方法之一[17]，通過對時域信號進行加窗處理得到時頻譜圖。但由于受海森堡不確定原理的約束[18]，其分辨率精度不能達到非穩態地震信號的要求。為了解決短時傅里葉變換窗口形狀固定不變的問題，Sinha等[19]提出連續小波變換方法。S變換巧妙地將短時傅里葉變換與小波變換相結合，兼具兩者的優勢，且其窗口寬度直接與頻率相聯系[20]。為了提高S變換的適用性和準確性，Moukadem等[21]通過增加控制高斯窗函數的參數而提出了廣義S變換，使其在高頻段具有較高時間分辨率和相對低的頻率分辨率，能更好地識別信號的高頻信息，此性質符合非穩態地震數據的動態衰減特性。因此，利用廣義S變換對地震記錄進行時頻分析，能獲得具有更高時頻分辨率和更好聚焦性的時頻分析結果[22-23]，從而準確刻畫地震數據頻譜隨時間軸的變化，實現逐點提取時變子波。

考慮到地震數據的非穩態特征，本文提出一種新的時變子波提取方法。該方法利用廣義S變換的優勢，將非穩態地震數據變換到時頻域，并基于自相關理論實現子波的逐點提取，克服了分段提取時變子波方法的局限性。同時，考慮到地震數據處理過程是全盲的，即子波通常是未知的，本文利用每一時刻提取的子波重構時變子波矩陣，而不是整道地震記錄僅提取一個時不變子波矩陣，并將其應用于反演模型中，最終實現非穩態地震記錄反射系數的盲反演。此外，該方法是由數據本身驅動的，對實際地震資料具有更強自適應性，有利于更精細地刻畫地層結構，提高地震資料反映薄層的能力。

1 方法原理

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

根據傳統褶積模型[24]，地震記錄s(t)可表示為

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

(1)

式中：w(t)為地震子波；r(t)為反射系數序列；n(t)為噪聲；t為時間；“*”為褶積算符。

基于反射系數白噪的假設條件，可認為反射系數的自相關函數為脈沖函數，即地震子波頻譜可從地震記錄的頻譜中獲取。因此，對式(1)等號兩邊同時取自相關，建立地震記錄自相關與子波自相關的關系式

(2)

式中“?”代表自相關算符。

根據Wiener-Khintchine定理[13，24]，自相關函數的傅里葉變換等于信號的能量譜，則式(2)變形為

F(s?s)=F(w?w)?F(s)2=F(w)2

(3)

式中F(·)表示傅里葉變換。

根據式(3)，基于穩態假設條件的子波估計方法通過對每一道地震記錄的能量譜取算術平方根后，再利用逆傅里葉變換得到一個時不變地震子波。由于時不變子波的假設條件忽略了地震資料的非穩態特征，將其用于地震資料反演時限制了反演結果的精度。而事實上，地震數據的頻譜是隨時間變化的，為了更好地表征非穩態地震數據的時變特點，利用廣義S變換對地震記錄進行分析，準確定位每一時刻的頻率信息，得到每一點的頻譜，再根據式(3)提取每一時刻的子波，最后沿地震記錄的時間軸得到每一采樣點的子波，而不是一整條地震道僅提取一個子波。

信號s(t)的傳統S變換數學表達式為

(4)

式中：t和f分別是時間和頻率；w(τ-t，f)是高斯窗函數，控制了窗口的位置。

傳統S變換雖然實現了多尺度分辨率表達，但仍受窗口本身形狀的影響。由于基本窗函數的固定形態限制了S變換的使用范圍，Moukadem等[21]提出引入控制參數重新定義高斯窗函數變化，則改進的S變換的數學表達式為

ST*(t，f)=

(5)

式中m、p、k、r均為控制參數，由信號本身決定最佳參數組合。

廣義S變換根據地震數據自身特點通過調節參數控制時窗，適應頻率的變化，使得時頻分析結果具有較好的聚焦性，有利于準確提取子波局部譜。圖1a和圖1b分別為一道合成地震記錄及廣義S變化后的時頻分析結果，可明顯地看出頻譜隨時間變化，再根據每一時刻的子波局部譜，利用逆傅里葉變換得到其對應采樣點的時變子波。

圖1 基于廣義S變換提取時變子波

1.2 基于時變子波矩陣的盲反射系數反演

為簡化反演模型，將式(1)中的褶積模型改寫為子波矩陣與反射系數向量乘積的形式

s=Wr+n

(6)

式中：s、r和n分別代表地震記錄、反射系數和噪聲的向量；W是由地震子波w(t)沿對角線平移組成的托普利茲矩陣，因此是時不變的。

該模型是基于穩態子波的假設條件，即認為子波在傳播過程中不會隨傳播距離的增加而變化。實際上，由于地層存在吸收衰減效應，地震子波會隨傳播時間不斷變化，因此，式(6)基于時不變子波矩陣的模型已不能準確描述實際地震資料。

為了表征子波在地層傳播的非穩態物理過程，將上述方法提取的時變子波沿對角線元素排列，進而得到重構的子波矩陣W*。值得注意的是，W*為時變子波矩陣，其中每一列代表通過該列傳播時間的子波，逐列取代W中的時不變子波。雖然該矩陣不再是托普利茲形式，但它具有代表地震子波非穩態特性的物理意義，允許子波隨時間而變化。

根據地層構造特點，假設反射系數序列具有稀疏性，基于稀疏約束策略利用式(6)反演反射系數存在不確定性和多解性[25]。為了得到稀疏的反射系數，在目標函數中加入L1范數約束目標函數[26]，考慮噪聲通常是隨機的且大多是非稀疏的，因此對噪聲向量施加L2范數約束，最終反演目標函數為

(7)

式中：W*是W的轉置；λ為正則化參數。

針對式(7)中稀疏正則化目標函數的求解問題，近年來研究人員提出了許多實用的算法。本文選取固定點迭代(FPC_AS)算法進行求解[27]，其優勢在于運算效率和穩定性方面都有明顯提高，并且針對大規模數據反問題也具有良好表現。

2 模型測試

通過一組非穩態模型測試，驗證本文所述方法的有效性和優越性。分別在無噪聲和含噪聲條件下與傳統基于時不變子波的反演方法進行對比，同時討論式(7)中參數的選擇對結果的影響。

2.1 時變子波估計

首先，考慮地層吸收衰減因素的影響，建立無噪聲條件下非穩態地震記錄模型(圖1a實線)。子波初始主頻為40Hz，且隨著時間增加主頻線性減小，傳播1s后子波主頻為30Hz。從利用廣義S變換得到的時頻分析結果(圖1b)可見，地震記錄的頻譜是隨傳播時間變化的，體現了非穩態地震數據的時變特性。因此，該時頻分析結果有利于提取每一時刻地震記錄的頻譜，從而實現逐點子波的提取。

圖1c左～圖1f左分別顯示沿圖1b中四條白色虛線在150、400、600和900ms處提取的地震記錄頻譜，再利用逆傅里葉變換得到其對應的時域子波(圖1c右～圖1f右的波形)。可明顯看出，隨著傳播時間增加，提取的子波主頻逐漸減小，此現象反映了子波在傳播過程中的動態衰減特性。沿著整道地震記錄做相同處理，即可得到每一時刻的頻譜及其時域子波。圖1a中虛線為利用提取的時變子波與已知反射系數重構的地震記錄結果，可見重構地震記錄與理論地震記錄模型十分吻合，表明利用本文方法提取時變子波的準確性。為了更全面評價該方法的準確性，還計算了重構地震記錄與理論地震記錄模型的相關系數c(c=0.95)。

為了進一步說明本文方法的優越性，利用S變換對該模型進行處理，得到圖2b所示的時頻分析結果。對比圖1b與圖2b時頻譜圖可知，基于廣義S變換的時頻分析結果具有更好的能量聚焦性，這將有利于時變子波的精確提取。同樣，圖2c～圖2f分別顯示沿圖2b中四條白色虛線在150、400、600和900ms處提取的地震記錄頻譜及其對應的時域子波。對比圖1與圖2中的波形圖，可看出時頻分析的精確度會直接影響提取時變子波的效果。圖2a中虛線為利用基于S變換提取的時變子波與已知反射系數重構的地震記錄結果，可見雖然基于廣義S變換方法的重構地震記錄與理論地震記錄模型基本吻合，但在振幅能量的恢復上有一定誤差，其重構地震記錄與理論模型的相關系數為0.88。因此，通過定性和定量的比較分析，本文所提方法都具有更好的可行性和更高精確性。

圖2 基于S變換提取時變子波

2.2 盲反射系數反演

基于上述模型，分別在無噪聲和含噪聲兩種情況下，對比測試本文方法與傳統基于時不變子波假設條件的反射系數反演方法。通過不含噪聲非穩態地震記錄(圖3a)和理論反射系數模型(圖3b)，得到基于時不變子波反演反射系數結果(圖3c)及其誤差(圖3d)、本文方法反演反射系數結果(圖3e)及其誤差(圖3f)。可見基于時變子波的反演結果與理論反射系數吻合程度明顯優于基于時不變子波的反演結果，通過誤差對比分析，利用本文方法得到的結果不僅改善了反射系數位置錯位的情況，并且最大程度地恢復了反射系數的相對振幅，因此其計算結果誤差更小，所得反射系數更精確。

在上述非穩態地震記錄模型中加入信噪比為15dB的隨機噪聲，分別基于時變子波和時不變子波進行反射系數反演實驗(圖4)。利用含噪聲非穩態地震記錄(圖4a)和理論反射系數模型(圖4b)，得到基于時不變子波反演反射系數的結果(圖4c)及其誤差(圖4d)、本文方法反演反射系數的結果(圖4e)及其誤差(圖4f)。對比圖3可見，由于噪聲的加入，會給反演結果帶來一定誤差，但基于時變子波反演得到的反射系數與理論模型仍具有較好一致性。分析、對比圖4c與圖4e，可知含噪情況下基于時不變子波的反演結果較差，圖4c中存在許多小脈沖，這些假反射系數會直接降低反演結果的分辨率，該現象也說明精確時變子波的提取對反演結果的重要性，同時展示了本文方法的有效性和穩定性。

圖3 無噪聲時反射系數反演結果

圖4 含噪聲時反射系數反演結果(信噪比為15dB)

2.3 參數設置

圖5 λ分別取0.001λmax、0.01λmax、0.05λmax和0.1λmax時的反演結果

3 實際資料處理

選取圖6a所示實際地震資料，道數為50，截取時窗范圍是0.6～1.6s，采樣間隔為2ms。圖6b左側為提取的時不變子波，右側為基于廣義S變換分別在0.8、1.1和1.4s三個時刻提取的時變子波。由圖可見，因地震記錄高頻成分吸收衰減得更快，故地震子波主頻逐漸降低，波形隨傳播時間增大而變寬，該現象符合子波在地層傳播時的動態衰減特征，也反映了實際地震資料的非穩態特性。利用傳統時不變子波進行反演時，因子波欠精確，使得反演結果有較大誤差，降低了反演剖面的分辨率。

圖6c的左側為第40道(圖6a紅色虛線)地震記錄，中間和右側分別為基于圖6b中時不變子波與時變子波的反褶積結果，可見基于時不變子波的反褶積結果中存在較強噪聲，并且不能較好地保護反射系數的相對振幅值(箭頭所指)。

圖6 實際地震資料處理效果

對比基于以上兩種方法得到的完整反射系數剖面(圖7)，可見基于時變子波所得剖面(圖7a)能更清晰地反映同相軸形態；相較于時不變子波的反演結果(圖7b)，本文方法反褶積結果橫向更穩定，分辨率也得到明顯提高(表現為同相軸增加，矩形框所示)，可更好地恢復地下地層特征。

圖7 基于時不變子波(a)和時變子波(b)反演的反射系數剖面

考察分別利用時不變子波(圖8a)和時變子波(圖8b)對二維地震資料處理實例、從圖8中矩形框提取的子波(圖9)，發現用本文方法提取的子波旁瓣能量小，頻帶更寬，高頻部分能量抬升，具有更高分辨率，使處理后剖面構造更清晰(箭頭所指)，進一步提高地震資料表征薄層真實細節的能力。

圖8 基于時不變子波(a)與時變子波(b)反演的剖面實例

圖9 時不變子波(a)及其頻譜(c)與時變子波(b)及其頻譜(d)的對比

4 結論

將基于廣義S變換提取的時變子波用于子波矩陣重構，通過求解L1范數稀疏約束問題，實現從非穩態地震數據中盲反演反射系數，數值模擬和實際資料處理結果均表明：

(1)基于廣義S變換的時頻分析結果可更好地表征地震數據的局部屬性，有利于更精細地分析地層結構。

(2)由于本文方法數據驅動的自適應性與實現的便捷性，逐點提取的時變子波精度較高，符合地震數據的時變特征，且具有一定的容噪能力，為后續反射系數反演提供了保證。

(3)相較于采用時不變子波的常規反演方法，基于廣義S變換的時變子波的盲反射系數反演方法可以獲得更高精度、更高分辨率的反射系數剖面，提高地震資料表征薄層的能力。

本文方法基于一維褶積模型，理論上旨在處理疊后地震數據，目前需逐道進行處理，如何充分考慮地震數據的橫向連續性是下一步的研究方向。