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

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

0 引言

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

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

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

1 方法原理

1.1 EMD的方法原理

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

(1)

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

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

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

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

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

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

1.2 地震記錄的振幅補(bǔ)償

因地層的吸收衰減作用，地震記錄的振幅譜有相對高頻的變化。設(shè)實(shí)際的地震記錄為S(t)，變換到頻域?yàn)椋?/p>

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

(2)

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

(3)

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

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

(4)

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

(5)

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

(6)

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

f=fLC,fLC+Vf,L,fHC

(7)

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

為了使本方法能保持相對振幅，需用輸入記錄的振幅譜予以標(biāo)定，以保證處理前、后的振幅譜的總能量不變。振幅補(bǔ)償?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分量進(jìn)行高分辨率處理；④將所有時間范圍內(nèi)高分辨率處理后的IMF分量疊加重構(gòu)；⑤為了減少分解、重構(gòu)等過程中引起的誤差，可根據(jù)實(shí)際情況對重構(gòu)的信號進(jìn)行平滑等修飾性處理；⑥重復(fù)步驟①～步驟⑤完成所有單道的處理，最終獲得高分辨率的二維地震剖面，用于地震解釋。同時需注意，不同區(qū)塊的地震資料差異較大，在選擇參數(shù)時應(yīng)根據(jù)實(shí)際的情況進(jìn)行試驗(yàn)后靈活處理，才可以最終達(dá)到最佳的處理效果。為了清晰地展示本算法，其流程圖如圖2所示。

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

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

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

為了驗(yàn)證該方法的可行性，設(shè)計了單道理論模型：用主頻為20 Hz的雷克子波(圖3(a))和6組數(shù)值為1.2的反射系數(shù)(圖3(b))進(jìn)行褶積合成信號(圖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))。利用本文方法對上述合成信號進(jìn)行提高分辨率處理，處理效果見圖5。由圖5(a)可以看出，本文EMD時頻分析高分辨率處理方法能夠有效地將一部分發(fā)生混疊的信號分開，從而使信號的分辨率得到提高。因該合成信號的采樣點(diǎn)數(shù)較少，筆者對處理前、后的振幅譜求包絡(luò)，得到合成記錄處理前后的頻譜包絡(luò)對比圖(圖5(b))。通過觀察高頻部分的振幅值有所提高，對應(yīng)說明能量得到一定程度補(bǔ)償。通過對比圖5可知，經(jīng)過本文方法的處理，該合成信號的主頻得到提高，高頻信號能量明顯得到提高，整個頻帶寬帶明顯拓寬。

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

理論剖面實(shí)驗(yàn)?zāi)Ｐ?，是?5個相同單道排列形成的剖面(圖6)。單道是由6個間距呈0.002 s遞增的反射系數(shù)(表1)與主頻為20 Hz的雷克子波褶積而成，采樣間隔為1 ms，采樣點(diǎn)數(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 實(shí)際資料應(yīng)用效果分析

3.1 實(shí)際資料剖面的分析

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

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

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

3.2 實(shí)際資料剖面的平均振幅譜分析

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

圖8 二維實(shí)際資料處理前后頻譜對比圖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 實(shí)際資料剖面的單道分析

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

4 結(jié)論

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