石戰(zhàn)戰(zhàn)，賀振華，文曉濤，唐湘蓉

（1.成都理工大學(xué)“油氣藏地質(zhì)及開發(fā)工程”國(guó)家重點(diǎn)實(shí)驗(yàn)室；2.成都理工大學(xué)信息工程學(xué)院）

一種基于EMD和GHT的儲(chǔ)層識(shí)別方法

石戰(zhàn)戰(zhàn)1，2，賀振華1，2，文曉濤1，2，唐湘蓉1，2

（1.成都理工大學(xué)“油氣藏地質(zhì)及開發(fā)工程”國(guó)家重點(diǎn)實(shí)驗(yàn)室；2.成都理工大學(xué)信息工程學(xué)院）

針對(duì)地震信號(hào)是具有一定帶寬的非平穩(wěn)信號(hào)，波形受這個(gè)帶寬內(nèi)所有頻率影響，提出結(jié)合經(jīng)驗(yàn)?zāi)B(tài)分解和廣義希爾伯特的儲(chǔ)層預(yù)測(cè)方法，并通過經(jīng)驗(yàn)?zāi)B(tài)分解將地震信號(hào)分解為頻率由高到低的一系列本征模函數(shù)，對(duì)分解后的每一個(gè)本征模函數(shù)分量作廣義希爾伯特計(jì)算。通過與已知井儲(chǔ)層發(fā)育情況對(duì)比，從中優(yōu)選出了對(duì)儲(chǔ)層響應(yīng)敏感的本征模函數(shù)分量進(jìn)行儲(chǔ)層預(yù)測(cè)。實(shí)際資料計(jì)算結(jié)果表明，結(jié)合經(jīng)驗(yàn)?zāi)B(tài)分解和廣義希爾伯特的儲(chǔ)層識(shí)別方法，具有較強(qiáng)的噪音壓制能力，能夠精確刻畫儲(chǔ)層形態(tài)。

經(jīng)驗(yàn)?zāi)B(tài)分解；本征模函數(shù)；廣義希爾伯特；非平穩(wěn)信號(hào)；儲(chǔ)層預(yù)測(cè)

0 引言

地震信號(hào)屬于非平穩(wěn)信號(hào)，其頻率成分隨時(shí)間變化[1］。傳統(tǒng)Fourier分析是建立在線性、平穩(wěn)信號(hào)的假設(shè)基礎(chǔ)之上，缺乏對(duì)信號(hào)時(shí)間和頻率的同時(shí)定位，不能有效分析信號(hào)局部特征[2］。Gabor[3］提出短時(shí)傅里葉變換（short time Fourier transform，STFT）或加窗傅里葉變換，STFT在整個(gè)頻率域上使用相同窗函數(shù)，不同頻率成分的分辨率相同。Gabor變換去掉了STFT窗函數(shù)必須是窄窗的要求，但其時(shí)間-頻率窗的高度和寬度不變，限制了其在非平穩(wěn)高頻和低頻信號(hào)上的應(yīng)用[4-5］。Vigner于 1932年提出Wigner分布（WD），Ville于 1948年提出 Wigner-Vile分布（WVD），它們都具有較好的時(shí)頻聚焦性，但對(duì)于多分量信號(hào)會(huì)產(chǎn)生嚴(yán)重的交叉干擾，使其譜分布難以解釋[4，6］。 Morlet于 1984 年提出具有多分辨率特性的小波變換，小波分析在頻域邊緣處能量泄露嚴(yán)重，小波基的選擇對(duì)于信號(hào)分析有很大影響[4-6］。Huang等于1998年提出經(jīng)驗(yàn)?zāi)B(tài)分解法（empirical modedecomposition，EMD），在機(jī)械振動(dòng)、結(jié)構(gòu)力學(xué)、天然地震等方面取得了初步應(yīng)用[7］。Flandrin和Gon?alvè[8］論述了 EMD 與小波分解的聯(lián)系與區(qū)別。Magrin-Chagnolleau和Baraniuk提出地震數(shù)據(jù)可用EMD 分析[1，9］。 筆者提出利用經(jīng)驗(yàn)?zāi)B(tài)分解法對(duì)信號(hào)進(jìn)行分解，再對(duì)分解后各分量進(jìn)行GHT計(jì)算。EMD分解的結(jié)果——本征模函數(shù)（IMF），它包含了信號(hào)從高到低不同頻率段的成分，是信號(hào)頻帶的一種自動(dòng)劃分，且隨信號(hào)本身的變化而變化，儲(chǔ)層信息也一定隱含在某些分量中[10］。這些隱含儲(chǔ)層信息分量的GHT則能更好地反映儲(chǔ)層的發(fā)育情況。

1 EMD方法原理

EMD方法分解思想是將信號(hào)分解為若干個(gè)相對(duì)平穩(wěn)的、互不相關(guān)的本征模態(tài)函數(shù)（intrinsic mode function，IMF）[11-15］，也稱為篩選過程（the sifting progress）[10］。 一個(gè)IMF分量必須具備以下 2個(gè)條件：①信號(hào)中的極值點(diǎn)（極大值和極小值點(diǎn)數(shù)之和）的數(shù)目必須和信號(hào)過0點(diǎn)的數(shù)目相等或最多只能相差1；②在信號(hào)的任意時(shí)間點(diǎn)上，信號(hào)的局部極大值確定的上包絡(luò)線和局部極小值定義的下包絡(luò)線平均值為0。

第一個(gè)限制條件近似于傳統(tǒng)的平穩(wěn)高斯過程關(guān)于窄帶的定義；第二個(gè)限制條件是一個(gè)創(chuàng)新想法，把傳統(tǒng)的全局限定變?yōu)榫植肯薅╗10-11，14-15］。

對(duì)于任意信號(hào)s（t），找出信號(hào)所有極值點(diǎn)，分別對(duì)極大值和極小值用三階樣條進(jìn)行插值，得到信號(hào)的上下包絡(luò)線 emax（t）和 emin（t）；求取上下包絡(luò)線的平均值，即可得到信號(hào)包絡(luò)的瞬時(shí)平均值m（t）

從原始信號(hào) s（t）中減去 m（t）得到兩者之差 h（t）

將 h（t）視為新的 s（t），重復(fù)以上步驟，直到 h（t）滿足 IMF 必備的 2 個(gè)條件時(shí)，則 h（t）為 s（t）的第一個(gè)IMF分量IMF1。

將視為新的 s（t），重復(fù)以上步驟，可以得到信號(hào) s（t）的第二個(gè)IMF2，第三個(gè)IMF3，……，如此重復(fù)直到最后一個(gè)序列rn不可再分解（極值點(diǎn)小于2個(gè)）時(shí)停止上述分解過程。此時(shí)，rn代表數(shù)據(jù)序列s（t）的某種趨勢(shì)。

IMF的2個(gè)限制條件只是一種理論上的要求，在實(shí)際篩選過程中，很難保證信號(hào)的局部均值絕對(duì)為0，篩選過程的停止準(zhǔn)則可以通過限制2個(gè)連續(xù)處理結(jié)果之間的標(biāo)準(zhǔn)差sd的大小來(lái)實(shí)現(xiàn)[11］

其中：sd為標(biāo)準(zhǔn)差；為求和符號(hào)；h（k-1）（t）和 h（k）（t）為信號(hào)減去包絡(luò)瞬時(shí)平均值得到新信號(hào)。

現(xiàn)以一個(gè)仿真信號(hào)為例說(shuō)明EMD分解。圖1a為原始仿真信號(hào)，圖1b、圖1c、圖1d、圖1e為該仿真信號(hào)的IMF1—IMF4分量，圖1f為原始信號(hào)殘余分量。從圖中可以看出，這4個(gè)IMF分量分別對(duì)應(yīng)信號(hào)從高到低的不同頻率段成分，而殘余分量代表信號(hào)的平均趨勢(shì)。若能從地震信號(hào)中分離出對(duì)油氣敏感的頻率成分，并對(duì)其進(jìn)行油氣檢測(cè)，檢測(cè)的成功率必然會(huì)提高[10］。

圖1 仿真信號(hào)的EMD分解Fig.1 Simulation signal and the EMD components

2 廣義希爾伯特原理

希爾伯特變換（Hilbert Transform，HT）是一種全通濾波器[16］，廣泛應(yīng)用于復(fù)數(shù)道分析，信號(hào) x（t）通過HT后，其負(fù)頻率成分作+90°相移，而正頻率作-90°相移。頻域HT表示如下

其中：X（ω）為輸入信號(hào) x（t）的 Fourier變換；hi（ω）為 x（t）頻域表示的 HT；sign（ω）為頻域符號(hào)函數(shù)；i為虛數(shù)單位。為了克服HT容易受噪音影響，羅毅于2003年提出加窗希爾伯特變換或廣義希爾伯特變換（Generalized Hilbert Transform，GHT）[17］。 信號(hào) x（t）的GHT用以下公式定義

現(xiàn)以正弦信號(hào)為例說(shuō)明GHT相對(duì)HT具有較好的去噪能力。圖2a為加噪音正弦信號(hào)；圖2b為通過HT計(jì)算獲得相位曲線；圖2c為通過GHT計(jì)算原始信號(hào)相位曲線；圖2d為筆者提出EMD分解獲得IMF3分量計(jì)算出的相位曲線。可以看出，通過GHT計(jì)算獲得的相位曲線優(yōu)于使用HT計(jì)算的結(jié)果；通過EMD分解后獲得的IMF3分量計(jì)算結(jié)果優(yōu)于直接使用GHT計(jì)算的結(jié)果；GHT所得到的曲線明顯光滑于HT所得到的結(jié)果。由此可見，GHT相對(duì)HT具有較好的去噪功能。地震信號(hào)通過GHT計(jì)算瞬時(shí)參數(shù)，預(yù)測(cè)結(jié)果精度必然提高。

圖2 仿真信號(hào)的GHT變換Fig.2 Simulation signal and generalized Hilbert transform

3 實(shí)例計(jì)算

為了驗(yàn)證該方法的有效性，筆者以我國(guó)某工區(qū)三疊系為例進(jìn)行計(jì)算。該區(qū)三疊系構(gòu)造受北東向斷裂帶和鹽體的塑性活動(dòng)控制，在擠壓應(yīng)力背景下形成北東向延伸的大型圈閉群。發(fā)育的儲(chǔ)層主要有三疊系上、中、下3個(gè)油組砂體[10］。筆者主要針對(duì)中油組進(jìn)行研究。

圖3a為原始地震數(shù)據(jù)過井剖面，well_1井在中油組鉆遇10 m油層，獲得高產(chǎn)工業(yè)油流，從剖面中雖然可以看出儲(chǔ)層位置（橢圓位置），但其振幅與周圍地層（圖中箭頭所示）差異不大，精確圈定儲(chǔ)層范圍有一定難度。圖3b、圖3c、圖3d為原始地震數(shù)據(jù)經(jīng)過計(jì)算EMD后獲得的IMF1，IMF2和IMF3分量數(shù)據(jù)體過井剖面。從圖中可以看出，IMF2分量中儲(chǔ)層邊界刻畫明顯（橢圓位置），儲(chǔ)層與周圍地層（圖中箭頭所示）振幅對(duì)比明顯；而IMF1（圖3a）和IMF3（圖3d）分量在相應(yīng)位置看不到儲(chǔ)層，說(shuō)明油氣信息主要集中在IMF2分量上。

對(duì)地震數(shù)據(jù)各分量計(jì)算GHT。圖4a為IMF2分量計(jì)算HT沿層切片，圖4b、圖4c、圖4d為各IMF分量計(jì)算GHT的沿層切片。從圖中可以看出，IMF2分量（圖4c）刻畫儲(chǔ)層邊界比其它分量清楚（圖示區(qū)域），比較 IMF2分量 HT計(jì)算結(jié)果（圖 4a）和 GHT計(jì)算結(jié)果（圖4c），可以看出IMF2分量計(jì)算GHT結(jié)果反映儲(chǔ)層分布范圍清晰。

圖3 實(shí)際地震資料EMD分解Fig.3 EMD of actual seismic data

圖4 各IMF分量GHT切片（歸一化）Fig.4 GHT slices of IMF components

4 結(jié)論

地震信號(hào)是具有一定帶寬的復(fù)合信號(hào)，通過EMD能夠分解出不同頻率段的分量，用井標(biāo)定后可優(yōu)選出對(duì)油氣響應(yīng)敏感的分量，而GHT具有較好的抑制噪聲能力，可提高儲(chǔ)層預(yù)測(cè)的成功率。但同時(shí)必須注意到EMD方法的端點(diǎn)效應(yīng)和GHT受窗函數(shù)影響較大，檢測(cè)結(jié)果必須結(jié)合測(cè)井、鉆井等其它資料，排除非地質(zhì)因素影響，提高解釋的可靠性和精度。

[1]Magrin-Chagnolleau I，Baraniuk R G.Empirical mode decomposition based time-frequency attributes[R］.69th Annual International Meeting， SEG，Expanded Abstracts，1999：1949-1952.

[2]劉麗娟，王山山.廣義S變換窗函數(shù)的分析和改進(jìn)[J］.巖性油氣藏，2007，19（2）：76-79.

[3]Gabor D.Theory of Communication[J］.J.IEEE，1946，93：429-457.

[4]張賢達(dá).現(xiàn)代信號(hào)處理：第二版[M］.北京：清華大學(xué)出版社，2003：368-395.

[5]劉喜斌，張寧，勾永峰，等.地震勘探信號(hào)時(shí)頻分析方法對(duì)比與應(yīng)用分析[J］.地球物理學(xué)進(jìn)展，2008，23（6）：743-753.

[6]張學(xué)濤，王祝文，原鏡海.利用時(shí)頻分析方法在陣列聲波測(cè)井中區(qū)分油水層[J］.巖性油氣藏，2008，20（1）：101-104.

[7]陳林，宋海斌.基于經(jīng)驗(yàn)?zāi)B(tài)分解的地震瞬時(shí)屬性提取[J］.地球物理學(xué)進(jìn)展，2008，23（4）：1 179-1 185.

[8]Flandrin P，Gon?alvès P.Empirical mode decompositions as datadriven wavelet-like expansions for stochastic processes[J］.InternationalJournalofWavelets，MultiresolutionandInformationProcesses，2004，2：477-496.

[9]Battista B M，Knapp C，McGee T，et al.Application of the empirical mode decomposition and Hilbert-Huang transform to seismic reflection data[J］.Geophysics，2007，72（2）：29-37.

[10]Wen X T，He Z H，Huang D J.Reservoir detection based on EMD and correlation dimension[J］.Applied Geophysics，2009，6（3）：70-76.

[11]何正嘉，訾艷陽(yáng)，張西寧.現(xiàn)代信號(hào)處理及工程應(yīng)用[M］.西安：西安交通大學(xué)出版社，2007：219-245.

[12]Huang N E，Shen Z，Long S R.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J］.Proceedings of the Royal Society.A，1998，454：903-995.

[13]李曉春，王祝文，岳崇旺，等.Hilbert-Huang變換在提取聲波測(cè)井信號(hào)儲(chǔ)集特性中的應(yīng)用[J］.巖性油氣藏，2009，21（1）：107-111.

[14]Flandrin P，Rilling G，Goncalves P.Empirical mode decomposition as a filter bank[J］.IEEE Signal Processing Letters，2004，11（2）：112-115.

[15]陳博濤，王祝文，丁陽(yáng)，等.Hilbert-Huang變換在陣列聲波測(cè)井信號(hào)時(shí)頻分析中的應(yīng)用[J］.巖性油氣藏，2010，22（1）：93-97.

[16]熊曉軍，賀振華，趙金明，等.一種基于GHT的裂縫檢測(cè)新方法[J］.石油地球物理勘探，2009，44（4）：442-444.

[17]Luo Y，Al-Dossary S，Marhoon M，et al.Acquisition，processing-Generalized Hilbert transform and its applications in Geophysics[J］.The Leading Edge，2003，22（3）：198-202.

[18]黨志敏，賀振華，黃德濟(jì).廣義希爾伯特變換在含噪信號(hào)邊緣檢測(cè)中的應(yīng)用效果分析[J］.油氣地球物理，2007，5（2）：25-28.

[19]賀振華，黃德濟(jì)，文曉濤.裂縫油氣藏地球物理預(yù)測(cè)[M］.成都：四川科技出版社，2007：138-141.

Reservoir detection based on EMD and GHT

SHI Zhan-zhan1，2， HE Zhen-hua1，2， WEN Xiao-tao1，2， TANG Xiang-rong1，2
（1.State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation， Chengdu University of Technology， Chengdu 610059， China；2.College of Information Engineering，Chengdu University of Technology，Chengdu 610059，China）

Seismic signal is a complexnon-stationarysignal with a definite frequency bandwidth,and the waveform is affected by all the frequency components in the band.So a new method for reservoir prediction combined empirical mode decomposition with generalized Hilbert transform is proposed.The empirical mode decomposition （EMD）technique adaptively decomposes the seismic signals into a series of intrinsic mode functions（IMF）,and every IMF component is carried out the GHT calculation.By comparing the results with reservoir development in fixed wells,the IMF components sensitive to the reservoir are optimized and used to carry out reservoir prediction.The computational result ofactual data shows that the method based on EMDand GHTcan be used todescribe the reservoir accurately.

empirical mode decomposition；intrinsic mode functions；generalized Hilbert transform；non-stationary signal；reservoir prediction

P315.01

A

1673-8926（2011）03-0102-04

2010-12-27；

2011-01-20

國(guó)家自然科學(xué)基金青年基金（編號(hào)：40904034）和國(guó)家自然科學(xué)基金“石油化工聯(lián)合基金”重點(diǎn)項(xiàng)目（編號(hào)：40839905）聯(lián)合資助。

石戰(zhàn)戰(zhàn)，1986年生，男，成都理工大學(xué)在讀碩士研究生，主要從事地震信號(hào)處理和儲(chǔ)層反演方面的研究工作。地址：（610059）四川省成都市成華區(qū)成都理工大學(xué)信息工程學(xué)院5301-2室。E-mail：shizhanzh@163.com

楊琦）