江雨濛 曹思遠(yuǎn)* 陳思遠(yuǎn) 蔡明俊 張家良

(①中國石油大學(xué)(北京)油氣資源與探測(cè)國家重點(diǎn)實(shí)驗(yàn)室，北京 102249；②中國石油大學(xué)(北京)地球物理學(xué)院，北京 102249；③中國石油大港油田公司，天津 300280)

0 引言

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

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

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

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

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

1 方法原理

1.1 基于廣義S變換的時(shí)變子波提取

根據(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è)條件，可認(rèn)為反射系數(shù)的自相關(guān)函數(shù)為脈沖函數(shù)，即地震子波頻譜可從地震記錄的頻譜中獲取。因此，對(duì)式(1)等號(hào)兩邊同時(shí)取自相關(guān)，建立地震記錄自相關(guān)與子波自相關(guān)的關(guān)系式

(2)

式中“?”代表自相關(guān)算符。

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

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

(3)

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

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

信號(hào)s(t)的傳統(tǒng)S變換數(shù)學(xué)表達(dá)式為

(4)

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

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

ST*(t，f)=

(5)

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

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

圖1 基于廣義S變換提取時(shí)變子波

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

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

s=Wr+n

(6)

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

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

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

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

(7)

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

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

2 模型測(cè)試

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

2.1 時(shí)變子波估計(jì)

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

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

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

圖2 基于S變換提取時(shí)變子波

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

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

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

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

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

2.3 參數(shù)設(shè)置

圖5 λ分別取0.001λmax、0.01λmax、0.05λmax和0.1λmax時(shí)的反演結(jié)果

3 實(shí)際資料處理

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

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

圖6 實(shí)際地震資料處理效果

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

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

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

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

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

4 結(jié)論

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

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

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

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

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