改进的EMD时频分析地震资料高分辨率处理方法及其应用研究

吴仕虎, 肖富森, 卢柃岐, 冉 崎, 张福宏, 杨广广

(1. 中国石油西南油气田分公司 勘探开发研究院，成都 610041；2. 成都理工大学，成都 610059；3. 成都理工大学 数学地质四川省重点实验室，成都 610059)

0 引言

随着全球油气勘探的不断发展，勘探难度也在不断增大，获得高保真、高分辨率的地震资料对复杂油气藏的勘探开发具有重要意义。地震信号是由激发的地震子波经过大地滤波后形成，由于地层介质滞弹性和非均质性，地震波在地层传播过程中受大地滤波、吸收等影响，使地表接收到的地震记录高频成分损失严重[1-2]，导致地震资料的分辨率、信噪比降低。为了使地震资料能满足储层研究和油藏描述的要求，很多地球物理学家提出并发展了地震资料高分辨率处理的方法，常见的有谱白化、反Q滤波、时频分析、多尺度联合分析等几类。

谱白化是通过展宽振幅谱来提高信号分辨率的方法，它不改变子波的相位谱，可以在时域中进行，也可以在频域中进行。有多位学者进一步提出了小波谱白化和Hilbert谱白化方法，小波变换克服了傅里叶变换不能分析信号局部特征的缺陷[3-6]；Hilbert谱白化可在时、频域内同时增强信号的局部细节信息，使地震剖面更为清晰和连续，且具有更高的信噪比。反Q滤波是一种补偿大地吸收衰减效应的技术，它不仅可以补偿振幅衰减和频率损失，而且还可以改善记录的相位特性,从而改善同相轴的连续性，提高弱反射波的能量和地震资料的信噪比、分辨率。早期的工作是依据Futterman数学模型提出的用级数展开作近似高频补偿的反Q滤波，进而提出振幅、相位同时进行稳定补偿的反Q滤波，并将其发展到Q值随时间或深度连续变化的情况，近年来发展出利用Futterman衰减模型同时对振幅衰减和速度频散进行校正，具有较高的稳定性和计算精度[7-11]。时频分析在时频域上分析非平稳信号，它克服了傅里叶分析时域与频域分离的缺陷，较准确地反映时间与频率信息。目前发展了基于短时傅里叶变换的地层吸收补偿技术、基于小波包分解的地层吸收补偿方法、可调节窗函数标准差的非对称广义S变换，对地震记录进行补偿和基于改进的广义S变换进行衰减补偿[12-17]。多尺度联合分析方法是利用测井、井间地震、VSP等地球物理技术，对地下同一目标地质体进行不同尺度性质的反映，通过它们之间的联合作用提高地震资料分辨率[18-21]。

上述方法对提高地震资料分辨率均有不错的效果，然而这些方法或适用于平稳信号、或需要假设条件、或需要提前获取Q值，局限性大，工作量大，成本高。为了以更经济的方式得到真实可靠的高分辨率的地震资料，笔者提出了一种适用性广、成本经济、合理有效地提高地震资料分辨率的技术方法，具有借鉴作用。本文方法综合考虑信号在时域与频域的特征，在时频域上对信号进行高分辨率处理，首先利用EMD技术把地震信号自适应分为有限个本征模态函数IMF；对每一个IMF分别进行高分辨率处理，然后叠加得到高分辨率的地震数据。理论资料与实际资料的处理结果表明，此方法具有较好的应用效果，为解决高频补偿问题提供了一条全新的思路。

1 方法原理

1.1 EMD的方法原理

经验模态分解(Empirical Mode Decomposition, EMD)是从信号本身出发，通过层层筛选，先将时间特征尺度小的高频IMF(Intrinsic Mode Function, IMF)分量分离出来，然后将时间特征尺度大的低频IMF分量分离出来，根据停止准则，最后得到一个近似单调的残余分量。为了去除叠加波的影响和使波形更加对称，EMD的分解过程采用“筛分”算法。对于一个实信号x(t)，按照流程图1的步骤进行分解，原信号x(t)被分解成n个IMFci(t)和一个信号余量rn(t)之和(式(1))。

(1)

图1 EMD的分解步骤Fig.1 Steps of EMD decomposition

其中每个IMF必须满足以下两个条件：

1)在整个数据集中，该信号极值点的数目与零交叉点的数目相等或者最多相差一个点。

2)在任意时间上，该信号由局部最大值点形成的上包络线和局部最小值形成的下包络线的均值为零(即上包络、下包络局部是对称的)。

IMF的第一个条件是为了满足高斯正态平稳过程的传统窄带要求；第二个条件把传统的全局性限定为局域性，目的是避免模态混叠，保证一个固有模态函数仅仅包含一个基本模式的振荡，去除因波形不对称而造成的波动。这样获得的IMF信号在任何时刻都只具有单一的频率成分[22]。

基于EMD分解的时频分析是一种基于经验的分解，而不是基于数学表达式的分解，它直接从信号本身产生基函数，没有交叉干扰项，因此具有很好的自适应性、正交性和完备性。除此之外，EMD时频分布还具有很好的时频聚集性和高分辨率，时频分布平面内的能量分布十分密集，可以快捷地给出有效的时频表示，突出时变信息。

1.2 地震记录的振幅补偿

因地层的吸收衰减作用，地震记录的振幅谱有相对高频的变化。设实际的地震记录为S(t)，变换到频域为：

S(f)=A(f)eiφ(f)

(2)

式中：A(f)为地震记录的振幅谱;φ(f)为相位谱。采用平滑滤波的方法，可以由A(f)估算出地震记录振幅谱的衰减曲线(即通常的地震子波的振幅谱)。该滤波器的平滑因子采用三角滤波算子(式(3))。

(3)

n=1,2,L…,2N-1(N为半滤波算子长)

式中:n为滤波算子样点序号。用H(n)对地震记录的振幅谱作平滑滤波，得到振幅谱的衰减曲线B(f)，如式(4)所示。

(4)

振幅补偿曲线C(f)为衰减曲线的倒数，即C(f)=1/B(f)。但是在实际的计算中需对衰减曲线加入适当的白噪，以保证振幅谱补偿曲线的稳定性。为了保证振幅补偿的不超过地震记录的有效频带范围，常对振幅补偿曲线作限频处理[23-24](式5)。

(5)

式中：fLC、fLP、fHP、fHC、fNq分别为地震信号的低截、低通、高通、高截、Nyquist频率；T是与白噪系数U有关的稳定系数，如式(6)所示。

(6)

为了从地震记录恢复反射系数振幅谱，还应对振幅谱补偿结果进行线性加权处理，加权因子为式(7)。

f=fLC,fLC+Vf,L,fHC

(7)

式中：K为高频增强系数，一般取大于1的数；Vf为频率采样间隔。

为了使本方法能保持相对振幅，需用输入记录的振幅谱予以标定，以保证处理前、后的振幅谱的总能量不变。振幅补偿的最终结果为式(8)中Y(f)。

Y(f)=A(f)·C(f)·

(8)

1.3 计算流程

笔者所述的地震信号高分辨率处理是以EMD方法为基础，针对地震信号振幅特性而展开的地震信号高分辨率处理，其中所述的地震信号为用于勘探地下岩层性质和形态的时域信号。具体的步骤：①获取地震剖面，对一道地震信号采用EMD方法自适应分解，得到固有模态函数IMF分量；②分析该地震信号、EMD分解得到的各个IMF分量的振幅谱，分别确定高分辨率处理的高截、高通、低通、低截阈值；③利用选择的阈值分别对各个IMF分量进行高分辨率处理；④将所有时间范围内高分辨率处理后的IMF分量叠加重构；⑤为了减少分解、重构等过程中引起的误差，可根据实际情况对重构的信号进行平滑等修饰性处理；⑥重复步骤①～步骤⑤完成所有单道的处理，最终获得高分辨率的二维地震剖面，用于地震解释。同时需注意，不同区块的地震资料差异较大，在选择参数时应根据实际的情况进行试验后灵活处理，才可以最终达到最佳的处理效果。为了清晰地展示本算法，其流程图如图2所示。

图2 基于EMD方法的高分辨率处理算法流程图Fig.2 Flow chart of high-resolution processing algorithm based on EMD method

2 理论模型应用效果分析

2.1 理论模型单道的效果分析

为了验证该方法的可行性，设计了单道理论模型：用主频为20 Hz的雷克子波(图3(a))和6组数值为1.2的反射系数(图3(b))进行褶积合成信号(图3(c))，并展示出该合成信号的振幅谱(图3(d))。在设计的过程中每组反射系数的时间距离遵循时间间距从0.121 s开始，以0.002为步长增加(表1)。

表1 每组反射系数之间的距离

图3 理论模型的设计Fig.3 Design of theory model(a)20 Hz的雷克子波；(b)一系列反射系数；(c)合成信号；(d)合成信号的振幅谱

图4 合成信号的IMF信号与IMF振幅谱Fig.4 IMF signal and IMF amplitude spectrum of synthesized signal(a)IMF的信号；(b)IMF的振幅谱

通过图3(c)，可以清晰地看出,该合成信号第1个波形完全没有分开，第2个波形几乎没有分开，从第3个波形开始，反射系数分开程度逐渐增大。对该合成信号用EMD方法分解得到5个IMF信号(图4(a))，其对应振幅谱见(图4(b))。利用本文方法对上述合成信号进行提高分辨率处理，处理效果见图5。由图5(a)可以看出，本文EMD时频分析高分辨率处理方法能够有效地将一部分发生混叠的信号分开，从而使信号的分辨率得到提高。因该合成信号的采样点数较少，笔者对处理前、后的振幅谱求包络，得到合成记录处理前后的频谱包络对比图(图5(b))。通过观察高频部分的振幅值有所提高，对应说明能量得到一定程度补偿。通过对比图5可知，经过本文方法的处理，该合成信号的主频得到提高，高频信号能量明显得到提高，整个频带宽带明显拓宽。

2.2 理论剖面模型的效果分析

理论剖面实验模型，是由35个相同单道排列形成的剖面(图6)。单道是由6个间距呈0.002 s递增的反射系数(表1)与主频为20 Hz的雷克子波褶积而成，采样间隔为1 ms，采样点数为1 000。

图5 高分辨率处理前、后的效果对比Fig.5 Comparison of the effects before and after high resolution processing(a)处理前后的信号对比图；(b)处理前后的振幅谱包络对比图

图6 高分辨率处理前、后的地震剖面Fig.6 High-resolution seismic section before and after processing(a)原始地震剖面；(b)高分辨率处理后的地震剖面

从图6(a)可以看出，在0 s～0.4 s之间，两个波形是出现复合的现象，波形叠加在一起无法分辨，而且从下到上，波形的复合现象是越来越严重，两个波形之间越来越难分辨。

理论剖面经过本文方法处理之后，如图6(b)所示，0.2 s～0.4 s之间的波形能够辨认，0.4 s～0.8 s之间的波形较之前分开明显，分辨率明显提高，0.8 s～1 s之间的两个波形完全分开。

3 实际资料应用效果分析

3.1 实际资料剖面的分析

为了进一步验证该方法适用性，将该方法应用于四川盆地川中地区高石梯井区的二维地震资料，该资料的目的层是灯影组，总共4 316道，每道采样点为2 450个，采样时间间隔为1 ms，实钻井位于第2 576道。按照本文的方法原理进行高分辨率处理，首先采用EMD方法对单道地震数据进行分解，然后对分解之后得到的IMF分量进行高分辨率处理，重构得到高分辨率的单道地震数据，最后逐道依次进行上述方法处理流程，重构得到高分辨率的二维地震资料。因该资料的道数、采点数较多，现从原二维地震资料中截取部分形成一个新的地震剖面进行效果的比较分析，截取第2 300 道～3 000道，共701道，采样时间2 s～2.45 s，共451个采样点，在形成的新剖面中钻井位于第277道。为了使呈现的图像更加清晰、美观，显示了第240 道～300道，采样时间0 s～0.4 s的剖面波形图，处理前、后的对比如图7所示。

图7 二维实际资料处理前后对比图Fig.7 Comparison of 2D actual seismic section before and after processin(a)实际资料的原始剖面；(b)实际资料的高分辨率处理后剖面

通过图7(a)与图7(b)对比可知，在0.03 s～0.05 s、0.07 s～0.09 s两个时间区间分辨出薄互层反射信号(如图7(a)和图7(b)中红色箭头所示)。高石梯地区该套地层岩性为砂泥岩薄层互层，经高分辨率处理后的地震剖面不仅保持了原始剖面的基本面貌，而且分辨率明显提高、剖面的成层性更好，一些弱化的层位信息得到加强，缺失的相关层位信息得到了显现。同时，可以清晰地看到，处理后的剖面上实钻井的灯影组顶部界面更清晰，更容易区分，灯影组下部的第1 波峰～2 波峰与实钻井的人工合成地震记录匹配度较好，反射波同相轴的连续性得到明显的改善，分辨率提高(如图7(a)和图7(b)中红色的矩形框所示)。这说明本文方法能够压缩子波的长度，提高剖面的分辨率，有利于进一步的地质体识别和解释。

3.2 实际资料剖面的平均振幅谱分析

图8为实际资料剖面处理前后的振幅谱对比，其中图8(a)是处理前地震资料剖面的平均振幅谱，图8(b)是处理后地震资料剖面的平均振幅谱。通过对比图8(a)和图8(b)可知，以主频为中心的优势频带提高约12 Hz，处理前地震资料的主频约为29 Hz，而处理后的资料主频提高到41 Hz左右，频带宽度由处理前的17 Hz～41 Hz扩展到处理后的17 Hz～65 Hz。处理前、后的效果图表明，本文方法能有效地提高地震资料的分辨率。

图8 二维实际资料处理前后频谱对比图Fig.8 Comparison of 2D actual seismic data spectrum before and after processing(a)原始剖面振幅谱；(b)处理后剖面振幅谱

图9 第277道处理前后的效果对比Fig.9 Comparison of the effects before and after processing(a)第277道处理前、后的信号；(b)第277道处理前、后的振幅谱

3.3 实际资料剖面的单道分析

抽取该二维地震资料的过井道，即第277道单道数据，用本文方法对其进行分析，处理前、后的效果如图9所示。观察图9(a)可知，处理后的整个波形趋势与原单道信号的波形基本保持一致，原来分开不明显的波形，经处理后波形有了明显的分开，信号的分辨率有了较大地提高。图9(b)显示了高分辨率处理前后振幅谱的对比，通过观察可知本文的方法提高了高频部分，拓宽了频谱的带宽。在原始剖面中，第277道单道数据高频振幅较弱，低频成分占主导地位，处理后该道数据高频信息显著增强，分辨率明显提高。

4 结论

针对大地滤波作用引起地震资料高频成分损失，导致分辨率降低的问题，笔者利用EMD自适应的时频处理方法，提出基于EMD的高分辨率地震资料处理方法，介绍了实现该方法的关键思路及详细过程。通过理论模型的测试和实际资料处理应用，经时域、频域的综合分析表明，该方法应用效果良好，具有理论的可行性和结果实现的可靠性。该方法的优势还在于能够保持原地震信号的主要特征不变的条件下，一些较弱或者缺失的层位信息得到了增强或者恢复，使地震资料的频带拓宽，高频成分提高，有效地提高地震资料的分辨率，有助于更好地识别薄层、薄互层、微小构造等，对目前越来越复杂的油气勘探具有借鉴意义。