2016, Vol. 31
时频分析方法(Cohen et al., 1995)广泛应用于地震沉积学与地震储层描述中, 尤其在层序地层旋回分析、岩性识别、频谱调谐确定储层厚度(陈学华等, 2009)、地层结构成像、储层参数预测及油气检测等中应用较多.目前在地球物理应用中, 常用的时频分析方法有四种：
1) 短时傅里叶变换(STFT);
2) 小波变换(CWT);
3) S变换(ST);
4) 瞬时频率.
由这四种方法的原理出发, 可以将时频分析方法分为二类, 第一类为时窗性质的时频分析方法, 例如短时傅里叶变换, 小波变化(Morlet et al., 1982; 高静怀等, 1997), S变换(Stockwell, 2007)等, 该类方法目前应用效果不错, 但主要缺陷有两点：第一, 该类算法用时间窗内的数学计算来代替样点的瞬时频率, 是一种近似.第二, 由于加窗计算, 其分辨率受海森伯格测不准准则限制, 即时间分辨率高, 频率分辨率就低, 频率分辨率高, 时间分辨率就低, 时频分辨率不可能同时达到高精度要求, 如图 1中的左图所示, 对于一个单周期的正弦信号, 通过S变换和小波变换, 都不能够达到所期望的时频局部性的要求.第二类为如1979年Taner提出的复数道瞬时频率方法(Taner et al., 1979), 该方法从频率计算的数学定义出发, 并利用希尔伯特变换构建复数道, 使得求解出的频率为正值, 从而其求取结果具有实际的物理意义, 是真频率, 但其有一个理论前提, 输入信号必须为平稳信号.如图 1中的右图所示, 考虑一个简单信号s(t)=β+cos(t), β=0, 0.5, 5, 很显然, 该信号的瞬时频率应为1/2π, 利用常规希尔伯特变换构建复数道求取瞬时频率, 可以发现, 只有在β=0时, 利用该方法求取的频率才是正确值1/2π, 而当β=0.5, 5时, 其结果表现为震荡状的畸变形态.可以预见, 对于一个简单信号, 直接应用希尔伯特变换求取频率都会出现错误, 那么若输入信号为复杂的地震信号, 其结果的准确性是值得商榷和怀疑的.
|
图 1 两类传统时频分析方法的缺陷 Figure 1 Drawbacks of two conventional time frequency analysis method |
希尔伯特黄变换(HHT)的出现为克服前述两类方法的缺陷提供了一种解决方案, 希尔伯特黄变换最初由美国工程院院士诺顿.黄提出(Huang et al., 1998), 随后有多位国内外学者对其进行了理论研究(Flandrin et al., 2004; Patrick et al., 2004; Peng et al., 2005; Wu et al., 2009), 在石油勘探领域, 也有部分学者在实际应用上进行了相应探索(段生全等, 2005; Battista et al., 2007; 皮红梅等, 2007; 葸晓宇等, 2007; 杨培杰等, 2007; 黄光南等, 2009; 刘庆敏等, 2009).希尔伯特黄变换实际上由两个子变换组成(如图 2), 即希尔伯特变换和黄变换, 其中黄变换又被称作经验模态分解(EMD).
|
图 2 希尔伯特黄变换 Figure 2 Hilbert-Huang transform |
希尔伯特变换(Huang et al., 1999)后的瞬时频率是具有实际物理意义的真实频率, 但其理论前提是输入信号为平稳信号.地震信号为典型的非平稳信号, 那么直接对其应用希尔伯特变换, 其结果显然是不准确的, 而经验模态分解过程则恰好是将一个非平稳信号转化为一系列平稳信号的过程, 其中的每个平稳信号, 称作为本征模态函数(IMF).假设输入信号为xori(t), 则经验模态分解(EMD)可分为两步完成.
第一步：EMD(xori(t))
(1) 原始函数xori(t), x(t)=xori(t).
(2) IMF1=Sift(x(t)), x(t)=x(t)－IMF1.
(3) 重复步骤(1) , (2) k次, 直至满足条件.
第二步：Sift(x(t))
(1) 确定函数x(t)的所有极值点, 极大值序列包络为e1(t), 极小值序列包络为e2(t).
(2) 求取均值包络
(3) h1(t)=x(t)－m1(t), 如果h1(t)满足固有模态函数所需条件, 则h1(t)即为第一阶固有模态, 否则将h1(t)作为新的输入函数, 重复步骤(1) , (2) .直至得到满足固有模态函数条件的hk(t).
由上述流程可以看出, 经验模态分解(EMD)是一种自适应变换, 但正是因为其自适应性, 缺乏约束条件, 使其不可避免的存在缺陷.这种缺陷称为“模态混叠”(Wu et al., 2009).EMD分解算法假设一个信号由多个模态组成, 一个模态描述一个单一的震动状态.这种假设条件类似于傅里叶变换, 假设信号由多个三角函数组成.而在经验模态分解过程中, 期望将这些单一的模态干净地分离出来, 传统经验模态分解方法(EMD)由于算法本身的局限性, 常常会在分离出来的平稳信号中包含多个模态, 造成模态混叠, 其结果会导致随后频谱计算错误.国内外学者对模态混叠缺陷提出了一些改进的方法.如禹丹江等采用FFT结合EMD(禹丹江和任伟新, 2005), 将信号先进行频谱分析, 根据频谱分析结果对EMD过程设置分解中每一次分解的特征时间尺度上限.该方法的主要局限是FFT本身是不适合分析非线性非平稳信号的.Wu和Huang提出了一种EEMD方法(Wu et al., 2009), Torres提出了一种CEEMD方法(Torres等, 2011), 这两种方法都利用高斯白噪声迭代求取各阶模态函数, EEMD方法在每次迭代时引入噪声, 在部分解决模态混叠现象的同时, 引入其他本不属于信号本身的模态, 从而干扰了信号模态的正确分解, CEEMD保证了信号分解的完备性, 可以较好地解决模态混叠效应.
CEEMD方法
CEEMD(完备集合EMD)算法流程(Torres et al., 2011)：
CEEMD(xori(t))
(1) 求取第一个固有模态函数, IMF1=
(2) 计算剩余量, r1(t)=xori(t)－IMF1.
(3) 剩余IMFS的计算采用如下公式：IMF(k+1) =
完备集合经验模态分解(CEEMD)算法基于传统的经验模态(EMD)算法, 其中一阶模态(IMF1)的求取与集合经验模态分解(EEMD)一致, 剩余阶模态采用对噪声求取EMD分解的第一项来完成, 完备集合经验模态分解(CEEMD)算法在输入信号上引入高斯白噪声, 通过高斯白噪声的扰动效应, 来达到消除模态混叠的作用, 同时其分解过程满足完备性, 不会引入新的非信号本身的模态.
图 3为CEEMD方法消除模态混叠效果展示, 其中输入信号采用合成数据, 该信号以20 Hz余弦信号作为背景, 在2 ms 至401 ms处线性叠加7 Hz, 30 Hz, 40 Hz余弦信号, 在1070 ms、1100 ms处叠加30 Hz雷克子波, 在1300 ms出叠加100 Hz雷克子波.图 3中左图为传统EMD分解效果, 右图为CEEMD分解效果, 从图中可以看出, EMD分解存在大量模态混叠现象, 可以看出IMF1分量未能单独分解出100 Hz的雷克子波, 而是伴随有许多低频分量, IMF2、IMF3亦存在高频低频混合的现象.显然, 得到的结果不利于信号的时频分析.CEEMD分解其模态混叠现象进一步得到缓解, 100 Hz雷克子波完全分解, IMF2, IMF3主要反映30 Hz雷克子波, IMF5主要反映20 Hz余弦背景信号.
|
图 3 CEEMD与EMD模态混叠消除效果对比 Figure 3 Mode mixing elimination comparison between CEEMD and EMD |
假设一个简单的线性叠加信号s(t)=cos4πt+cos40πt+cos200πt, 该信号由频率为2 Hz, 20 Hz, 100 Hz的三个余弦信号线性叠加而成.图 4为分别利用基于CEEMD的希尔伯特黄变换(HHT), 小波变换(CWT), 短时傅里叶变换(STFT)的得到的时频谱.从图中可以看出, STFT方法不能将这三个信号完全区分开, CWT和HHT均能有效区分三个频率, 但HHT纤细谱线代表了更高的分辨率.图 5为另外一个合成数据的例子, 分析信号为一频率阶跃信号s(t)=
|
图 4 频率为2, 20, 100的余弦线性叠加信号的时频谱对比图 Figure 4 Comparison of time frequency spectrum of linear stack cosine signal at 2 Hz, 20 Hz, 100 Hz |
|
图 5 频率阶跃信号的时频谱对比图 Figure 5 Comparison of time frequency spectrum of frequency step signal |
采用两种时频分析方法对图 6是所示的地震道进行时频分析, 图 6中和图 6右分布为小波变换和基于CEEMD的希尔伯特黄变换的结果.从图中可以看出, 小波变换和基于CEEMD的黄变换的分频结果基本形态上保持一致, 但CEEMD的线状谱代表了更高的分辨率, 小波变换的结果更趋于CEEMD后平滑后的效果.图 7是分别利用基于CEEMD的黄变换, 短时傅里叶变换, 小波变换对一分频20 Hz的时间切片进行时频分析的结果.可以看出, CEEMD的分频结果具有高分辨率的优势.
|
图 6 实际地震道HHT与CWT时频谱对比 Figure 6 Time frequency spectrum comparison between HHT method and CWT method for real seismic trace |
|
图 7 三种分频算法20 Hz平面图对比 Figure 7 Slice comparison of three time frequency methods at 20 Hz |
图 8为一个基于CEEMD的HHT变换与S变换在油气检测上的例子.该数据为准噶尔盆地的碳酸盐与碎屑岩过渡带, 其中流体类型丰富, 油、气、水及干层均存在.通过对该工区多口井的频谱分析, 发现20 Hz能够代表该工区的油层响应, 因此对比显示了20 Hz的S变换与20 Hz的CEEMD希尔伯特黄变换的联井剖面图.图中共有4口井, 其中对比显示了3口, 从井段上的油层干层上看, 其中的红圈代表油层, 绿圈代表干层, 可以看出, 无论是S变换还是基于CEEMD的HHT变换, 都能够反应匹配油层、干层信息, 从HHT变换的峰值频率剖面图来看, 在油层为低频响应, 干层为高频响应, 与井上的显示一致, 且通过井上分层对比发现, HHT变换能更精确地描述油层、干层的位置, 也就是说HHT的不连续行代表了更高的分辨率.
|
图 8 CEEMD希尔伯特黄变换与S变换联井剖面油气检测比较 Figure 8 Oil&Gas detection comparison between Hilbert-Huang transform based on CEEMD and S transform at coupled well profile |
图 9是图 8同一剖面下一口单井的放大显示图, 对比的是基于EMD的HHT变换与CEEMD的HHT变换的20 Hz单频图.图中的黄圈代表了井上的两个层段, 其中o1_t, o1_b代表了油层, d1_t, d1_b代表了干层, 可以看出, 基于CEEMD的HHT变换与基于EMD的HHT变换在干层上表现一致, 但基于CEEMD的HHT变换在油层有着更加明显的高振幅显示.由于CEEMD能够克服大部分的“模态混叠”效应, 其在油气检测上比基于EMD的HHT变换具有更高的准确度.
|
图 9 基于EMD的HHT变换与基于CEEMD的HHT变换在油气检测中的对比 Figure 9 Oil and Gas detection comparison between HHT based on CEEMD and HHT based on EMD |
4.1 合成模型数据以及实际数据的应用效果表明, HHT时频分析方法能够正确反映信号的时频特性, 与小波变换(CWT), 短时傅里叶变换(STFT), S变换(ST)等常规时频方法相比较而言, HHT时频分析方法具有更高的分辨率.
4.2 CEEMD算法克服了原EMD算法中大量存在的“模态混叠”现象, 使得HHT时频分析方法提供更加精确的时频分辨能力.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [1] | Battista B M, Knapp C, McGee T, et al.2007. Application of the empirical mode decomposition and Hilbert-Huang transform to seismic reflection data[J]. Geophysics, 72 (2) : H29–H37. DOI:10.1190/1.2437700 |
| [2] | Chen X H, He Z H, Huang D J, et al.2009. Low frequency shadow detection of gas reservoirs in time-frequency domain[J]. Chinese J. Geophys. (in Chinese), 52 (1) : 215–221. |
| [3] | Cohen L.1995. Time-frequency Analysis[M]. Englewood Cliffs NJ: Prentice-Hall . |
| [4] | Duan shengquan, He Zhenhua, Huang Deji.32. Application of the Hilbert-Huang transform to the analysis of seismic signal[J]. Journal of Chengdu University of Technology (Science & Technology Edition)(in Chinese), 4 (396) : 400. |
| [5] | Flandrin P, Rilling G, Goncalves P.2004. Empirical mode decomposition as a filter bank[J]. IEEE Signal Processing Letters, 11 (2) : 112–114. DOI:10.1109/LSP.2003.821662 |
| [6] | Gao Jinghuai, Wange Wenbing, Zhu Guangming.1997. Wavelet Transform and instantaneous attributes analysis[J]. Chinese J. Geophys. (in Chinese), 40 (6) : 821–832. |
| [7] | Huang Guangnan. 2009. Hilbert Huang transform and applications of it in seismic data analysis and processing[D]. Qingdao: Ocean University of China. |
| [8] | Huang N E, Shen Z, Long S R, et al.1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analyses[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 454 (1971) : 903–995. DOI:10.1098/rspa.1998.0193 |
| [9] | Huang N E, Shen Z, Long S R.1999. A new view of nonlinear water waves: The Hilbert spectrum[J]. Annual Review of Fluid Mechanics, 31 : 417–457. DOI:10.1146/annurev.fluid.31.1.417 |
| [10] | Huang N E, Wu Z H, Long S R, et al.2009. On instantaneous frequency[J]. Advances in Adaptive Data Analysis, 1 (2) : 177–229. DOI:10.1142/S1793536909000096 |
| [11] | Liu Qingming, Yang Wuyang, Li Shuping, et al.2009. Application of the Hilbert-Huang transform to seismic data processing[J]. Northwestern Seismological Journal (in Chinese), 31 (3) : 211–216. |
| [12] | Morlet J, Arens G, Fourgeau E, et al.1982. Wave propagation and sampling theory—Part II: : sampling theory and complex waves[J]. Geophysics, 47 (2) : 222–236. DOI:10.1190/1.1441329 |
| [13] | Patrick F, Rilling G. 2004. Detrending and denoising with empirical mode decomposition[C].// EUSIPCO, 12th European Signal Processing Conference. |
| [14] | Peng Z K, Tes P W, Chu F L.2005. An improved Hilbert-Huang transform and its application in vibration signal analysis[J]. Journal of Sound and Vibration, 286 (1-2) : 187–205. DOI:10.1016/j.jsv.2004.10.005 |
| [15] | Pi Hongmei, Liu Cai, Wang Dian.2007. Using Hilbert-Huang transform to pick up instantaneous parameters of seismic signal[J]. Oil Geophysical Prospecting (in Chinese), 42 (4) : 424–430. |
| [16] | Stockwell R G.2007. A basis for efficient representation of the S-transform[J]. Digital Signal Processing, 17 (1) : 371–393. DOI:10.1016/j.dsp.2006.04.006 |
| [17] | Taner M T, Koehler F, Sheriff R E.1979. Complex seismic trace analysis[J]. Geophysics, 44 (6) : 1041–1063. DOI:10.1190/1.1440994 |
| [18] | Torres M E, Colominas M A, Schlotthauer G, et al. 2011. A complete ensemble empirical mode decomposition with adaptive noise[C].// IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Prague: IEEE, 4144 -4147. |
| [19] | Vakman D.1996. On the analytic signal, the Teager-Kaiser energy algorithm, and other methods for defining amplitude and frequency[J]. IEEE Transactions on Signal Processing, 44 (4) : 791–797. DOI:10.1109/78.492532 |
| [20] | Wu Z H, Huang N E.2009. Ensemble empirical mode decomposition: A noise-assisted data analysis method[J]. Advances in Adaptive Data Analysis, 1 (1) : 1–41. DOI:10.1142/S1793536909000047 |
| [21] | Xi Xiaoyu, Liu Hong.2007. The application of HHT method in study of seismic cycle[J]. Journal of Jilin University (Earth Science Edition) (in Chinese), 37 (3) : 624–628. |
| [22] | Yang Peijie, Yin Xingyao, Zhang Guangzhi.2007. Seismic signal time-frequency analysis and attributes extraction based on HHT[J]. Progress in Geophysics (in Chinese), 22 (5) : 1585–1590. DOI:10.3969/j.issn.1004-2903.2007.05.037 |
| [23] | 陈学华, 贺振华, 黄德济, 等.2009. 时频域油气储层低频阴影检测[J]. 地球物理学报, 52 (1) : 215–221. |
| [24] | 段生全, 贺振华, 黄德济.2005. HHT方法及其在地震信号处理中的应用[J]. 成都理工大学学报(自然科学版), 32 (4) : 396–400. |
| [25] | 高静怀, 汪文秉, 朱光明.1997. 小波变换与信号瞬时特征分析[J]. 地球物理学报, 40 (6) : 821–832. |
| [26] | 黄光南. 2009. 希尔伯特黄变换及其在地震资料分析处理中的应用[D]. 青岛: 中国海洋大学. |
| [27] | 刘庆敏, 杨午阳, 李书平, 等.2009. 希尔伯特-黄变换在地震资料处理中的应用[J]. 西北地震学报, 31 (3) : 211–216. |
| [28] | 皮红梅, 刘财, 王典.2007. 利用Hilbert-Huang变换提取地震信号瞬时参数[J]. 石油地球物理勘探, 42 (4) : 424–430. |
| [29] | 葸晓宇, 刘洪.2007. HHT方法在研究地震旋回体中的应用[J]. 吉林大学学报(地球科学版), 37 (3) : 624–628. |
| [30] | 杨培杰, 印兴耀, 张广智.2007. 希尔伯特-黄变换地震信号时频分析与属性提取[J]. 地球物理学进展, 22 (5) : 1585–1590. DOI:10.3969/j.issn.1004-2903.2007.05.037 |