改進的EMD時頻分析地震資料高分辨率處理方法及其應(yīng)用研究

吳仕虎, 肖富森, 盧柃岐, 冉 崎, 張福宏, 楊廣廣

(1. 中國石油西南油氣田分公司 勘探開發(fā)研究院，成都 610041；2. 成都理工大學(xué)，成都 610059；3. 成都理工大學(xué) 數(shù)學(xué)地質(zhì)四川省重點實驗室，成都 610059)

0 引言

隨著全球油氣勘探的不斷發(fā)展，勘探難度也在不斷增大，獲得高保真、高分辨率的地震資料對復(fù)雜油氣藏的勘探開發(fā)具有重要意義。地震信號是由激發(fā)的地震子波經(jīng)過大地濾波后形成，由于地層介質(zhì)滯彈性和非均質(zhì)性，地震波在地層傳播過程中受大地濾波、吸收等影響，使地表接收到的地震記錄高頻成分損失嚴(yán)重[1-2]，導(dǎo)致地震資料的分辨率、信噪比降低。為了使地震資料能滿足儲層研究和油藏描述的要求，很多地球物理學(xué)家提出并發(fā)展了地震資料高分辨率處理的方法，常見的有譜白化、反Q濾波、時頻分析、多尺度聯(lián)合分析等幾類。

譜白化是通過展寬振幅譜來提高信號分辨率的方法，它不改變子波的相位譜，可以在時域中進行，也可以在頻域中進行。有多位學(xué)者進一步提出了小波譜白化和Hilbert譜白化方法，小波變換克服了傅里葉變換不能分析信號局部特征的缺陷[3-6]；Hilbert譜白化可在時、頻域內(nèi)同時增強信號的局部細(xì)節(jié)信息，使地震剖面更為清晰和連續(xù)，且具有更高的信噪比。反Q濾波是一種補償大地吸收衰減效應(yīng)的技術(shù)，它不僅可以補償振幅衰減和頻率損失，而且還可以改善記錄的相位特性,從而改善同相軸的連續(xù)性，提高弱反射波的能量和地震資料的信噪比、分辨率。早期的工作是依據(jù)Futterman數(shù)學(xué)模型提出的用級數(shù)展開作近似高頻補償?shù)姆碤濾波，進而提出振幅、相位同時進行穩(wěn)定補償?shù)姆碤濾波，并將其發(fā)展到Q值隨時間或深度連續(xù)變化的情況，近年來發(fā)展出利用Futterman衰減模型同時對振幅衰減和速度頻散進行校正，具有較高的穩(wěn)定性和計算精度[7-11]。時頻分析在時頻域上分析非平穩(wěn)信號，它克服了傅里葉分析時域與頻域分離的缺陷，較準(zhǔn)確地反映時間與頻率信息。目前發(fā)展了基于短時傅里葉變換的地層吸收補償技術(shù)、基于小波包分解的地層吸收補償方法、可調(diào)節(jié)窗函數(shù)標(biāo)準(zhǔn)差的非對稱廣義S變換，對地震記錄進行補償和基于改進的廣義S變換進行衰減補償[12-17]。多尺度聯(lián)合分析方法是利用測井、井間地震、VSP等地球物理技術(shù)，對地下同一目標(biāo)地質(zhì)體進行不同尺度性質(zhì)的反映，通過它們之間的聯(lián)合作用提高地震資料分辨率[18-21]。

上述方法對提高地震資料分辨率均有不錯的效果，然而這些方法或適用于平穩(wěn)信號、或需要假設(shè)條件、或需要提前獲取Q值，局限性大，工作量大，成本高。為了以更經(jīng)濟的方式得到真實可靠的高分辨率的地震資料，筆者提出了一種適用性廣、成本經(jīng)濟、合理有效地提高地震資料分辨率的技術(shù)方法，具有借鑒作用。本文方法綜合考慮信號在時域與頻域的特征，在時頻域上對信號進行高分辨率處理，首先利用EMD技術(shù)把地震信號自適應(yīng)分為有限個本征模態(tài)函數(shù)IMF；對每一個IMF分別進行高分辨率處理，然后疊加得到高分辨率的地震數(shù)據(jù)。理論資料與實際資料的處理結(jié)果表明，此方法具有較好的應(yīng)用效果，為解決高頻補償問題提供了一條全新的思路。

1 方法原理

1.1 EMD的方法原理

經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)是從信號本身出發(fā)，通過層層篩選，先將時間特征尺度小的高頻IMF(Intrinsic Mode Function, IMF)分量分離出來，然后將時間特征尺度大的低頻IMF分量分離出來，根據(jù)停止準(zhǔn)則，最后得到一個近似單調(diào)的殘余分量。為了去除疊加波的影響和使波形更加對稱，EMD的分解過程采用“篩分”算法。對于一個實信號x(t)，按照流程圖1的步驟進行分解，原信號x(t)被分解成n個IMFci(t)和一個信號余量rn(t)之和(式(1))。

(1)

圖1 EMD的分解步驟Fig.1 Steps of EMD decomposition

其中每個IMF必須滿足以下兩個條件：

1)在整個數(shù)據(jù)集中，該信號極值點的數(shù)目與零交叉點的數(shù)目相等或者最多相差一個點。

2)在任意時間上，該信號由局部最大值點形成的上包絡(luò)線和局部最小值形成的下包絡(luò)線的均值為零(即上包絡(luò)、下包絡(luò)局部是對稱的)。

IMF的第一個條件是為了滿足高斯正態(tài)平穩(wěn)過程的傳統(tǒng)窄帶要求；第二個條件把傳統(tǒng)的全局性限定為局域性，目的是避免模態(tài)混疊，保證一個固有模態(tài)函數(shù)僅僅包含一個基本模式的振蕩，去除因波形不對稱而造成的波動。這樣獲得的IMF信號在任何時刻都只具有單一的頻率成分[22]。

基于EMD分解的時頻分析是一種基于經(jīng)驗的分解，而不是基于數(shù)學(xué)表達(dá)式的分解，它直接從信號本身產(chǎn)生基函數(shù)，沒有交叉干擾項，因此具有很好的自適應(yīng)性、正交性和完備性。除此之外，EMD時頻分布還具有很好的時頻聚集性和高分辨率，時頻分布平面內(nèi)的能量分布十分密集，可以快捷地給出有效的時頻表示，突出時變信息。

1.2 地震記錄的振幅補償

因地層的吸收衰減作用，地震記錄的振幅譜有相對高頻的變化。設(shè)實際的地震記錄為S(t)，變換到頻域為：

S(f)=A(f)eiφ(f)

(2)

式中：A(f)為地震記錄的振幅譜;φ(f)為相位譜。采用平滑濾波的方法，可以由A(f)估算出地震記錄振幅譜的衰減曲線(即通常的地震子波的振幅譜)。該濾波器的平滑因子采用三角濾波算子(式(3))。

(3)

n=1,2,L…,2N-1(N為半濾波算子長)

式中:n為濾波算子樣點序號。用H(n)對地震記錄的振幅譜作平滑濾波，得到振幅譜的衰減曲線B(f)，如式(4)所示。

(4)

振幅補償曲線C(f)為衰減曲線的倒數(shù)，即C(f)=1/B(f)。但是在實際的計算中需對衰減曲線加入適當(dāng)?shù)陌自?，以保證振幅譜補償曲線的穩(wěn)定性。為了保證振幅補償?shù)牟怀^地震記錄的有效頻帶范圍，常對振幅補償曲線作限頻處理[23-24](式5)。

(5)

式中：fLC、fLP、fHP、fHC、fNq分別為地震信號的低截、低通、高通、高截、Nyquist頻率；T是與白噪系數(shù)U有關(guān)的穩(wěn)定系數(shù)，如式(6)所示。

(6)

為了從地震記錄恢復(fù)反射系數(shù)振幅譜，還應(yīng)對振幅譜補償結(jié)果進行線性加權(quán)處理，加權(quán)因子為式(7)。

f=fLC,fLC+Vf,L,fHC

(7)

式中：K為高頻增強系數(shù)，一般取大于1的數(shù)；Vf為頻率采樣間隔。

為了使本方法能保持相對振幅，需用輸入記錄的振幅譜予以標(biāo)定，以保證處理前、后的振幅譜的總能量不變。振幅補償?shù)淖罱K結(jié)果為式(8)中Y(f)。

Y(f)=A(f)·C(f)·

(8)

1.3 計算流程

筆者所述的地震信號高分辨率處理是以EMD方法為基礎(chǔ)，針對地震信號振幅特性而展開的地震信號高分辨率處理，其中所述的地震信號為用于勘探地下巖層性質(zhì)和形態(tài)的時域信號。具體的步驟：①獲取地震剖面，對一道地震信號采用EMD方法自適應(yīng)分解，得到固有模態(tài)函數(shù)IMF分量；②分析該地震信號、EMD分解得到的各個IMF分量的振幅譜，分別確定高分辨率處理的高截、高通、低通、低截閾值；③利用選擇的閾值分別對各個IMF分量進行高分辨率處理；④將所有時間范圍內(nèi)高分辨率處理后的IMF分量疊加重構(gòu)；⑤為了減少分解、重構(gòu)等過程中引起的誤差，可根據(jù)實際情況對重構(gòu)的信號進行平滑等修飾性處理；⑥重復(fù)步驟①～步驟⑤完成所有單道的處理，最終獲得高分辨率的二維地震剖面，用于地震解釋。同時需注意，不同區(qū)塊的地震資料差異較大，在選擇參數(shù)時應(yīng)根據(jù)實際的情況進行試驗后靈活處理，才可以最終達(dá)到最佳的處理效果。為了清晰地展示本算法，其流程圖如圖2所示。

圖2 基于EMD方法的高分辨率處理算法流程圖Fig.2 Flow chart of high-resolution processing algorithm based on EMD method

2 理論模型應(yīng)用效果分析

2.1 理論模型單道的效果分析

為了驗證該方法的可行性，設(shè)計了單道理論模型：用主頻為20 Hz的雷克子波(圖3(a))和6組數(shù)值為1.2的反射系數(shù)(圖3(b))進行褶積合成信號(圖3(c))，并展示出該合成信號的振幅譜(圖3(d))。在設(shè)計的過程中每組反射系數(shù)的時間距離遵循時間間距從0.121 s開始，以0.002為步長增加(表1)。

表1 每組反射系數(shù)之間的距離

圖3 理論模型的設(shè)計Fig.3 Design of theory model(a)20 Hz的雷克子波；(b)一系列反射系數(shù)；(c)合成信號；(d)合成信號的振幅譜

圖4 合成信號的IMF信號與IMF振幅譜Fig.4 IMF signal and IMF amplitude spectrum of synthesized signal(a)IMF的信號；(b)IMF的振幅譜

通過圖3(c)，可以清晰地看出,該合成信號第1個波形完全沒有分開，第2個波形幾乎沒有分開，從第3個波形開始，反射系數(shù)分開程度逐漸增大。對該合成信號用EMD方法分解得到5個IMF信號(圖4(a))，其對應(yīng)振幅譜見(圖4(b))。利用本文方法對上述合成信號進行提高分辨率處理，處理效果見圖5。由圖5(a)可以看出，本文EMD時頻分析高分辨率處理方法能夠有效地將一部分發(fā)生混疊的信號分開，從而使信號的分辨率得到提高。因該合成信號的采樣點數(shù)較少，筆者對處理前、后的振幅譜求包絡(luò)，得到合成記錄處理前后的頻譜包絡(luò)對比圖(圖5(b))。通過觀察高頻部分的振幅值有所提高，對應(yīng)說明能量得到一定程度補償。通過對比圖5可知，經(jīng)過本文方法的處理，該合成信號的主頻得到提高，高頻信號能量明顯得到提高，整個頻帶寬帶明顯拓寬。

2.2 理論剖面模型的效果分析

理論剖面實驗?zāi)Ｐ?，是?5個相同單道排列形成的剖面(圖6)。單道是由6個間距呈0.002 s遞增的反射系數(shù)(表1)與主頻為20 Hz的雷克子波褶積而成，采樣間隔為1 ms，采樣點數(shù)為1 000。

圖5 高分辨率處理前、后的效果對比Fig.5 Comparison of the effects before and after high resolution processing(a)處理前后的信號對比圖；(b)處理前后的振幅譜包絡(luò)對比圖

圖6 高分辨率處理前、后的地震剖面Fig.6 High-resolution seismic section before and after processing(a)原始地震剖面；(b)高分辨率處理后的地震剖面

從圖6(a)可以看出，在0 s～0.4 s之間，兩個波形是出現(xiàn)復(fù)合的現(xiàn)象，波形疊加在一起無法分辨，而且從下到上，波形的復(fù)合現(xiàn)象是越來越嚴(yán)重，兩個波形之間越來越難分辨。

理論剖面經(jīng)過本文方法處理之后，如圖6(b)所示，0.2 s～0.4 s之間的波形能夠辨認(rèn)，0.4 s～0.8 s之間的波形較之前分開明顯，分辨率明顯提高，0.8 s～1 s之間的兩個波形完全分開。

3 實際資料應(yīng)用效果分析

3.1 實際資料剖面的分析

為了進一步驗證該方法適用性，將該方法應(yīng)用于四川盆地川中地區(qū)高石梯井區(qū)的二維地震資料，該資料的目的層是燈影組，總共4 316道，每道采樣點為2 450個，采樣時間間隔為1 ms，實鉆井位于第2 576道。按照本文的方法原理進行高分辨率處理，首先采用EMD方法對單道地震數(shù)據(jù)進行分解，然后對分解之后得到的IMF分量進行高分辨率處理，重構(gòu)得到高分辨率的單道地震數(shù)據(jù)，最后逐道依次進行上述方法處理流程，重構(gòu)得到高分辨率的二維地震資料。因該資料的道數(shù)、采點數(shù)較多，現(xiàn)從原二維地震資料中截取部分形成一個新的地震剖面進行效果的比較分析，截取第2 300 道～3 000道，共701道，采樣時間2 s～2.45 s，共451個采樣點，在形成的新剖面中鉆井位于第277道。為了使呈現(xiàn)的圖像更加清晰、美觀，顯示了第240 道～300道，采樣時間0 s～0.4 s的剖面波形圖，處理前、后的對比如圖7所示。

圖7 二維實際資料處理前后對比圖Fig.7 Comparison of 2D actual seismic section before and after processin(a)實際資料的原始剖面；(b)實際資料的高分辨率處理后剖面

通過圖7(a)與圖7(b)對比可知，在0.03 s～0.05 s、0.07 s～0.09 s兩個時間區(qū)間分辨出薄互層反射信號(如圖7(a)和圖7(b)中紅色箭頭所示)。高石梯地區(qū)該套地層巖性為砂泥巖薄層互層，經(jīng)高分辨率處理后的地震剖面不僅保持了原始剖面的基本面貌，而且分辨率明顯提高、剖面的成層性更好，一些弱化的層位信息得到加強，缺失的相關(guān)層位信息得到了顯現(xiàn)。同時，可以清晰地看到，處理后的剖面上實鉆井的燈影組頂部界面更清晰，更容易區(qū)分，燈影組下部的第1 波峰～2 波峰與實鉆井的人工合成地震記錄匹配度較好，反射波同相軸的連續(xù)性得到明顯的改善，分辨率提高(如圖7(a)和圖7(b)中紅色的矩形框所示)。這說明本文方法能夠壓縮子波的長度，提高剖面的分辨率，有利于進一步的地質(zhì)體識別和解釋。

3.2 實際資料剖面的平均振幅譜分析

圖8為實際資料剖面處理前后的振幅譜對比，其中圖8(a)是處理前地震資料剖面的平均振幅譜，圖8(b)是處理后地震資料剖面的平均振幅譜。通過對比圖8(a)和圖8(b)可知，以主頻為中心的優(yōu)勢頻帶提高約12 Hz，處理前地震資料的主頻約為29 Hz，而處理后的資料主頻提高到41 Hz左右，頻帶寬度由處理前的17 Hz～41 Hz擴展到處理后的17 Hz～65 Hz。處理前、后的效果圖表明，本文方法能有效地提高地震資料的分辨率。

圖8 二維實際資料處理前后頻譜對比圖Fig.8 Comparison of 2D actual seismic data spectrum before and after processing(a)原始剖面振幅譜；(b)處理后剖面振幅譜

圖9 第277道處理前后的效果對比Fig.9 Comparison of the effects before and after processing(a)第277道處理前、后的信號；(b)第277道處理前、后的振幅譜

3.3 實際資料剖面的單道分析

抽取該二維地震資料的過井道，即第277道單道數(shù)據(jù)，用本文方法對其進行分析，處理前、后的效果如圖9所示。觀察圖9(a)可知，處理后的整個波形趨勢與原單道信號的波形基本保持一致，原來分開不明顯的波形，經(jīng)處理后波形有了明顯的分開，信號的分辨率有了較大地提高。圖9(b)顯示了高分辨率處理前后振幅譜的對比，通過觀察可知本文的方法提高了高頻部分，拓寬了頻譜的帶寬。在原始剖面中，第277道單道數(shù)據(jù)高頻振幅較弱，低頻成分占主導(dǎo)地位，處理后該道數(shù)據(jù)高頻信息顯著增強，分辨率明顯提高。

4 結(jié)論

針對大地濾波作用引起地震資料高頻成分損失，導(dǎo)致分辨率降低的問題，筆者利用EMD自適應(yīng)的時頻處理方法，提出基于EMD的高分辨率地震資料處理方法，介紹了實現(xiàn)該方法的關(guān)鍵思路及詳細(xì)過程。通過理論模型的測試和實際資料處理應(yīng)用，經(jīng)時域、頻域的綜合分析表明，該方法應(yīng)用效果良好，具有理論的可行性和結(jié)果實現(xiàn)的可靠性。該方法的優(yōu)勢還在于能夠保持原地震信號的主要特征不變的條件下，一些較弱或者缺失的層位信息得到了增強或者恢復(fù)，使地震資料的頻帶拓寬，高頻成分提高，有效地提高地震資料的分辨率，有助于更好地識別薄層、薄互層、微小構(gòu)造等，對目前越來越復(fù)雜的油氣勘探具有借鑒意義。