经验模态分解法在地震资料处理中的应用

李世涛 张勇 范存辉 梁亮

(1.中海油能源发展股份有限公司上海采油技术服务分公司，上海 200030;2.西南石油大学油气藏地质及开发工程国家重点实验室，成都 610500)

小波变换在石油勘探中已被广泛应用，它本质上是一种自适应加窗富氏变换，要求小波窗内的信号短时平稳，没有摆脱富氏变换的局限性;其次，小波基的多样性及零相位不适应非零相位子波生成的信号，即小波基的有限长度还会造成信号能量的泄漏:因此，基于小波变换的时频分析很难适应非平稳时间序列的时频分析。为此，本次研究应用Huang等人(1996)提出的经验模态分解法(Empirical Mode Decomposition，简称 EMD)中的固有模态函数(Intrinsic Mode Function，简称IMF)作为基函数，对非平稳时间序列的时频分析方法进行深入研究，把该研究成果应用于地震和测井资料联合处理解释之中。

1 经验模态分解法时频分析

经验模态分解是对一种非线性、非平稳信号的新的时频分析方法。它能够对数据进行线性化、平稳化处理，并在处理过程中保留信号自身特性;它利用非平稳信号的统计特征尺度随时间信号变化的局部时变特征，从原信号中提取若干个固有模态函数(IMF)和一个残余量，用多个固有模态函数表征信号的局部特征，用残余分量表征信号缓慢变化的趋势。固有模态函数需满足2个条件:(1)整个数据集合中，极值点数目和过0点数目必须相等或最多相差一个;(2)在任何一点局部区间，由局部最大值和最小值形成的包络均值为0。下面介绍非平稳信号经验模态分解算法。

(1)EMD分解步骤，分解过程采用以下“筛分”算法。

第1步，对信号x(t)求取所有的局部极值点，用样条函数分别拟合所有的局部极大值点和极小值点，形成上下包络线bm(t)和bn(t)。

第2步，计算上下包络线的均值:

并将原始序列x(t)减去包络均值:

此时判别h1，1(t)是否满足IMF的2个条件。如果不满足，将h1，1(t)作为原始数据重复上述处理过程，直至序列

满足 IMF的条件。这样，就得到 IMF的第一个分量:

第3步，由信号x(t)分离出第一个分量c1(t)，得剩余量:

将r1(t)作为一个新的原始序列，重复上述步骤，依次求出第2，3，…，n个 IMF ci(t)，当剩余量rn(t)成为一个单调函数或小于预定值时分解结束。

这样，原始信号x(t)可表示为:

(2)筛分终止条件及样条函数的选择。在实际IMF求取时，IMF的第2个条件很难满足。为此，应用Huang等人提出的收敛准则，定义:

当Is≤ξ时，则终止筛分。对局部求取的上下极值点用样条插值函数插值组成包络。插值函数的边界端点处理亦是EMD分解的一个难题。本次研究运用自适应外推法进行边界处理。

2 希尔伯特变换和瞬时谱计算

由式(6)可知，信号x(t)经EMD分解得到的每一个IMF的ci(t)均为单分量。现对每一个单分量进行富氏变换:

式中:P表示柯西积分主值。由此可得解析信号:

Am(t)、θ(t)为瞬时振幅和相位，其计算公式为:

由瞬时相位可求得瞬时频率:

将振幅值用时间和频率表示，可得到希尔伯特幅值谱:

将Hm(ω，t)对时间积分，得边间谱:

由希尔伯特幅值谱和边间谱，可给出信号序列的平稳度函数:

对应时间域有:

式中:

3 应用实例

3.1 地震资料的EMD分解

图1展示了地震道xi(t)，EMD分解全过程包括IMF单分量ck(t)及分解时产生的逐级上、下极值包络bkn(t)。

图1 地震道EMD分解全过程

图2 图1中EMD分解后的IMF k

图2 展示了图1中EMD分解后的IMFk包络bkn和单分量ck结果 。图3显示的是地震道x(t)，通过EMD分解得到的IMF分量ck(t)及IMF的第2个条件筛分准则中的近似包络值bkn(t)。IMF分量ck(t)的频率随着k的增大逐渐降低，波长逐渐增大。各个分量ck(t)包含有不同的时间尺度及分辨率的信号特征，这种特征与信号自身的特点是自适应的。其EMD分解后有6个IMF单分量。以下逐级IMF分量ck(t)的值及包络平均值bkn(t)逐次依0.5的几何级数衰减，为方便展示，全部作归一化处理。图3表明，具有高频成分的主分量c1(t)的平稳性较差，具有低频成分的平稳性较好。

图3 图1中EMD分解后局部放大的IMF k

3.2 地震资料EMD分解IMF主分量分析剖面效果

图4 是某测线偏移剖面图，图5展示了该地震剖面EMD分解后得到的主分量、次分量及第3分量剖面。图6是某测线井段附近剖面的EMD分解主分量剖面与小波分解剖面图，展示了主分量剖面与分解的1级包络剖面的关系:在包络剖面上绝对幅度越大表示该处的平稳性越差、地层的非白噪性越强;在主分量剖面上，几处低频高幅处是储层及含油气层，而小波分解剖面上没有这种特征。

图4 某测线偏移剖面

图5 某测线EMD单分量IMF剖面

图6 井段剖面EMD分解与小波分解剖面

4 地震资料EMD分解后单分量IMF k的时频分析

地震资料x(t)在经过EMD分解后获得近似平稳的固有模态函数IMFk，个数比较少只有6～7个，很难直接应用这么少量的IMF作时频分析。因此采用基于小波基的时频分析方法对逐个IMFk作时频分析的有益尝试，取得了一定成效。

图7是地震道EMD及连续小波变换CWT分解后所作的瞬时Hilbert幅值谱对比图。两者对比表明EMD幅值谱较CWT幅值谱能量强的特点。图8是地震道EMD分解后计算得到的时域和频域的平稳度曲线。

由图7—图9可知，EMD分解后所作的瞬时Hilbert幅值谱比连续小波变换CWT分解后的谱的分辨率高，并与地震道波形匹配自然。

图7 瞬时Hilbert幅值谱对比

图8 地震道EMD分解的平稳度曲线

图9 地震剖面对应的EMD分解Hilbert 谱

5 结语

本次研究对经验模态分解法在地震数据处理领域中的应用作了探索论证，并取得一定的应用效果。运用经验模态分解法和小波变换2种方法进行试验比较，得到4个方面的结论。

(1)在信号分频处理中，基于小波变换的时频分析应用方便灵活，具有广泛的应用性，优于常规富氏变换的分频处理技术;信号的小波分解是一个褶积过程，但不利于非平稳数据的时频分析。

(2)在地震数据x(t)的EMD分解中，由于分解得到的IMFk一般很少，而且只有前几个分量在地震数据时频分析的频段内，因此一般较难用于时频分析。而EMD的第一个分量c1(t)具有较宽的频率成分，因此，与地震解释配合，EMD方法的主分量可用于地震属性分析。

(3)由于测井资料与地震数据的相似特性，可在重构测井曲线，储层参数反演方面有良好的应用前景。

(4)在对地震资料的EMD分解作时频分析时，插值函数的边界端点处理是EMD分解的一个难题。本次研究采用自适应外推法作为边界条件，但需要进一步研究模态函数IMFk之间的分形插值方法，使少量的IMFk插成较密的IMFk簇，提高地震资料EMD分解法的时频分析效果。

［1］ Huang E N，Wu C M，Long R S，et al.A Confidence Limit for the Empirical Mode Decomposition and Hilbert Spectral Analysis［C］//Proceedings Royal Society of London.2003:2317-2345.

［2］杨世锡，胡劲松，吴昭同，等.基于高次样条插值的经验模态分解方法研究［J］.浙江大学学报(工学版)，2004，3(3):267-270.

［3］程磊，瞿伟廉.基于Hilbert-Huang变换理论的非平稳数据处理［J］.建筑科学与工程学报，2007，24(1):26-30.

［4］王祝文，刘菁华，聂春燕，等.Hilbert-Huang变换在提取阵列声波信号动力特性中的应用［J］.地球物理学进展，2008，23(2):450-455.

［5］ 俞寿朋.宽带 Ricker子波［J］.石油地球物理勘探，1996，31(5):605-615.

［6］田芳，高静怀，陈文超，等.广义S变换与薄互层地震响应［J］.地球物理学报，2003，46(4):526-532.

［7］盖广洪.经验模态分解的一种改进算法［J］.西安交通大学学报，2004，38(11):1199-1202.