應用時變子波的盲反射系數(shù)反演

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

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

0 引言

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

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

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

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

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

1 方法原理

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

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

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

(1)

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

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

(2)

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

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

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

(3)

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

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

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

(4)

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

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

ST*(t，f)=

(5)

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

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

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

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

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

s=Wr+n

(6)

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

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

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

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

(7)

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

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

2 模型測試

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

2.1 時變子波估計

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

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

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

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

2.2 盲反射系數(shù)反演

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

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

圖3 無噪聲時反射系數(shù)反演結果

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

2.3 參數(shù)設置

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

3 實際資料處理

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

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

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

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

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

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

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

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

4 結論

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

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

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

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

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