江雨濛,曹思遠(yuǎn),陳思遠(yuǎn),馬敏瑤,鄭 鐸,黃 芳,曹國明

(1.中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室,北京102249;2.中國石油大學(xué)(北京)地球物理學(xué)院,北京102249;3.中國電子科技集團(tuán)第二十九研究所,四川成都610036;4.中國石油大港油田公司,天津300280)

地震波在地下介質(zhì)中傳播時,由于地層吸收衰減和噪聲干擾等因素的影響,子波高頻成分能量衰減更快,因此實際采集的地震數(shù)據(jù)往往是非穩(wěn)態(tài)的,且地震子波的頻譜形態(tài)會隨傳播時間變化。由于地震子波是反射地震勘探的重要參數(shù)[1-2],也是地震資料處理和解釋的基礎(chǔ),在很多環(huán)節(jié)中都起著關(guān)鍵作用,因此地震子波提取方法的精度直接影響后續(xù)地震資料處理及解釋的精度和準(zhǔn)確性。傳統(tǒng)的地震子波提取方法大多基于時不變子波的假設(shè)條件,忽略了地震數(shù)據(jù)的非穩(wěn)態(tài)特征,使得子波估計結(jié)果存在誤差,進(jìn)而限制了后續(xù)反演結(jié)果的分辨率。為了滿足高精度地震勘探的需求,如何從非穩(wěn)態(tài)地震數(shù)據(jù)中準(zhǔn)確提取時變子波成為近年來研究的熱點之一[3-4]。

針對上述問題,前人提出了基于分段處理的時變子波估計方法[5-10],即首先將地震記錄劃分為若干段,并假設(shè)每一段內(nèi)地震子波是穩(wěn)態(tài)的,進(jìn)而結(jié)合不同方法提取每一段的時不變子波。此類方法依賴于分段方法和分段長度的選擇,且每段提取的平均意義上的子波與實際子波存在一定的誤差,使得后續(xù)反演結(jié)果不精確[11]。近年來,為了消除時窗長度對提取子波的影響,提高子波估計精度,王蓉蓉等[12-13]和張漫漫等[14]相繼提出基于時頻譜模擬估計時變混合相位子波,實現(xiàn)逐點提取子波。在此基礎(chǔ)上,李婧等[15]利用廣義S變換改進(jìn)了該方法并將提取的子波應(yīng)用于地震反演。此類方法利用經(jīng)驗公式擬合時變子波振幅,要求子波振幅譜滿足類似雷克子波的假設(shè)條件,限制了子波振幅譜的形態(tài),因此實際應(yīng)用存在局限。姚振岸等[16]通過點譜模擬提取廣義地震子波繼而實現(xiàn)非穩(wěn)態(tài)基追蹤反演。ZHANG等[17]利用局部譜提取技術(shù)求取時變子波并應(yīng)用于地震反演,該方法基于自相關(guān)技術(shù),要求反射系數(shù)是白噪的,但大多數(shù)反射系數(shù)并不符合這一假設(shè)條件。

時頻分析技術(shù)能刻畫時間與瞬時頻率之間的非平穩(wěn)關(guān)系,是分析非穩(wěn)態(tài)數(shù)據(jù)的重要工具之一[18-20]。S變換[21-23]將短時傅里葉變換和小波變換相結(jié)合,既滿足了多分辨率的特點,又避免了小波變換的容許條件,其實現(xiàn)過程具有較好的便捷性與適應(yīng)性。同步擠壓變換[24-26]是一種較新的具有較高分辨率的時頻分析方法,在原有時頻分析基礎(chǔ)上按照某種關(guān)系重排壓縮,獲得更高精度更好聚焦性的時頻譜圖。TAO等[27]提出了一種基于二階自適應(yīng)同步擠壓S變換(ASST2)的時頻分析方法,通過改進(jìn)瞬時頻率的計算公式實現(xiàn)高精度時頻定位,其時頻分析結(jié)果能量聚焦性更強(qiáng),有利于精細(xì)刻畫地震記錄每一時刻的頻率分量。該方法已經(jīng)成功地應(yīng)用于地震處理中的譜分解和面波壓制環(huán)節(jié)。

本文提出了一種基于二階自適應(yīng)同步擠壓S變換(ASST2)的時變子波提取方法,該方法從非穩(wěn)態(tài)地震數(shù)據(jù)中直接逐點提取子波,避免了誤差累計,同時,對于實際地震資料具有更好的適應(yīng)性。首先利用ASST2技術(shù)對地震記錄進(jìn)行譜分解,獲取每一時刻的局部譜;然后通過建立頻譜統(tǒng)計參數(shù)與加權(quán)指數(shù)型子波(FWE)[28]之間的數(shù)學(xué)關(guān)系,根據(jù)每一時刻的頻譜信息得到對應(yīng)的時變子波及其解析式,而不是每一條地震道提取一個時不變子波;最終得到基于時變子波矩陣的非穩(wěn)態(tài)地震記錄模型,為后續(xù)反演提供理論基礎(chǔ)。本文方法規(guī)避了傳統(tǒng)子波估計方法要求地震記錄穩(wěn)態(tài)或分段穩(wěn)態(tài)的假設(shè)條件,從時頻域逐點提取子波,充分考慮了地震資料的非穩(wěn)態(tài)特性。此外,該方法是數(shù)據(jù)驅(qū)動的,克服了常規(guī)方法對地震子波形態(tài)的限制條件,同時精確的數(shù)學(xué)推導(dǎo)為本文方法提供了理論依據(jù),提高了估計子波的精度,為后續(xù)實現(xiàn)高精度地震資料處理及解釋提供了保障。

1 方法理論

1.1 ASST2時頻分析方法

時頻分析方法是研究非穩(wěn)態(tài)信號的有效工具,可以準(zhǔn)確定位每一時刻的頻率成分,獲取地震記錄每一采樣點的子波振幅譜,進(jìn)而實現(xiàn)逐點提取子波,更精細(xì)地刻畫地震記錄的非穩(wěn)態(tài)特性。

ASST2的數(shù)學(xué)表達(dá)式[27]為:

(1)

同步擠壓變換在實現(xiàn)過程中可以加入閾值來控制譜系數(shù)的利用,在一定程度上可以降低原有時頻表示對噪聲的敏感程度。ASST2將自適應(yīng)S變換與同步擠壓變換相結(jié)合,通過改進(jìn)瞬時頻率和自適應(yīng)時窗,實現(xiàn)時頻分析結(jié)果能量聚焦并保證定位準(zhǔn)確[27]。相較于基于短時傅里葉變換和小波變換的同步擠壓變換時頻分析方法,ASST2具有更高的時頻分析精度,即使在有噪聲的情況下,也能夠有效地分辨出地震數(shù)據(jù)中頻率成分的變化。此外,該方法是數(shù)據(jù)驅(qū)動的,可根據(jù)信號自身特點自適應(yīng)調(diào)節(jié)窗口,有利于分析具有動態(tài)衰減特性的非穩(wěn)態(tài)地震數(shù)據(jù)。ASST2已經(jīng)成功地應(yīng)用于面波去除,說明了該方法的可行性和有效性[27]。

1.2 基于時頻譜提取時變子波

對于一個任意振幅譜w(f),統(tǒng)計參數(shù)(中心頻率fm和方差σ2)是用來定量描述振幅譜形狀及帶寬的重要參數(shù),其定義為:

(2)

為了更好地表示實際子波特征,HU等[28]構(gòu)建了加權(quán)指數(shù)型子波(FWE)表達(dá)式,該式對參數(shù)要求相對簡單,實際應(yīng)用相對靈活,波形豐富可廣泛用于多種形態(tài)的振幅譜[29]。加權(quán)指數(shù)型子波頻率域解析式如下:

(3)

式中:a是振幅歸一化系數(shù);n是控制對稱性的對稱指數(shù);f0是控制帶寬的特征頻率。將公式(3)代入公式(2),建立FWE公式的理論參數(shù)(n和f0)與振幅譜統(tǒng)計參數(shù)(fm和σ2)之間的解析關(guān)系:

(4)

由公式(4)可知,FWE中的理論參數(shù)(n和f0)可以由頻譜的中心頻率fm和方差σ2唯一確定,也即,對于任意頻譜,可以利用公式(4)和公式(3)進(jìn)行定量表征,再做反傅里葉變換得到最終時變子波。由于能量譜(即振幅譜的平方)對噪聲具有一定壓制作用,因此,公式(4)中通過能量譜得到的統(tǒng)計參數(shù)具有較好的抗噪性,進(jìn)而提高了子波估計的穩(wěn)定性。同時,公式(4)的推導(dǎo)過程沒有進(jìn)行近似,避免了誤差累計,最大程度地保證子波提取的精度。

根據(jù)傳統(tǒng)的褶積模型[30],地震記錄的矩陣形式可表示為地震子波矩陣W與反射系數(shù)向量的乘積,即:

s=Wr+n

(5)

式中:s,r和n分別代表了地震記錄、反射系數(shù)和噪聲的向量;W是由時不變地震子波w(t)組成的托普利茲矩陣。這一模型建立在子波穩(wěn)定的假設(shè)條件上,即假設(shè)地震子波在傳播過程中不會隨傳播距離的增加而變化。實際傳播過程中,由于地層吸收衰減效應(yīng),地震子波會隨著波的傳播不斷變化,因此,公式(5)不能準(zhǔn)確描述實際地震數(shù)據(jù)。

考慮實際地震數(shù)據(jù)的非穩(wěn)態(tài)特性,利用ASST2可以得到地震記錄每一時刻的頻率成分,即每一時刻的頻譜,根據(jù)公式(4)箭頭左側(cè)公式得到每一時刻頻譜的統(tǒng)計參數(shù)(fm和σ2),再利用公式(4)箭頭右側(cè)統(tǒng)計參數(shù)(fm和σ2)和加權(quán)指數(shù)型子波理論參數(shù)(f0和n)之間的解析關(guān)系計算f0和n,最終由公式(3)唯一確定該時刻子波表達(dá)式。將提取的時變子波沿對角線排列,并逐列取代W中的時不變子波w(t),建立時變子波矩陣W*。由此,可以將公式(5)中穩(wěn)態(tài)褶積模型改寫為基于時變矩陣的非穩(wěn)態(tài)地震記錄模型:

s=W*r+n

(6)

考慮到地震數(shù)據(jù)的非穩(wěn)態(tài)特性,本文提出的基于高精度時頻分析的時變子波提取方法的處理流程如圖1所示,根據(jù)公式(2)至公式(4)得到每一時刻子波的解析表達(dá)式,精確的推導(dǎo)過程為該方法的可行性提供了理論依據(jù),同時提高了時變子波提取的準(zhǔn)確性。值得注意的是,本文方法實現(xiàn)了從非穩(wěn)態(tài)地震記錄中直接提取時變子波,最大程度地避免了誤差累計,與常規(guī)方法相比,本文方法摒棄分段平穩(wěn)和子波振幅的假設(shè)條件,提高了子波提取方法的適用性,更有利于精細(xì)刻畫實際地震資料的時變特征。

圖1 基于ASST2時頻分析提取時變子波方法流程

2 模型數(shù)據(jù)測試

考慮地層吸收衰減因素的影響,建立如圖2a所示的非穩(wěn)態(tài)地震記錄模型,地震子波初始主頻為40Hz,隨著時間增加子波主頻不斷下降。圖2a中藍(lán)色實線為不含噪聲情況下的合成地震記錄(理論地震記錄),圖2b為對圖2a原始合成記錄利用ASST2得到的時頻譜圖。由圖可見,地震記錄的頻譜從淺至深逐漸變化,體現(xiàn)了地震記錄的非穩(wěn)態(tài)特性。同時,基于ASST2得到的高精度時頻分析結(jié)果能量聚焦性較好,有利于精確提取地震記錄的局部譜和實現(xiàn)時變子波的準(zhǔn)確估計。從圖2b中分別提取每一時刻的頻譜,并通過公式(2)至公式(4)得到每一時刻對應(yīng)的地震子波。圖2c至圖2f分別為沿著圖2b中白色虛線提取的189,192,450,800ms時刻時變子波振幅譜和時變子波(圖中紅色虛線所示)與相應(yīng)理論子波振幅譜和理論子波(圖中藍(lán)色實線所示)的對比結(jié)果。如圖中所示,利用本文方法提取的子波不僅在頻率域與理論振幅十分吻合,反變換到時域亦是如此,證明了方法的可行性和準(zhǔn)確性。對比圖2c和圖2d,上下相鄰地層之間雖然存在相互影響,但本文方法仍然可以準(zhǔn)確地提取每一時刻的子波,為后續(xù)高精度反演提供保障。隨著時間的增加,提取的子波主頻逐漸降低,頻帶變窄,相應(yīng)時域波形變寬,說明該方法真實地反映了地震子波在傳播過程中的動態(tài)衰減特性。圖2a 中的紅色虛線為根據(jù)公式(6)利用所提取的時變子波構(gòu)建時變子波矩陣W*與已知反射系數(shù)乘積重構(gòu)的地震記錄,與理論地震記錄吻合度較好。重構(gòu)地震記錄與理論地震記錄的相關(guān)系數(shù)c=0.99,進(jìn)一步驗證了本文方法所提子波的精確性。

為了測試該方法的抗噪聲能力,在上述非穩(wěn)態(tài)地震記錄模型中加入信噪比為8dB的隨機(jī)噪聲,測試結(jié)果如圖3所示。圖3a中藍(lán)色實線為含有噪聲的非穩(wěn)態(tài)地震記錄(理論地震記錄),圖3b為對圖3a理論地震記錄利用ASST2得到的時頻譜圖。對比圖2b和圖3b 的時頻分析結(jié)果可以發(fā)現(xiàn),由于ASST2算法具有一定的抗噪聲能力,因此即使在有噪聲的情況下,仍然可以得到高分辨率、高精度的時頻分析結(jié)果,進(jìn)而保證逐點提取子波的有效性。為了保證振幅譜統(tǒng)計參數(shù)(fm和σ2)計算的準(zhǔn)確性,在信噪比較低時,我們選取頻譜的有效頻帶范圍(0～100Hz)進(jìn)行計算。同樣,從圖3b中分別提取每一時刻的子波,圖3c 至圖3f分別為189,192,450,800ms時刻提取的時變子波振幅譜和時變子波(圖中紅色虛線所示)與理論子波振幅譜和理論子波(圖中藍(lán)色實線所示)對比結(jié)果。如圖中所示,在含有噪聲的情況下,利用本文方法所得時變子波和振幅都與理論子波基本吻合。對比圖2和圖3可知,所提取的時變子波并沒有因為噪聲的加入而發(fā)生劇烈變化,波形保持了很好的穩(wěn)定性,這是因為本文方法利用能量譜(即頻譜的平方)來估計時變子波的理論參數(shù),對噪聲有一定的壓制作用。此外,利用所提取的時變子波構(gòu)建時變子波矩陣與已知反射系數(shù)乘積重構(gòu)的地震記錄,如圖3a 中的紅色虛線所示,與理論地震記錄仍具有較好的一致性,其中,重構(gòu)地震記錄與理論地震記錄相關(guān)系數(shù)c=0.97。定性和定量分析發(fā)現(xiàn),即使是在含有中等噪聲的情況下,本文方法也可以提供準(zhǔn)確可靠的時變子波估計結(jié)果,由此表明該方法具有較好的抗噪能力和穩(wěn)定性。

圖2 基于ASST2時頻分析提取時變子波a 理論地震記錄(藍(lán)色實線)和重構(gòu)地震記錄(紅色虛線); b ASST2時頻譜; c 189ms時刻的子波振幅譜(上)和子波(下); d 192ms時刻的子波振幅譜(上)和子波(下); e 450ms時刻的子波振幅譜(上)和子波(下); f 800ms時刻的子波振幅譜(上)和子波(下)

圖3 含噪聲情況下基于ASST2時頻分析提取時變子波a 含噪聲的理論地震記錄(藍(lán)色實線)和重構(gòu)地震記錄(紅色虛線); b ASST2時頻譜; c 189ms時刻的子波振幅譜(上)和子波(下); d 192ms時刻的子波振幅譜(上)和子波(下); e 450ms時刻的子波振幅譜(上)和子波(下); f 800ms時刻的子波振幅譜(上)和子波(下)

為驗證本文方法的優(yōu)越性,利用圖2a所示地震記錄對基于自相關(guān)逐點提取子波的方法[17]進(jìn)行了處理測試,結(jié)果如圖4所示。圖4a中藍(lán)色實線為不含噪聲情況下的合成地震記錄(理論地震記錄),圖4b為利用短時傅里葉變換得到的時頻譜圖。對比圖2b和圖4b可以很直觀地感受到ASST2在對地震信號時頻分析中所具有的優(yōu)勢,其時頻譜圖品質(zhì)更好、能量聚焦性更好,由于圖2b中能量團(tuán)集中在較小的時間范圍內(nèi),有助于時變子波的精準(zhǔn)提取。圖4c至圖4f 分別為189,192,450,800ms時刻提取的時變子波振幅譜和時變子波(圖中紅色虛線所示)與理論子波振幅譜和理論子波(圖中藍(lán)色實線所示)的對比結(jié)果。如圖中所示,基于自相關(guān)法逐點提取的子波與理論子波存在較大的差異,且易受到相鄰反射系數(shù)的影響(圖4c和圖4d)。利用所提時變子波與反射系數(shù)重構(gòu)地震記錄的結(jié)果如圖4a紅色虛線所示,由于子波估計的不準(zhǔn)確導(dǎo)致合成地震記錄與理論值存在明顯的誤差,幅值也不能得到很好地恢復(fù),其中,重構(gòu)地震記錄與理論地震記錄相關(guān)系數(shù)c=0.87。因此利用該子波進(jìn)行后續(xù)反演處理,其結(jié)果必然不準(zhǔn)確。

圖4 基于自相關(guān)法提取時變子波a 理論地震記錄(藍(lán)色實線)和重構(gòu)地震記錄(紅色虛線); b 短時傅里葉時頻譜; c 189ms時刻的子波振幅譜(上)和子波(下); d 192ms時刻的子波振幅譜(上)和子波(下); e 450ms時刻的子波振幅譜(上)和子波(下); f 800ms時刻的子波振幅譜(上)和子波(下)

3 實際資料應(yīng)用

圖5顯示了實際地震資料單道處理結(jié)果(第11道),其中,采樣間隔為2ms,截取的時窗范圍為0.5～1.5s。分別應(yīng)用自相關(guān)法提取時不變子波、基于自相關(guān)法提取時變子波[17]以及本文方法提取時變子波。對于實際地震資料而言,無法直接從波形判斷子波提取的準(zhǔn)確性,為此根據(jù)所得子波構(gòu)建子波矩陣再利用L1范數(shù)稀疏約束正則化問題反演反射系數(shù)[31],進(jìn)一步對所提取的子波進(jìn)行評價。圖5a為第11道地震記錄,圖5b至圖5d中紅色虛線分別為基于自相關(guān)方法提取時不變子波、基于自相關(guān)方法提取的時變子波和本文方法所提子波的單道反演結(jié)果。對比基于3種方法提取子波的反演結(jié)果和利用測井?dāng)?shù)據(jù)計算得到的反射系數(shù)(圖5b至圖5d中藍(lán)色實線所示)可見,采用本文方法提取的時變子波的反演結(jié)果與基于測井?dāng)?shù)據(jù)計算的結(jié)果匹配性最好。分別計算反演結(jié)果與測井?dāng)?shù)據(jù)得到的反射系數(shù)的相關(guān)系數(shù),其中,采用時不變子波時,相關(guān)系數(shù)為0.39,采用自相關(guān)法提取的時變子波時,相關(guān)系數(shù)為0.42;而采用本文方法提取的時變子波時,相關(guān)系數(shù)為0.56。對實際資料的處理結(jié)果與理論結(jié)果一致,時變子波更符合地震信號的非穩(wěn)態(tài)特性。對比圖5中反射系數(shù)反演結(jié)果可知,子波提取的準(zhǔn)確性直接影響反射系數(shù)結(jié)果的精度,前兩種方法的反褶積結(jié)果存在由于子波的誤差引起的假反射系數(shù),利用本文方法處理后的噪聲較弱,同時相鄰反射系數(shù)更加清晰、層位置更容易分辨(圖5d箭頭所示)。

圖5 第11道地震記錄(a)和基于自相關(guān)方法提取時不變子波(b)、自相關(guān)方法提取時變子波(c)、本文方法提取子波(d)的反演結(jié)果(紅色虛線)以及基于測井資料計算的反射系數(shù)(藍(lán)色實線)

圖6為分別利用上述3種方法提取的子波矩陣?；谧韵嚓P(guān)法提取的時不變子波(圖6a)其波形不隨時間變化,且存在一定的旁瓣,基于自相關(guān)法提取的時變子波(圖6b)雖然其波形隨傳播時間逐漸變寬,符合地下傳播的動態(tài)衰減特性,但存在較多的旁瓣(圖中箭頭所示),會引起調(diào)諧效應(yīng)。與前兩種方法提取子波結(jié)果相比,本文方法得到的時變子波(圖6c)不僅可以反映地震資料的時變特征,同時其旁瓣能量較弱,避免了由于旁瓣能量過強(qiáng)而引起的假象。

圖6 自相關(guān)法所提取的時不變子波矩陣(a)、自相關(guān)法所提取的時變子波矩陣(b)和本文方法所提取的時變子波矩陣(c)

圖7為上述3種方法處理后的振幅譜,其中藍(lán)色實線為處理前地震記錄的振幅譜,紅色虛線為處理后的振幅譜。由圖7可知,應(yīng)用3種方法對實際地震資料處理后,振幅譜的有效頻帶均得到拓寬,而本文方法在保持地震資料信息的同時,實現(xiàn)了高頻、低頻信息同時恢復(fù)。根據(jù)已有研究[32],相對于高頻信息,低頻信息的恢復(fù)難度更大,而低頻信息對相對頻寬的貢獻(xiàn)較大,能較好地降低地震子波旁瓣效應(yīng),且在地震反演中起到非常重要的作用。圖8為分別利用上述3種方法提取子波后反演的二維地震剖面,道數(shù)為50,可見,與基于時不變子波處理后的剖面(圖8b)相比,基于時變子波的反褶積后的剖面質(zhì)量得到改善(圖8c和圖8d),分辨率較高?；跁r變子波的反演結(jié)果可以分辨原始數(shù)據(jù)中不可分辨的薄層,比如圖8a 中單個同相軸的反射(圖中箭頭所指位置),在基于本文方法的反演結(jié)果中呈現(xiàn)為兩個反射同相軸。對比圖8c和圖8d發(fā)現(xiàn),本文方法效果最佳,不僅增強(qiáng)了同相軸能量(圖中矩形框區(qū)域),而且呈現(xiàn)了更多的薄層細(xì)節(jié),同時提高了反演結(jié)果的分辨率,進(jìn)一步驗證了低頻信息在地震反演中所起的重要作用。定性和定量分析可知,相較于常規(guī)子波提取方法,本文方法具有較好的可行性和優(yōu)越性,能夠進(jìn)一步提高地震資料反映薄層真實細(xì)節(jié)的能力,有利于后續(xù)地震資料的精細(xì)解釋。

圖7 處理前、后地震記錄振幅譜a 基于時不變子波處理結(jié)果; b 基于自相關(guān)法時變子波處理結(jié)果; c 基于本文方法處理結(jié)果

圖8 實際地震剖面(a)、基于時不變子波的反演剖面(b)、基于自相關(guān)法提取時變子波的反演剖面(c)以及基于本文方法的反演剖面(d)

4 結(jié)論

針對地震資料的非穩(wěn)態(tài)特性,本文提出了基于二階自適應(yīng)同步擠壓S變換的時變子波提取方法,實現(xiàn)了從地震數(shù)據(jù)中直接提取時變子波,并建立了反映子波動態(tài)衰減特性的非穩(wěn)態(tài)地震數(shù)據(jù)褶積模型。相較于常規(guī)子波提取方法,本文方法通過推導(dǎo)局部譜的統(tǒng)計參數(shù)與加權(quán)指數(shù)型子波之間的數(shù)學(xué)關(guān)系,準(zhǔn)確得到時變子波的解析式,最大程度避免了誤差累計。模型數(shù)據(jù)測試結(jié)果證明了本文方法的可行性和有效性,即使在信噪比較低的情況下,也能夠保證子波提取的精度。實際資料應(yīng)用進(jìn)一步驗證了該方法能夠適應(yīng)地震資料的時變特性,準(zhǔn)確估計地震子波,有利于提高地震資料分辨薄層細(xì)節(jié)的能力和后續(xù)地震資料解釋的準(zhǔn)確性。

本文方法是基于子波相位時不變提出的,而實際地震資料中子波相位也會隨傳播時間變化,如何將時變子波相位提取與本文方法相結(jié)合是我們下一步的研究方向。