基于CEEMDAN的地震信号高分辨率时频分析方法

尚平萍,李 鹏,杨安琪,陈学国

(1.中国石油化工股份有限公司石油物探技术研究院,江苏南京211103;2.中国石油化工股份有限公司胜利油田分公司勘探开发研究院,山东东营 257015)

近年来,随着油气勘探的难度不断增大,时频分辨率低的传统谱分解技术,越来越难以满足高分辨率地震数据处理和解释的需求,高分辨率地震时频分析技术的研究与应用成为一种趋势。目前常用的地震时频分析技术主要包括短时傅里叶变换、小波变换、S变换和魏格纳-威利分布(Wigner-Ville distribution,WVD)变换等[1],但这些技术均存在一定不足:短时傅里叶变换的时窗长度固定,确定时窗长度后无法调节时间和频率分辨率[2-3];连续小波变换克服了短时傅里叶变换时窗固定的限制,具有多分辨率特性,信号特征刻画效果好,但受Heisenberg不确定性原理影响,时频分辨率有限[4];S变换是小波变换和短时傅里叶变换的组合,它的窗函数形状固定,且仍受Heisenberg测不准原理的限制[5-7];魏格纳-威利分布变换虽然能够在时频域很好地刻画单一频率成分的信号,但复杂信号会产生严重的交叉干扰,影响WVD变换结果的解释,尽管相关的平滑改进方法能较好地抑制交叉项,但是以降低信号的局部时频精度为代价的[8]。可见,上述方法虽然具有简单快速的优点,但是都存在一定的局限性。

1998年,HUANG等[9]提出了一种适用于非线性非平稳信号的时频分析新方法—希尔伯特-黄变换(HHT),在信号处理研究中受到广泛关注,该方法首先对非线性非平稳信号进行EMD,获得固有模态函数(intrinsic mode function,IMF)分量后,再进行希尔伯特变换(HT),从而得到分辨率较高的Hilbert时频谱。然而,当信号中存在间断信号或噪声干扰时,EMD会产生难以克服的模态混叠现象,破坏了每个IMF分量所蕴含的物理意义,降低了分解的准确性,使得基于EMD的地震数据处理方法难以推广。

为克服EMD的模态混叠效应,WU等[10]提出了集合经验模态分解(ensemble empirical mode decomposition,EEMD)方法,使信号在不同尺度上具有连续性,但是该方法是以增加计算成本为代价来提高集合平均次数以降低重构误差[11]的,所以它不是一个完备的分解方法。针对EEMD方法所存在的问题,TORRES等[12]提出了CEEMDAN方法,该方法可获得更好的模态频谱分离结果,而且在较低的计算成本之下可重构出精确的原始信号。

本文以CEEMDAN代替EMD并将其与希尔伯特变换(HT)相结合,形成了一种适应于非线性非平稳信号的高分辨率时频分析方法。该方法对CEEMDAN得到的每个IMF分量进行HT谱分析,计算其瞬时振幅和瞬时频率,然后生成瞬时频谱,用于地震时频分析。对合成记录和实际数据,采用不同的时频分析方法得到的时频谱进行对比,以证明基于CEEMDAN的时频分析方法的准确性和稳定性。该方法兼具较高的时间分辨率和频率分辨率,可获得更精确的时频谱,因此作为一种地震资料精细解释的有力工具,为后续的含油气检测提供了有力支撑。

1 方法原理

希尔伯特-黄变换(HHT)主要包括Huang变换和Hilbert变换两个部分[13],Huang变换又称作EMD[14-16]。我们先利用EMD获得数目有限的IMF分量,再利用瞬时频率方法和HT变换获得信号的时频谱。EMD将信号分解为许多单分量信号的和[17],在该过程中,IMF分量频率与分离的先后顺序有关,分离时间越早IMF分量频率越高,反之分离时间越迟则IMF分量频率越低[18],原信号的最高频率成分存在于最早得到的IMF分量中。由于经过了多通带滤波处理,每个IMF分量都是稳态的。对各IMF分量进行HT变换,可获得三瞬谱:瞬时相位谱、瞬时振幅谱和瞬时频率谱。

1.1 HHT存在的模态混叠问题及改进方法

尽管与传统时频分析方法相比HHT在求取地震记录的瞬时参数时,可获得更准确的结果,但当信号不连续或者存在噪声时,EMD会导致模态混叠现象,这基本是不可避免的,此时每个IMF分量所蕴含的物理意义都受到影响,导致求取的瞬时参数不理想,解释结果存在误差,制约了其实际应用范围。

为了克服或者降低模态混叠效应,一些学者提出了不同的解决方案。FLANDRIN等[19]设定了截止频率以分离间歇信号,但是门槛值的选择具有很强的主观性,故难以获得较好的效果。EEMD是对EMD的一种改进,其核心思想在于分解原始信号之前加入频谱分布均匀的高斯白噪声[11]。EEMD的基本思想是先在原始信号中添加高斯白噪声,再对其整体进行EMD,由于高斯白噪声的频谱分布均匀,这使得不同尺度的信号成分自动映射到相应的参考尺度上,对分解出的分量进行平均以此抵消前期加入的白噪声。EEMD算法中有两个重要的参数需要预先设定:添加高斯白噪声的幅值和集合平均次数,HUANG等[9]依靠经验确定这两个参数,既不方便也不客观。作为一种噪声辅助分析方法,EEMD在一定程度上克服了EMD的模态混叠效应,改善了分解效果,但同时引入了其它问题,如残余噪声、分量解非唯一、白噪声幅值和集合平均次数不能自适应地获取等。针对上述问题,YEH等[20]提出了一种改进的EEMD方法,即互补集合经验模态分解(complete ensemble empirical mode decomposition,CEEMD)方法,可消除IMF分量的残余噪声。

在CEEMD方法中,原始信号中加入成对的白噪声,可产生两个IMF分量集合。加入噪声后的信号如下:

(1)

式中:S为原信号;N为白噪声;M1、M2分别为加入正、负成对噪声后的信号。

1.2 自适应噪声的总体集合经验模态分解方法原理

自适应噪声的总体集合经验模态分解(CEEMDAN)[12]本质上也属于EMD,是对其改进后得到的一种变换形式[21]。与EEMD方法类似,该方法也属于噪声辅助分析方法。但不同之处在于:CEEMDAN方法克服了模态混叠效应,可精确地重构出原始信号,并获得了较好的模态分离谱,同时提高了计算效率。

首先定义算子Ej(·)为给定信号通过EMD得到的第j阶模态分量;ωi(n)表示高斯白噪声,其中i=1,2,…,I;εk系数表示在每个阶段的信噪比,其中k=0,…,K。设x(t)为目标信号,CEEMDAN方法包括以下6个步骤[22]。

1) 第一阶固有模态分量IMF1(t)的求取方法与EMD方法相同,基于不同的噪声,重复EMD分解过程I次,计算集合平均值,目标信号x(t)的IMF1(t):

(2)

2) 当k=1时,计算一阶残差r1(t):

(3)

3) 对r1(t)添加经EMD分解后的噪声分量E1[ωi(t)],再进行一次EMD分解,并定义集合平均值为第二阶固有模态分量IMF2(t):

(4)

4) 当k=2,…,K时(其中,K是模态总体数量),计算k阶残差rk(t):

(5)

5) 提取rk(t)+εkEk[ωi(t)]的IMF分量,计算它们的集合平均值,得到目标信号的第k+1阶固有模态分量IMF(k+1)(t):

(6)

6) 重复上述步骤,当残差不能再被分解时,得到的结果即为最终残差R(t):

(7)

目标信号x(t):

(8)

上式表明,CEEMDAN方法可以精确地重构原始信号。

EEMD方法中所添加的高斯白噪声的幅值和集合平均次数不能自适应地获取,针对这一问题,CEEMDAN方法首先设置期望的信号分解相对误差e,考虑到EMD方法得到的第一个IMF分量不一定能够准确描述信号中的高频信息,故可以通过抑制低频信号来获取高频信号,该过程是通过添加白噪声实现的。根据相关准则[12]确定白噪声的幅值范围:

(9)

式中:ε=σn/σs为加入白噪声幅值标准差σn与原信号幅值标准差σs的比值;ξ=σh/σs为信号中高频成分幅值标准差σh与原始信号幅值标准差σs的比值。因此,公式(9)可表示为:

(10)

由公式(10)可知加入的白噪声的幅值范围。设定e值之后,再根据HUANG等[9]提出的统计规律可得到集合平均次数N:

(11)

2 合成信号时频谱分析

分别采用EMD、EEMD及CEEMDAN 3种方法对合成信号进行测试,以验证CEEMDAN方法的优势。图1a是各组分信号及其合成信号,自上而下分别是10Hz(0～1s)和20Hz(1～2s)的余弦信号、100Hz的Morlet子波(1.8s处)、50Hz的间歇性余弦信号以及合成信号。图1b是采用EMD方法分解合成信号得到的结果,可以发现EMD方法将合成信号分解成6个IMF分量,分量存在明显的模态混叠效应。

图1c是添加10%的高斯白噪声且集合平均次数为500时利用EEMD方法分解合成信号得到的结果(图中只列出前6个IMF分量),可以看出,EEMD方法已在很大程度上降低了模态混叠效应,IMF1和IMF2两个分量中分别恢复出了100Hz的高频Morlet子波和50Hz的余弦信号;然而IMF2中掺杂了20Hz的余弦信号,并且IMF3也混合了20Hz和10Hz的余弦信号,这说明仍存在模态混叠效应,但是相比EMD方法得到的结果,其模态混叠效应已经有了明显降低。

CEEMDAN方法采用了与EEMD方法相同的高斯白噪声,但集合平均次数为100,分解结果如图1d所示。图中只列出了前6个IMF分量,可以看出,IMF1完全恢复了100Hz的高频子波,基于IMF2、IMF3和IMF4可完全提取50Hz、20Hz和10Hz的简谐波分量。采用CEEMDAN方法得到的结果几乎完全克服了模态混叠效应的影响。相较于EEMD方法,CEEMDAN方法得到的分解结果更准确。

从图2中可以看出,采用CEEMDAN方法重构得到的信号误差小,几乎可以忽略。可见,CEEMDAN方法可以精确地重构原始信号。

图1 各组分信号及合成信号和应用EMD、EEMD、CEEMDAN方法分解结果a 各组分信号及合成信号; b EMD方法分解结果; c EEMD方法分解结果(前6个分量); d CEEMDAN方法分解结果(前6个分量)

图2 采用CEEMDAN方法得到的重构信号

为了验证CEEMDAN方法的时频高分辨率特性,基于图1a中的合成信号进行了试算,采用不同时频分析方法获得的时频谱如图3所示。图3a、图3b和图3c分别为采用STFT(short-time Fourier trans-form)、CWT(continuous wavelet transform)和广义S变换方法得到的时频谱,可以看出,采用STFT方法得到的时频谱分辨率单一,能量聚焦性低;采用CWT和广义S变换方法得到的时频谱在低频端频率分辨率高,在高频端时间分辨率高,但因受Heisenberg测不准原理限制,时间和频率分辨率不能同时达到最佳。图3d、图3e和图3f分别为采用EMD、EEMD和CEEMDAN方法得到的时频谱,这3种方法都克服了传统方法的不足,得到了较为精确的结果,但是采用EMD方法得到的时频谱出现了严重的模态混叠现象;采用EEMD方法得到的时频谱模态混叠现象虽有改善,但仍存在;采用CEEMDAN方法得到的时频谱几乎不存在模态混叠现象,各频率分量都清晰可辨。可见CEEMDAN方法得到的时频谱最为精细和准确。

图3 对合成信号采用不同时频分析方法得到的时频谱a SFFT; b CWT; c 广义S变换; d EMD; e EEMD; f CEEMDAN

3 实际应用

将CEEMDAN方法应用于某工区的实际地震资料处理,以检验该方法的实际应用效果,尤其是该方法的时频高分辨率特性。图4对比了分别采用STFT、CWT、广义S变换及CEEMDAN方法得到的时频谱结果,可以看出相较于其它3种时频分析方法,采用CEEMDAN方法得到的时频谱时间分辨率和频率分辨率明显提升,时频谱主频的变化特征清晰,分辨率高,这有助于描述地质细节,也证明了该方法的有效性。从图5可以看出,该工区应用CEEMDAN方法后得到的地震剖面的分辨率明显提高,地质细节刻画清晰。

图4 某工区实际单道地震记录及采用不同时频分析方法得到的时频谱a 实际单道地震记录; b STFT; c CWT; d 广义S变换; e CEEMDAN

图5 某工区应用CEEMDAN方法前(a)、后(b)的地震剖面

研究表明,地震波的吸收衰减性质主要取决于岩石骨架、孔隙度和孔隙流体性质,当地层含油气时,地层对高频成分吸收增强,导致高频成分快速衰减,利用这一特性可以预测地层的含油气性。如图6的钻测井数据所示,红色表示含油层,橙色代表连井的砂体分布,1井2砂组油层厚度8.6m,1井1砂组含油水层厚度15.2m;2井1砂组油层厚度14.5m，有效厚度6.7m。图7中黑色表示频谱变化斜率属性高值,红色表示频谱变化斜率属性较高值,黄色表示频谱变化斜率属性较低值,绿色表示频谱变化斜率属性低值,对比图7a和图7b可以看出,应用CEEMDAN方法后砂组预测结果与测井数据更匹配,高频信息更丰富。将给定频率低频分量的频谱变化斜率减去应用CEEMDAN方法后的频谱变化斜率,得到的吸收属性剖面如图8所示,利用该属性可较好地预测储层的含油气性。图8中黑色表示吸收属性高值,红色表示该属性的较高值,黄色表示该属性的较低值,绿色表示该属性的低值,可以看出1井1砂组含油气性较差,2井1砂组含油气性较好,这与测井数据吻合,应用CEEMDAN方法后得到的剖面高频信息更丰富,具有较高的分辨率,有利于提升薄油层的预测效果。

图6 某工区1井、2井储层对比

图7 某工区原始地震剖面频谱变化斜率属性剖面(a)及应用CEEMDAN方法后频谱变化斜率属性剖面(b)

图8 某工区连井线吸收属性剖面

4 结束语

EMD作为HHT的核心技术,存在严重的模态混叠效应,信号分析的精度不高。基于EMD的改进算法CEEMDAN,不仅克服了模态混叠效应,而且能够精确地重构原始信号,同时能够自适应地确定EEMD算法中需预先确定的两个重要参数:添加高斯白噪声的幅值和集合平均次数。

本文以CEEMDAN代替EMD并与HT相结合,提出了一种适用于非线性非平稳信号的高分辨率时频分析方法。信号频谱分析和实际应用结果表明,与STFT、CWT、广义S变换方法等传统的时频分析方法相比,本文方法能获得较高的时间和频率分辨率,得到更精准的信号频谱。因此基于CEEMDAN的时频分析方法可作为地震资料精细解释的有力工具,应用于精细沉积微相的分析以及薄层储集体的刻画。