一種基于EMD和GHT的儲層識別方法

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

（1.成都理工大學“油氣藏地質及開發(fā)工程”國家重點實驗室；2.成都理工大學信息工程學院）

一種基于EMD和GHT的儲層識別方法

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

（1.成都理工大學“油氣藏地質及開發(fā)工程”國家重點實驗室；2.成都理工大學信息工程學院）

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

經驗模態(tài)分解；本征模函數(shù)；廣義希爾伯特；非平穩(wěn)信號；儲層預測

0 引言

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

1 EMD方法原理

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

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

對于任意信號s（t），找出信號所有極值點，分別對極大值和極小值用三階樣條進行插值，得到信號的上下包絡線 emax（t）和 emin（t）；求取上下包絡線的平均值，即可得到信號包絡的瞬時平均值m（t）

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

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

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

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

其中：sd為標準差；為求和符號；h（k-1）（t）和 h（k）（t）為信號減去包絡瞬時平均值得到新信號。

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

圖1 仿真信號的EMD分解Fig.1 Simulation signal and the EMD components

2 廣義希爾伯特原理

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

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

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

圖2 仿真信號的GHT變換Fig.2 Simulation signal and generalized Hilbert transform

3 實例計算

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

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

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

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

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

4 結論

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

[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］.巖性油氣藏，2007，19（2）：76-79.

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

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

[5]劉喜斌，張寧，勾永峰，等.地震勘探信號時頻分析方法對比與應用分析[J］.地球物理學進展，2008，23（6）：743-753.

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

[7]陳林，宋海斌.基于經驗模態(tài)分解的地震瞬時屬性提取[J］.地球物理學進展，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]何正嘉，訾艷陽，張西寧.現(xiàn)代信號處理及工程應用[M］.西安：西安交通大學出版社，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變換在提取聲波測井信號儲集特性中的應用[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]陳博濤，王祝文，丁陽，等.Hilbert-Huang變換在陣列聲波測井信號時頻分析中的應用[J］.巖性油氣藏，2010，22（1）：93-97.

[16]熊曉軍，賀振華，趙金明，等.一種基于GHT的裂縫檢測新方法[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］.油氣地球物理，2007，5（2）：25-28.

[19]賀振華，黃德濟，文曉濤.裂縫油氣藏地球物理預測[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

國家自然科學基金青年基金（編號：40904034）和國家自然科學基金“石油化工聯(lián)合基金”重點項目（編號：40839905）聯(lián)合資助。

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

楊琦）