地震信号线性与非线性时频分析方法对比

黄昱丞 郑晓东 栾 奕 杨廷强

(①中国石油勘探开发研究院油气地球物理研究所，北京 100083； ②香港中文大学理学院，香港 999077)

1 引言

时频分析是处理地震信号等非平稳信号的常用工具，现已广泛应用于沉积相带、岩性、储层和流体的检测分析。其基本思想是将一维时域地震信号表达为时频联合函数的形式，在时频二维空间进行描述，与传统的Fourier分析相比，在分析信号局部统计特征方面具有明显的优势。

时频分析方法在地球物理勘探领域的应用可上溯至二十世纪八、九十年代。1982年，Morlet等[1]提出用Gauss包络不同波长的零相位余弦子波描述层状介质下地震波传播与散射机制，通过构造复函数小波，在时频域开展地震信号变化特征的研究，其瞬时谱的思想直接触发了小波变换(WT)的萌芽。Okaya等[2]利用短时Fourier变换(STFT)对炮集数据的地表振动耦合引起的扫频假象进行时频滤波。Chakraborty等[3]则将连续小波变换(CWT)和匹配追踪分解(MPD)引入地震信号处理，获得了比STFT更高精度的时频谱，引起业界广泛关注。Grubb等[4]将离散小波变换用于地震勘探信号的分解与重构。李宏兵等[5]基于正演模型将WT用于含气砂岩的地震衰减检测。Partyka等[6]将时频分析正式引入地震资料解释领域，利用分频剖面成功地实现了地层厚度估算和可视化，形成了定性与半定量解释技术——谱分解(Spectral Decomposition)。Castagna等[7]提出基于地震瞬时谱直接检测烃类的技术。Liu等[8-10]则提出了基于Morlet子波和Ricker子波的MPD算法与谱均衡方法，拓宽频带，提高了数据分辨率，并将其用于河道识别与储层预测。Wang[11，12]将MPD应用于地震资料精细解释与烃类检测，并发展了多道MPD算法。Sinha等[13]提出了对CWT逆变换做Fourier变换以直接获得时频谱的时频连续小波变换，省去了CWT中尺度图与时频谱之间的繁琐转化。Pinnegar等[14]、高静怀等[15]、陈学华等[16，17]分别在S变换(ST)基础上提出了窗函数变化可控的参数化ST方法，统称为广义S变换(GST)，并将各自方法应用于地震相分析、沉积旋回的识别和油气检测。徐阳等[18]利用GST和离散WT进行面波干扰压制，同时可以减少对有效波的损伤。Li等[19]将平滑的Wigner-Ville分布(SWVD)法用于碳酸盐岩储层的地震资料衰减分析，Wu等[20]采用重排的平滑伪Wigner-Ville分布(RSPWVD)方法获得更为精确的时频谱用于低频阴影检测与储层展布刻画。Han等[21]利用第三代Hilbert-Huang变换(HHT)对地震信号进行谱分解，显示了自适应数据驱动类方法在地震数据分析中的有效性。曹思远等[22]利用改进的HHT结合小波包变换，消除了虚假信号分量，反映了信号的真实频率特性。薛雅娟等[23]基于集合经验模态分解方法(EEMD)结合WT进行本征模态函数的优选，研究了地震信号的衰减特征并提出了一套优化处理流程。刘晗等[24]应用同步挤压S变换进行地震信号时频分析，分频解释效果优于传统方法。Khonde等[25]较为系统地总结了谱分解方法在地球物理勘探领域的发展和应用现状，包括在三维地震数据解释[26]、直接烃类检测[27]、水合物勘探[28]、煤系地层岩浆岩侵入机理[29]等方面的应用。

总体而言，时频分析方法根据核函数表达形式的不同可分为线性、双线性和非线性三大类。线性方法包括Gabor展开[30]、STFT[31]、CWT[32]、ST[33]和GST[34-36]等。线性方法计算速度快，受限于不确定准则，时频分辨率有限，李振春等[37]曾对线性方法做过系统总结。双线性方法，如Wigner-Ville分布(WVD)[38]、伪Wigner-Ville分布(PWVD)、平滑伪Wigner-Ville分布(SPWVD)[39]等，归于Cohen类时频分布范畴。其中，WVD具有所有时频分析方法中最高的时频聚集性能，但受交叉项干扰，分辨率极差，后续发展的加窗改进方法则以牺牲部分分辨率为代价压制交叉项。除以上两类方法之外的都归为非线性方法，包括基于数据驱动的自适应分解方法，如自适应STFT[40，41]，根据信号特性灵活变化时窗长度，计算效率很高，分辨率较STFT有很大提高。Huang等[42]提出的经验模态分解(EMD)及其改进版本的集合经验模态分解(EEMD)、完全集合经验模态分解(CEEMD)[43，44]的核心均在于经验模态分解，尽管时频分辨率非常高，并采用加噪提纯技术消除模态混叠现象，但仍然存在边界效应等问题，且抗噪性较差，计算效率也不高。近年提出的融入EMD和CWT的同步挤压小波变换[45]，并以此发展出的一大类基于不同线性方法(STFT、ST)的“后处理”方法——同步挤压变换(SST)[46，47]，兼具较高的分辨率和计算效率并可实现信号重构，然而对低信噪比信号分析鲁棒性较差。基于信号稀疏表示的非线性方法，如包括基于Gabor原子的匹配追踪分解[48，49]、基于Chirplet变换(CT)[50，51]的自适应分解(Chirplet Decomposition，CD)[52]、指数调频小波变换[53]等，这类方法若匹配适当可以实现信号的稀疏分解，分辨率很高且免于交叉项干扰，但结果常常受控于时频原子库的架构，易出现原子弥散问题，参数较多，不易控制，且计算效率极低。还有衍生自CT的多项式调频小波变换、样条调频小波变换、泛谐波调频小波变换[54，55]、广义线性调频小波变换(GLCT)[56]等，这些新方法或通过挤压时频谱带宽，或通过调频核函数与时频脊线迭代拟合的方式，往往能够获得较高精度的信号时频表征，但一定程度上也受窗函数影响。近十年来还发展了如约束最小二乘谱分析[57]、基于反褶积的STFT[58]、自适应稀疏时频分解[59]、最大熵Wigner-Ville分布联合算法[60]等方法。

时频分析是一种将一维时域信号映射到二维时频空间描述和表达的方法，在信号局部或瞬态特征刻画方面具有独特的优势。目前使用的时频分析方法众多，使用条件和适用范围存在差异。Cas-tagna等[61]认为在地震解释领域时频分析方法没有对错之分，而是针对具体应用合不合适的问题。因此，比较它们之间的异同，分析参数选择的影响，对合理应用时频分析方法有实际价值。因此，系统地比较不同时频分析方法有助于改善地震时频分析在储层预测中的应用效果。

2 基本方法分析与参数选择

本文选取了包括线性方法中的STFT、CWT、ST、GST，双线性方法中的WVD、PWVD和SPWVD，以及非线性方法中的HHT、基于STFT的SST、CD和GLCT，总计三类11种方法，梳理各种方法的适用条件，研究参数选取的原则。上述各种时频分析方法的基本公式与控制参数如表1所示。

2.1 线性方法

线性方法是指信号变换满足线性叠加原理的时频分析方法[31]，包括STFT、WT、ST、GST等。

STFT[62]的本质是对信号加窗分段进行滑动谱分析，并假定信号在时窗内平稳。其时窗函数常用Gauss窗、Hamming窗等。采用Gauss窗可以获得信号最小的时宽和带宽乘积，此时STFT亦称作Gabor变换[30]。STFT的概念与定义所反映的物理意义简单明确，计算效率很高，但受固定窗函数影响导致分辨率有限，时窗长度是唯一制约时频表征的参数。实际应用中，根据输入信号特点和不同的分析需求权衡时间分辨率和频率分辨率。

表1 时频分析方法归类及定义汇总

WT在STFT基础上对时频网格进行改造，具有多分辨率特性[63]。由于地震信号分析中常用Gabor小波或Morlet小波，这类母小波不存在尺度函数，故一般以存在信息冗余的连续小波变换(CWT)的形式出现。在CWT表达式中，大尺度对应低频端，时间分辨率低，频率分辨率高；反之，小尺度对应高频端，时间分辨率高，频率分辨率低。为分析的直观性，一般需要将尺度转换为频率，它反映的是一个频段而非频点。良好的时频局域特性是CWT的优势所在，但若母小波选择不当，包括类型和峰值频率等因素，应用效果会大受影响。

ST继承和发展了STFT和CWT局部化的优点[33]，信号的ST是一种非严格意义上的CWT，若将ST看做采用特定窗函数的STFT，则相比Gauss窗函数，ST的窗函数相当于将尺度参数σ(标准差)表达为频率绝对值的倒数σ=1/|f|。与CWT类似，既满足了多分辨特性，同时又不需要从尺度到频率的转换，物理意义更加直观明了，计算效率很高，同时存在逆变换，可以实现信号滤波与重构。

虽然ST是小波变换的改进，但对信号局部时频特性的分析能力仍有待提升，因为ST采用固定的母小波，其时间窗形态随频率线性变化，使该方法在实际应用中仍受到一定限制。因此，众多学者提出了若干形态各异的时窗作为ST的基函数，统称为GST。本文采用陈学华[64]等提出的GST，即引入两个参数λ和p对ST的核函数进行改造，在ST的基础上，改写尺度参数σ=1/(λ|f|p)，其中λ>0、p>0。参数λ和p调节窗函数形状和变化趋势。λ和p越大，GST高频端时间分辨率越高，频率分辨率越低，且变化趋势对指数参数p更敏感，但调节过程中应避免高频失真的情况。因此，GST通过改变窗函数随频率变化的形态进一步改善分辨率，在实际应用中更为灵活。

总之，线性方法均受Heisenberg不确定性原理限制，无法同时在时域和频域获得最高的能量聚集性能。

2.2 双线性方法

双线性方法属于非线性方法范畴，满足二次叠加原理[31]。这类方法从能量角度刻画信号的时频分布，相比STFT的谱图或CWT的尺度图等线性方法取模值平方的方式，双线性方法更适于描述地震信号能量随时间频率的变化。为了避免负频率带来的交叉项干扰，双线性方法一般是把基于实信号Hilbert变换获得的解析信号作为输入项以获取单边时频谱。

所有的双线性方法中，WVD[38]因具有良好的边缘特性和极高的时频聚集性能而成为最基本、应用也最广泛的时频分布。解析信号z(t)的WVD定义为信号瞬时相关函数关于时延参数τ的Fourier变换。WVD对单分量平稳信号或者线性调频信号具有最高的时频分辨率，但其严重缺陷在于分析多分量信号时面临交叉项的干扰，也正由此派生了众多围绕抑制交叉项的双线性方法。

WVD可看做是窗函数为脉冲函数δ(t)的双线性分布，为了抑制交叉项，对时延参数τ加窗，实现对频率轴的平滑，由此衍生出PWVD[39]。本质上PWVD是以牺牲一定的时频聚集性能为代价，对WVD的结果进行了低通滤波。时窗越长，交叉项抑制越明显，但(频率)分辨率也越低。若τ时窗的窗宽为1，则退化为WVD。

事实上，若考虑同时对时间变量t和时延参数τ加窗平滑，亦即同时在时间和频率域进行平滑，可进一步抑制交叉项，从而出现了SPWVD[39]，但从一重积分变为二重积分的代价是计算量明显增加。具体而言，时间变量窗口越宽，每次循环参与计算的样点数越多，耗时越长。若时间变量窗口宽度为1，则退化为PWVD。

2.3 非线性方法

非线性方法包括：基于数据驱动的HHT、时频原子匹配追踪自适应分解以及同步挤压变换等参数化方法，它们都属于广义非线性方法范畴。

HHT是Huang等[42]提出的基于数据驱动的时频分析方法，其算法实现分为两步：①经验模态分解(EMD)，将数据信号分解为有限个本征模态函数(IMF)和背景趋势值或残差；②Hilbert谱分析，对这些IMF进行Hilbert变换求取瞬时频率和瞬时振幅，对瞬时频率重排，将信号映射为二维时频谱，亦称Hilbert谱。HHT可以获得高精度线谱，但边界效应、模态混叠和薄弱的数理基础成为HHT推广应用的困难。为避免模态混叠问题，本文采用集合经验模态分解(EEMD)[43]算法筛选IMF，该算法参数主要包括参与计算的白噪声集合数目NE和白噪声标准差Nstd。NE越大，理论上IMF越稳定，但耗时越长。对强时变信号，取较大的Nstd为宜，但一般不超过30%。

SST源自于Daubechies等[45]提出的一种“后处理”时频分析方法——同步挤压小波变换(Synchro-squeezed Wavelet Transform，SWT)。SWT的基本思想借鉴了HHT的EMD方法，将信号看做多个本征函数的叠加，在CWT基础上，进一步“挤压”瞬时频带，获得更精确的时频谱。基于小波系数对时间求导，初步估计瞬时频率。这样，将时间—尺度域转化为时间—频率域。沿用这一思路，将挤压操作与STFT、ST等方法结合，相继出现了同步挤压短时Fourier变换[65]、同步挤压S变换[46]等方法。将这类方法统称为“同步挤压变换”(Synchro-squeezing Transform，SST)。SST具有极高的时频聚集性，运算效率较高，但基于微分算子的挤压操作会引起数值不稳定，抗噪性低。本文采用基于STFT的SST进行信号分析，其优势在于时频网格均匀性和中低频段的信号时频特征刻画的精度和稳定性[66]。

CD源于Mann等[50]提出的Chirplet变换(CT)，在Gabor原子基础上加入线性调频相位项，实现在时频空间对非平稳信号的一阶逼近。基于MPD[48]的思想，利用预定义的归一化四参数Gaussian Chirplet原子库[67]，对实信号进行Hilbert变换得到解析信号，通过不断迭代匹配搜索最佳参数集合，满足时频原子在信号残差方向上投影最大化。经过若干次迭代满足判定条件，匹配终止并构建无交叉项WVD谱。CD在时频原子库选择恰当的情况下，可以实现信号的稀疏匹配，时频聚集性能极高，并免于交叉项干扰。本文使用的算法[68]需要预定义匹配分解的Chirplet原子个数NA，原子分解个数M控制CD分解拟合的精细程度，算法耗时正比于NA，故分量越多、变化越复杂的信号，处理耗时越长。

GLCT由Yu等[56]提出，其本质是在线性调频小波变换(Linear Chirplet Transform，LCT)的基础上进行调频斜率旋转优化。通过在每个时频样点上计算不同的调频斜率值对应的LCT，挑选其投影能量最大值作为该点的时频响应。实际操作中，需要对时频空间离散化，Chirplet原子旋转角度划分的稠密度N是关键参数。理论上N值越大，拟合局部时频特征越准确。GLCT在处理强时变、多分量信号方面具有一定优势，与传统方法相比，GLCT可以自适应拟合多分量非平稳信号时频特征，具有较高的时频分辨率，计算效率优于稀疏匹配类算法，且表征多分量信号时频谱时不受交叉项干扰。不足之处在于其本身受窗函数影响(默认采用Gauss窗)，且N的选取存在主观性，随着N值增大，计算量也显著增加，一般不宜超过5。

3 时频分辨率分析

分别采用上述方法对单分量谐波调频信号、多分量指数调频信号以及Ricker子波合成地震信号进行分析，前两类信号给出基于瞬时相位变化率定义的瞬时频率理论值作为参考，信号长度均为1024样点，采样间隔1ms。关于时频分辨率的优劣一般都是来自时频谱视觉效果上的直接比对，定量分析指标很少，但基于信息熵的时频聚集性检测则可以为此提供一个参考的定量指标。

Baraniuk等[69]提出利用Rényi熵估计信号所含信息量与时频谱复杂程度的思路。信号时频谱TFRs(t,f)的Rényi熵定义为

(1)

式中：α为Rényi熵的阶次，一般取3； 对数底b一般取2即可。Rényi熵是Shannon熵的广义形式，Shannon熵是Rényi熵在α→1时的极限，Rényi熵更具普适性。Rényi熵值越小，说明信号在时频域的能量复杂度越低，时频聚集性越高；反之，则复杂度越高，能量分布越均匀，聚集性越差。需要注意的是，聚集性和分辨率是两个概念，采用三阶Rényi熵仅作为刻画时频分析方法能量聚集性高低的一种定量手段，聚集性高并不说明分辨率一定高，也不一定反映信号真实的局部频率特征。

3.1 模型1：谐波调频信号

谐波调频信号的瞬时频率具有周期振荡特征(图1)，时频谱上呈现谐波形态，是一类典型的非平稳信号，即

s(t)=cos{2π[sin(30t)+60t]t}

(2)

信号的瞬时频率为

IF(t)=30cos(30t)+60

(3)

利用该信号考察各类方法对振荡调频瞬时频率曲线的时频聚集性能。由于瞬时频率振荡变化较为剧烈，故采用较短时窗处理。由图1可见，线性方法的分辨率普遍较低，WVD因交叉项的干扰出现频率假象。STFT、PWVD、SPWVD和GLCT等方法可以较准确地反映信号瞬时频率变化特征，但聚集性略差。Chirplet原子分解呈现时频特征的一阶逼近，故在拟合信号瞬时频率曲线峰谷点处误差较大。无噪声情况下，对于这种固定周期震荡调频信号的时频表征，HHT和SST具有最高的时频聚集性能。

图1 谐波调频信号、真实瞬时频率以及不同方法获得的时频谱

方法STFTCWTSTGSTWVDPWVDSPWVDHHTSSTCDGLCT熵4.98855.10514.89125.07694.72914.52154.62670.04621.71213.22114.7242

表2表明HHT和SST两种方法的时频聚集性是最高的，同时，其余双线性、非线性方法的时频聚集性能也都高于线性方法，这与时频谱时频分辨率一致。但WVD时频谱明显受交叉项干扰，Rényi熵却低于所有线性方法，也说明Rényi熵测度的定量判据肯定了WVD的高聚集性，但忽视了其交叉项干扰导致的低分辨率事实。

3.2 模型2：指数调频信号

参考Sinha等[13]的模型，信号由两个具有不同指数调频因子的扫频信号叠加而成，即

(4)

信号的两个瞬时频率分量为

(5)

两分量在0.6s前相互紧贴，随后瞬时频率变化趋势逐渐分开，信号后半段频率变化剧烈。

利用该信号有效频带范围远离Nyquist频率，考察各类方法的多分辨性能，并检验低频端时域分辨率，对各类方法性能提出了较高要求(图2)。由图2可见线性方法中只有GST可以较好地区分开两个分量并兼顾高低频端的分辨率。双线性方法受交叉项干扰严重，或者也因窗函数导致高频段分辨率降低。SST受窗函数影响，长时窗导致高频端分辨率较低。Chirplet原子分解因为原子匹配度不高，一阶拟合信号时频分量存在较大误差。HHT和GLCT两种方法可以高分辨地显示信号的两个扫频分量，但HHT在高频段存在IMF线谱振荡失稳的现象，而GLCT存在Chirplet时频原子旋转交叠的问题。总体而言，GST和HHT对该信号时频刻画效果最佳。

图2 指数调频信号、真实瞬时频率以及不同方法获得的时频谱

方法STFTCWTSTGSTWVDPWVDSPWVDHHTSSTCDGLCT熵3.87403.49374.01574.05104.79724.11414.02880.99062.26422.75364.3537

从Rényi熵角度看，HHT表现出了最高的时频聚集性能，GLCT算法本身在时频面每个样点能量都是非零的，故而Rényi熵较高，也说明了GLCT处理此类指数调频信号存在不足。

3.3 模型3：合成地震信号

参考Castagna等[7]的模型，给出不同峰值频率和时延的零相位Ricker子波合成的地震信号。图3给出了各单频子波及合成信号，可知有十二个零相位Ricker子波参与信号合成，包括两个40Hz子波、七个30Hz子波、两个20Hz子波和一个10Hz子波。该信号考察各类方法的多分辨性能，检验高频端频域分辨率亦即薄层(相邻子波)分辨率。

图3 合成地震信号、单频子波以及不同方法获得的时频谱

由图3可见，线性方法中STFT不能兼顾高低频，CWT、ST、GST等方法在低频端和高频端总体聚集性不高。相对而言，CWT、GST高频端可以区分相邻Ricker子波。SST呈现稀疏表示的特征，但对薄层的分辨率并不高。双线性方法中SPWVD表现最好，高低频分辨率大致保持，并具有较高的薄层分辨率，WVD、PWVD受交叉项干扰。HHT时频聚集性能极高，但从图中不能清晰区分子波各分量，分辨率不高。参数选取恰当时，CD可以提取出大部分子波分量，对低频的分辨能力较为突出，但薄层仍不能分辨。相比之下，GLCT则具有类似SPWVD的优异表现。

由于信号的稀疏性，HHT和SST的Rényi熵值均小于零，故这两种方法的时频聚集性最高，然而，两者对合成地震信号时频谱的刻画均出现较大程度的失真。

表4 合成地震信号不同方法时频谱的Rényi熵

4 计算效率分析

在XW9400型HP工作站(64位、AMD Opte-ron 2.8GHz双核处理器、RAM 8GB)上处理以上三个模型信号，计算时间如图4所示。所有的线性方法以及双线性方法中的WVD和PWVD处理1024个样点信号的耗时都在10-2～10-1s数量级，其余大部分非线性方法耗时都在1s以上，高出一到两个数量级。大部分方法计算效率受窗函数长度影响显著，如SPWVD在处理模型2时采用较宽的双重时窗上处理时间较多。此外，SST耗时主要体现在频带挤压操作中，但在非线性方法中效率相对较高；HHT在加噪提纯与IMF的筛选中耗时过长；CD处理时间正比于匹配的原子个数，若一阶拟合三个模型信号，分别预定义了10个、8个和12个原子，这种基于稀疏匹配的贪婪算法成为最耗时方法；GLCT耗时正比于原子旋转角度划分的稠密度N。总体而言，线性方法的计算效率远高于大部分非线性方法，非线性方法中诸如CD等方法是基于MPD算法思想，虽然对某些信号的时频估计效果较好，但在实际应用中计算成本较高。

图4 不同时频分析方法三个模型信号处理效率对比

一般而言，地震信号记录长度很有限，所以在单道试算时，时间频度表达式中被忽略的不同系数最高次幂和各低次幂项造成的时间消耗不容忽视。总体而言，各类方法的计算效率优劣对比明显(一到两个数量级)，实际计算耗时会因为海量的地震道数(如涉及叠前或三维地震数据的计算)而陡增，影响算法的实际应用，这也是应当考虑的现实问题。

5 抗噪性能分析

参考Yu等[56]的抗噪性能分析，对不同方法获得的时频谱|TFR(t,f)|进行瞬时峰值频率检测。对每一个样点t，瞬时峰值频率为

(6)

式中fc为信号的Nyquist频率。通过与理论瞬时频率曲线IF(t)相比较可计算IFe(t)的均方误差

(7)

以模型1谐波调频信号为例，对信号加入不同级别(无噪声，10dB，0dB)高斯白噪声，通过对比各类方法估计的瞬时峰值频率的均方误差考察其稳健性。WVD因交叉项干扰严重，而CD本身的不稳定性较高，两种方法检测结果均明显偏离真实值，故未参与测评。

不同的方法随着噪声级别增大(图5a)，瞬时频率估计值的标准偏差都不同程度地增大，其中CWT、SPWVD和GLCT的稳健性最高，是低信噪比情况下抗噪效果较好的方法；由不同方法的抗噪性对比(图5b)可见STFT、ST、GST、SST、PWVD和HHT随噪声级别增大估计偏差迅速增加，尤以GST(9.33%)、ST(7.61%)、HHT(6.82%)和SST(6.56%)受噪声干扰严重。对高信噪比(无噪声、加10dB白噪)数据，线性方法中的STFT(0.08%)和非线性方法中的GLCT(0.05%～0.09%)具有最高的稳健性；对低信噪比(0dB)数据，CWT(0.44%)、SPWVD(0.57%)和GLCT(1.84%)是首选方法。

综上所述，各类方法性能表现按“最低<低<一般<较高<高”排列分级，总结如表5所示。

图5 不同噪声级别下不同方法估计偏差对比(a)和不同方法在不同噪声级别下的稳健性变化(b)

方法窗函数交叉项聚集性低频分辨率薄层分辨率计算效率抗噪性STFT固定无低低低高较高CWT多分辨无低较高较高高高ST多分辨无低较高一般高低GST多分辨无低高高高低WVD无有高低低高-PWVD固定有较高高一般高一般SPWVD固定有较高高高一般高HHT无无高较高一般低低SST固定无高高较高较高低CD无无较高高低最低-GLCT固定无较高较高高一般高

6 应用实例分析

本文实际数据所在研究区位于塔里木盆地塔中地区，目的层段是上奥陶统良里塔格组碳酸盐岩，埋深约5000m，溶蚀孔隙极为发育，非均质性强[19]。工区三维地震数据采样间隔4ms，截取3000～4500ms时窗，层位mfs与sb3为目的层段顶、底界面，资料信噪比较低，主频不足20Hz，常规剖面(图6)难以识别储层形态与含油气特征。

选取STFT、CWT、SPWVD和GLCT作为代表，进行效果对比。图7给出了四种方法5Hz分频剖面，STFT剖面中反射轴不连续，也不能区分目的层顶、底；CWT低频端时窗较宽，因而时间分辨率很低，难以分辨目的层响应；SPWVD是基于能量的双线性分布，其浅层强反射掩盖了目的层信息，顶、底区分不开；GLCT显示了沉积层较多细节，目的层台地内顶底面反射被区分开并可以连续追踪，浅层地层反射连续性也更好。图8为四种方法获得的35Hz层内均方根振幅切片，在GLCT和SPWVD切片上台地边界(黄色箭头所示)以及潮汐水道(虚线圆框所示)更为清晰，台地内部水道(白色箭头所示)和泻湖沉积区前缘与内部断裂体系(黑色箭头所示)细节也更丰富，而STFT、CWT的切片效果相当，不能明显观察到这些地质特征。总体上，尽管信噪比较低，两类较稳健的非线性时频分析方法仍可以挖掘数据体中丰富的低频信息和平面细节，储层几何特征与空间展布也得到精细刻画。

图6 实际数据剖面

图7 实际数据5Hz不同方法分频剖面

图8 层间均方根振幅35Hz不同方法分频切片

7 结论

本文系统分析了用于地震勘探领域的11种时频分析方法，基于三个模型信号和碳酸盐岩储层实际地震数据比较了各类方法的时频分辨率、计算效率和抗噪性能。结论如下：

(1)线性类方法受限于不确定准则，时频聚集性能普遍较低，但计算效率很高，且不存在交叉项干扰；另外，STFT、CWT的稳健性较高，GST、GLCT的薄层分辨率较高； 非线性类方法时频聚集性能普遍较高，除SPWVD和GLCT外，抗噪性普遍较差；双线性方法中WVD、PWVD计算速度快但受交叉项干扰较严重，其余非线性类方法计算效率普遍不高。

(2)除了信号本身特征外，窗函数普遍影响各种方法的时频分辨率，而特定方法、特定的参数集合也会影响时频分辨率、计算效率和抗噪性能，如HHT(EEMD)参与计算的白噪信号的数目、Chirplet原子分解法中预定义的Chirplet原子数目、GLCT中时频面角度划分稠密度等参数对时频分辨率和计算耗时影响较大，前两种方法的实际应用中计算成本较高。

(3)碳酸盐岩储层实际地震资料分析显示，稳健性较高的非线性方法获得的分频剖面上可以更方便地追踪储层顶、底反射界面，分频切片上可以更清晰地反映潮汐水道和礁滩相带展布等特征。

(4)总体而言，低信噪比数据分析宜采用线性类方法，高信噪比数据宜采用时频分辨率更高的非线性类方法，有可能获得高精度的地层阻抗与局部频率变化信息。