基于SBL-WVD的地震高分辨率時頻分析

紀永禎,張渝悅,朱立華,李 博

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.中國石油大學(北京),北京102249)

地震波在地下傳播時受吸收衰減等因素的影響,因而記錄到的地震數(shù)據(jù)通常具有非穩(wěn)態(tài)的特征,即地震數(shù)據(jù)的頻率成分隨時間變化而變化[1-2]。時頻分析方法可以將時域非穩(wěn)態(tài)地震數(shù)據(jù)變換到時頻域,利用變換獲得的頻率維度描述局部頻譜變化。不同的頻譜變化通常可以指示不同的地質(zhì)構(gòu)造和油氣儲層特征等[3-4]。因此,時頻分析方法是地震數(shù)據(jù)處理和解釋中常用的重要方法技術(shù)之一,被廣泛應(yīng)用于儲層描述和油氣檢測[5-6]、地質(zhì)構(gòu)造分析[7-8]、異常體識別[9-10]。

在過去幾十年間,時頻分析方法被廣泛研究并應(yīng)用于多種信號分析。在地震數(shù)據(jù)時頻分析中,最常用的工具包括短時傅里葉變換(STFT)、連續(xù)小波變換(CWT)、Gabor變換、S變換及其變體[11-14]。然而,由于窗口函數(shù)的影響,這些線性時頻分析工具通常無法兼顧時間分辨率和頻率分辨率,這種限制使這些工具在特定的情況下效果不佳[15]。魏格納威利分布(WVD)是一種雙線性的時頻分析方法,可以獲得較高分辨率的時頻分布結(jié)果[16]。然而,應(yīng)用于多頻率成分并存的地震數(shù)據(jù)時,WVD會受到交叉項干擾,且交叉項通常表現(xiàn)為高度震蕩,會嚴重影響時頻分布的精度和分辨率。一種減少甚至消除交叉項干擾的方法是利用匹配追蹤等算法將地震數(shù)據(jù)分解為獨立的子波分量,然后求取子波分量的時頻分布,最后重建地震數(shù)據(jù)的時頻分布[17-18]。此類方法有利于提高基于WVD時頻分析結(jié)果的分辨率,但其對地震數(shù)據(jù)的分解需要人工控制參數(shù)的介入,容易引起人為誤差,并且分辨率不足以識別薄層。

貝葉斯學習算法是TIPPING[19]提出的一種基于核函數(shù)的回歸分類算法,在多個領(lǐng)域得到了成功的應(yīng)用:WILLIAMS等[20]將稀疏貝葉斯學習算法用于視覺追蹤問題;SILVA等[21]運用貝葉斯學習算法進行了文本分類工作;THAYANANTHAN等[22]將模型推廣到多元輸出回歸問題;DEMIR等[23]將其運用于圖像處理中;NYEO等[24]將該算法用于動態(tài)光散射問題研究中;YUAN等[25]將其運用于疊后地震數(shù)據(jù)的反射系數(shù)反演;JI等[26]在頻率域疊前反演中初次引入該算法,獲得了較好的效果。在子波分解領(lǐng)域,貝葉斯學習算法較匹配追蹤、正交匹配追蹤等具有更多優(yōu)勢。匹配追蹤類算法屬于貪婪類算法,按順序選取構(gòu)建的子波庫中的原子,而貝葉斯學習算法可以通過自動相關(guān)判別選取更準確的分解子波表達式,同時可以避免“結(jié)構(gòu)性誤差”[27]。在含噪聲情況下,由于匹配追蹤類算法需要人為選定正則化參數(shù)(盡管可以參考信噪比),會引入人為誤差,而貝葉斯學習算法可以通過參數(shù)的自動學習確定最優(yōu)參數(shù),無需人為設(shè)定[28]。本文通過引入基于貝葉斯學習的子波分解算法來改善子波分解時的抗噪性、分辨率和準確度,并利用更接近于地震子波的雷克子波構(gòu)建子波庫,提出了一種兩步高分辨率時頻分析方法。在輸入端利用貝葉斯學習算法將獲得更好的子波分解結(jié)果,而基于此的WVD時頻分析方法也將具有更高的準確度和分辨率。文章首先介紹了時變雷克子波庫的構(gòu)建和非穩(wěn)態(tài)地震數(shù)據(jù)的正演模型,然后介紹了基于貝葉斯學習算法的子波分解方法,最后介紹了基于子波分解結(jié)果和WVD的地震數(shù)據(jù)時頻分析方法,通過與常規(guī)方法的比較,證實了本文方法在抗噪性、分辨率和準確度上的優(yōu)勢。

1 方法理論

1.1 構(gòu)建時變雷克子波庫

雷克子波是地震勘探中最常用的子波之一,采用不同主頻、相位和中心時差的雷克子波作為構(gòu)建子波庫的樣本。由于非穩(wěn)態(tài)地震數(shù)據(jù)具有時變特征,構(gòu)建子波庫時也需要考慮時變特點。首先,相同主頻和相位、不同中心時差的雷克子波可以寫成褶積矩陣的形式:

(1)

式中:W為(n+m-1)×n的矩陣,n代表子波采樣點數(shù),m代表反射系數(shù)采樣點數(shù);wn代表子波在第n個采樣點的值。將不同主頻和相位的雷克子波褶積矩陣組成如下子波庫:

(2)

1.2 基于時變雷克子波庫構(gòu)建頻率域地震褶積模型

若想將時變子波庫與非穩(wěn)態(tài)地震數(shù)據(jù)聯(lián)系起來,需要修改傳統(tǒng)的地震褶積模型,將反射系數(shù)的定義拓展為時變子波分解系數(shù),將非穩(wěn)態(tài)地震數(shù)據(jù)寫成由不同主頻、相位和中心時差的雷克子波組成的時變子波庫與時變子波分解系數(shù)的線性矩陣形式:

(3)

1.3 基于貝葉斯學習算法獲取時變子波分解系數(shù)的估計

將時變子波庫看成是稀疏核,采用貝葉斯學習算法獲取子波分解系數(shù)的估計。具體步驟如下:采用常用假設(shè),即非穩(wěn)態(tài)地震數(shù)據(jù)中含有噪聲n,符合零均值高斯分布,方差為σ2,那么非穩(wěn)態(tài)地震數(shù)據(jù)的似然函數(shù)可以寫成:

(4)

引入層狀高斯分布作為每一個待估計子波分解系數(shù)的先驗分布,對待估計子波分解系數(shù)進行稀疏約束[19]。先驗信息如下:

(5)

在貝葉斯框架下,待估計的子波分解系數(shù)的后驗概率可以寫成先驗信息和似然函數(shù)的乘積[15]:

(6)

式中:

(7)

(8)

式中:Q=σ2I+GH-1GH。將(8)式估計的h和μ代回(7)式獲得最終的參數(shù)估計結(jié)果。

由(8)式估計h,需采用基于最大期望算法的估計方法,該方法可以在算法進行中自適應(yīng)地改變估計時的中間參數(shù),有利于算法收斂并提高估計結(jié)果在含噪時的準確性。最大期望算法采用如下公式求解、估計h和σ2[17]。

(9)

其中,M是當前非零反射系數(shù)的個數(shù),h與待估計的子波分解系數(shù)的均值和方差有關(guān),σ2與數(shù)據(jù)匹配程度和稀疏度有關(guān)。至此,我們就能通過貝葉斯學習算法獲得子波分解系數(shù)的估計。

1.4 利用WVD獲取地震數(shù)據(jù)時頻分布

(10)

式中:Wseismic是時變地震數(shù)據(jù)的時頻譜;‖Φp‖是與第p個非零系數(shù)對應(yīng)的雷克子波的L2范數(shù);WVD(·)代表常規(guī)魏格納威利分布的計算結(jié)果。

2 應(yīng)用實例

2.1 模擬數(shù)據(jù)實驗

設(shè)計一個含噪聲模型,測試本文方法的子波分解和抗噪能力。選取5個雷克子波作為組成無噪數(shù)據(jù)的參考子波,然后加入與信號能量一致的隨機噪聲,即信噪比為1。其中雷克子波的主頻分別是50,30,50,40,30Hz;時間位置分別是50,65,90,105,150ms。與基于Gabor變換和匹配追蹤算法的時頻分析方法進行對比。圖1a給出了組成模擬數(shù)據(jù)的雷克子波。圖1b給出了加入噪聲的模擬數(shù)據(jù),其中,藍色曲線為加噪后的數(shù)據(jù),紅色曲線為無噪聲數(shù)據(jù),圖1b 是將圖1a中子波進行線性疊加獲得的。在圖1b 中,由于子波的時間位置較近,無噪聲數(shù)據(jù)展現(xiàn)了很多子波調(diào)諧效應(yīng),噪聲也使得模擬地震數(shù)據(jù)更加復(fù)雜。圖1c顯示了本文方法獲得的子波分解結(jié)果,與圖1a相比,子波分解結(jié)果非常準確。圖1d給出了利用匹配追蹤方法獲得的子波分解結(jié)果,由圖1d可見,其子波分解結(jié)果不理想,單個子波的估計結(jié)果與參考子波不同且分解子波的個數(shù)也多于參考子波個數(shù),產(chǎn)生了冗余信息。圖1e對比了無噪數(shù)據(jù)(藍色曲線)、本文方法重構(gòu)數(shù)據(jù)(黑色曲線)和匹配追蹤方法重構(gòu)數(shù)據(jù)(紅色曲線)。由圖1e可見,本文方法獲得的重構(gòu)數(shù)據(jù)與原始無噪數(shù)據(jù)更加匹配,雖然匹配追蹤算法的匹配程度也非常好,但其利用的子波并非完全準確,利用了更冗余的子波來重建原始數(shù)據(jù)的過程會導(dǎo)致后續(xù)時頻分布不準確和分辨率下降。圖2給出了本文方法、匹配追蹤方法和Gabor變換方法得到的時頻分布結(jié)果,圖中白色十字標記了參考的準確時頻分布位置。由圖2可見,基于Gabor變換方法得到的時頻分布結(jié)果分辨率較低(圖2c),其分辨率一定程度上取決于窗函數(shù)的選取,但由于無法同時兼顧時間和頻率分辨率,即使調(diào)節(jié)窗函數(shù),也很難獲得高分辨率的結(jié)果。對比本文方法(圖2a)和匹配追蹤算法(圖2b)得到的時頻分布結(jié)果可以看出,本文方法在有噪聲的情況下仍然較好地獲得了地震數(shù)據(jù)的時頻分布結(jié)果;而匹配追蹤算法獲得的時頻分布受到子波分解和噪聲的影響,時頻聚焦性和準確度與參考結(jié)果有一定差異。由此可見,本文方法的抗噪性能強,子波分解的準確度和分辨率高。

圖1 模擬數(shù)據(jù)及子波分解和重構(gòu)結(jié)果a 模擬數(shù)據(jù)所用的雷克子波; b 加噪聲的模擬地震數(shù)據(jù)(信噪比=1); c 基于本文方法獲得的子波分解結(jié)果; d 基于匹配追蹤算法獲得的子波分解結(jié)果; e 無噪聲地震數(shù)據(jù)與重構(gòu)數(shù)據(jù)的對比

圖2 3種方法得到的時頻分布結(jié)果a 本文方法; b 匹配追蹤方法; c Gabor變換方法

2.2 實際數(shù)據(jù)實驗

實際數(shù)據(jù)實驗中,首先選取井旁地震數(shù)據(jù)進行實驗對比,根據(jù)測井數(shù)據(jù)的解釋結(jié)果,在500ms附近有一含氣儲層,表現(xiàn)為高頻衰減的特征。圖3a、圖3b和圖3c分別給出了本文方法、匹配追蹤算法和Gabor變換方法獲得的井旁地震道時頻譜。對比圖3b和圖3c中400～500ms處紅色框內(nèi)的時頻譜可以看到,框中展現(xiàn)了兩個時頻聚焦的位置,體現(xiàn)了高頻衰減的特征。從圖3a可以清晰地看到,在紅色框中實際存在4個時頻聚焦的位置,在接近500ms時,時頻分布展現(xiàn)了幅度最大的從高頻向低頻移動的特征?；趫D3b和3c判斷,含氣儲層存在于400～500ms處,但從圖3a可以推斷,含氣儲層應(yīng)該出現(xiàn)在480～500ms位置。由此可見,本文方法有效提高了時頻分析的分辨率和精度。另外,在時間深度800ms以下的部分,本文方法也體現(xiàn)了更高的分辨率,一些小層展現(xiàn)得比較好。

圖3 3種方法獲得的井旁地震道時頻譜a 本文方法; b 匹配追蹤方法; c Gabor變換方法

圖4a、圖4b和圖4c分別給出了實際地震數(shù)據(jù)采用本文方法、匹配追蹤算法和Gabor變換方法獲得的30Hz單頻切片。對比3種方法得到的單頻切片可見,基于子波分解類的時頻分析方法較基于時窗的Gabor方法分辨率高。對比圖4a和圖4b可以看出,圖4b中的結(jié)果出現(xiàn)了比較明顯的“掛面條”現(xiàn)象,受噪聲影響,各道數(shù)據(jù)的時頻分析結(jié)果不夠穩(wěn)定,且分辨率比圖4a低,時頻聚焦性也略遜于圖4a。實際數(shù)據(jù)測試結(jié)果證明了采用本文方法得到的時頻分析結(jié)果在分辨率、精度和抗噪性上具有一定優(yōu)勢,有利于提高后續(xù)基于時頻分析結(jié)果的處理和解釋成果的質(zhì)量。然而本文方法也具有一定的局限性,從圖4a可以看出,本文方法獲得的時頻分布的連續(xù)性仍有提高的空間。

圖4 3種方法獲得的實際數(shù)據(jù)30Hz單頻切片a 本文方法; b 匹配追蹤方法; c Gabor變換方法

3 結(jié)論

本文提出了一種基于貝葉斯學習算法和魏格納威利分布的兩步高分辨率時頻分析方法。首先利用貝葉斯學習算法的高精度信號分解能力,獲得地震數(shù)據(jù)的子波分解結(jié)果,同時在一定程度上壓制了隨機噪聲,再利用魏格納威利分布的高分辨率時頻分析能力獲得地震數(shù)據(jù)的時頻分布結(jié)果。模擬數(shù)據(jù)和實際地震數(shù)據(jù)測試結(jié)果表明,本文方法與常規(guī)基于Gabor變換和匹配追蹤的時頻分析方法相比,具有更高的精度和分辨率。需要注意的是,本文方法采用的子波庫為雷克子波庫,在地震子波與雷克子波差距較大時,該方法獲得的時頻分析精度將受到一定影響。下一步可以通過改進構(gòu)建的雷克子波庫,提高方法的適用性。