专利名称：基于短时分数阶傅里叶变换的地震信号时频分解方法
技术领域：
本发明涉及一种基于短时分数阶傅里叶变换(Short Time Fractional FourierTransform, STFrFT)的地震信号时频分解方法。
背景技术：
在地震勘探中，当地震波在地下介质中传播时，其传播路径、振动强度和波形将随所穿过介质的弹性性质和几何形态发生复杂的变化。因此，地面将接收到经不同路径传播的P波、S波和较大振幅的发散面波以及各种噪声等信号成分，它们不仅到达时间不同，而且运动学和动力学特征也不同，并且经过了多次反射、折射和透射，加之地下介质对不同频率成分的吸收衰减也有差异。因此，地震信号是典型的非平稳信号，其频谱成分及信号的各种统计特性随时间发生显著变化。在地震资料处理中，如何应用信号处理方法，检测地震信 号的振幅、频率和相位等多种信号特征参数的不稳定变化，对于准确地确定地下界面反射 波的到达时间、获取地震反射界面的位置、提高地震资料分辨率、精细刻画地下的地质构造等，具有非常重要的指导意义。地震谱分解是指对地震道进行连续时频分析获得地震的频谱(振幅谱及相位谱)。谱分解技术不仅可以提高地震资料对薄储层的解释预测能力，而且能从常规宽频地震数据体中提取出更丰富的地质信息，提高地震资料对特殊地质体的解释识别能力。因此，这项新的地震属性分析技术一经推出便引起业界广泛关注，很快就成为地震勘探技术发展的热点。1947年，R. K. Potter等首次提出了一种实用的时频表示方法-短时Fourier变换(short-time Fourier transform, STFT,又名 time-dependent Fourier transform,或windowed Fourier transform),并将其绝对值的平方称为“声音频谱图’，此即为后来所说的谱图(spectrogram)。20世纪60年代中期,Cohen发现众多的时频分布只是Wigner-Ville分布的变形，可以用统一的形势表示，习惯称之为Cohen类时频分布。之后，Jones等提出的数据自适应最优核时频分布等。这类方法实际上是从提高分辨率和抑制交叉干扰项的目的出发。1982年，小波变换线性时频表示，其创造性思想是由法国地球物理学家J. Morlet提出的，后经其他几位法国学者的再塑造，使之成了一种基础坚实、应用广泛的信号分析工具。1996年，R. G. Stockwell等人在小波变换的基础上，提出了 S变换。2007年，YanghuaWang发表文章，“通过匹配追踪方法的地震时频谱分解”。地震道可以被分解成一系列子波，这里的子波是通过匹配追踪算法和时频信号进行匹配计算得到的，每个子波选择的重复过程需要在大的并且冗余的时频信号字典中进行。2008年，陈学华等人对S变换改进后得到一种新的广义S变换，用于分析和补偿地震信号的高频成分，得到了精细的地震层序识别剖面。实际资料处理表明，它对砂岩油气藏的成层特征的分析具有分辨率高和高信噪比的优点。地震信号的广义S变换时频谱分解可作为地震层序识别和解释的重要补充。
以上这些方法都有它的局限性和适应性，都受实际地震资料和地下地质条件的制约，谱分解技术也不例外。与本发明相关的现有技术包括时频特征表示的任务是描述信号的频谱含量在时间上的变化情况，研究并了解时变频谱在数学和物理方面的概念，最终目的是建立一种分布，以便能在时间和频率域上同时表示信号的能量或者强度，并对其进行分析和处理。时频表示方法按照时频联合函数的不同分为线性表示和双线性时频表示两种。基于在时间和频率均局域化的基本函数(亦称“时频原子”或“原子”)分解的线性方法，包括短时傅利叶变换、小波变换等。双线性时频表示也称作二次型时频表示，主要有Cohen类时频分布和仿射类(Affine)双线性时频分布，其中最著名的是Wigner-Ville分布。此类方法均在分辨率和交叉项干扰之间取得某种折衷。以下将介绍现有相关的时频分解方法短时傅里叶变换、魏格纳分布。
I :短时傅里叶变换尽管傅里叶变换及其离散形式DFT已经成为信号处理，尤其是时频表示中最常用的工具，但是在信号处理，尤其是非平稳信号处理过程中，特别是地震信号，人们经常需要对信号的局部频率以及该频率发生的时间段有所了解。由于标准傅里叶变换只在频域有局部分析的能力，而在时域内不存在局部分析的能力，故Dennis Gabor于1946年引入短时傅里叶变换(Short-Time Fourier Transform) 短时傅里叶变换的基本思想是把信号划分成许多小的时间间隔，用傅里叶变换分析每个时间间隔，以便确定该时间间隔存在的频率。设信号f(x)，并假定该信号在一个以时间τ为中心且范围有限的窗口函数gU-τ)内是稳定的，这样窗口函数的傅里叶变换就定义为短时傅里叶变换
_3] f{co,T) = ^J\x)gT{x)e~,xo,dx( -l)如果gT (x)为方波函数
I r I~r,x G It
=水I(1-2)
、0,其他其中| τ|表示Ιτ的长度。其中R表示整个实轴。从上述公式很容易看出，为了分析信号f(x)在时刻τ的局部频域信息，上式实质上是对函数f(x)加上窗函数g,(x)。显然，窗口的长度| τ I越小，则越能够反映出信号的局部频域信息。2:魏格纳一威利分布魏格纳(Wigner)分布是在定性方面不同于频谱图的一些分布的原型，发现它的长处和短处已经成为这个领域研究的主要动向。魏格纳分布(Wigner Distribution,简称WD)是由威利(Ville)引进信号分析的，大约在魏格纳发表论文15年以后(1959年)。维尔对魏格纳分布给出了一个似乎合理的论证，并根据特征函数方法推导得出了魏格纳分布。值得注意的是，同样类型的推导大约在同一时间Moyal也使用了。魏格纳威利变换(WVD)定义如下信号S⑴和它的频谱S ( ω )的WVD是
IF(/,iy) =—(t--T)s(t + —r)e jreWr(13)= 士 JY 一
Ιπ·122(1-4)这两个表达式的等价性通过用频谱表示信号，很容易验证。魏格纳威利分布是作为信号的双线性表示的，因为信号在其计算中两次出现。考虑公式(1-3)和(1-4)中不含有任何的窗函数，因此避免短时傅立叶变换时间分辨率与频率分辨率相互牵制的矛盾，它的时间一带宽达到了测不准原理给出的下界。但是魏格纳一威利分布本质不是线性的，即两信号和的WVD分布并不等于每一个信号的WVD分布之和。令 x(t) = X1 (t)+X2 (t),则:(^Ω) = J[X1 (t + ~) + x2{t + —)]*[X1 {t~—) + x2 (t — —)]e ,n'dr
J 2 2 2 2 (1-5)= ^ (/,n) + ^,(/,n) + 2ReF, ..J/,Ω)]( !_6)Κφ2Ι14『Χι+τ2( ，Ω)^Χια) X2 (t)的互 WVD，称之为交叉项。由公式(1-6)可以看到有时魏格纳分布在时间和频率上把这些值置于两个信号的中间；有时这些值又处在时频平面和所期望的成分争夺位置。因此产生了交叉项。交叉项极大地干扰了时频分布，同时也抑制了二次型时频分布的推广。正如上所述魏格纳威利分布对于单分量信号有很好的能量聚集性，但对于多分量信号，由于其固有的双线性性质，使得这种时频分解存在严重的交叉项，对于地震信号中一些较弱的有用分量检测是不利的，短时傅里叶变换(Short-Time Fourier Transform)能够避免交叉项的出现，但是在低信噪比地震信号条件下分解信号频谱效果不理想。
发明内容
为了克服现有技术的上述缺点，本发明提供了一种基于短时分数阶傅里叶变换的地震信号时频分解方法，选用高斯窗函数，通过自适应调整窗口的宽度，获得较高的频率分辨率，同时利用分数阶核函数角度旋转的特点，在最优旋转角时得到地震信号的最佳能量聚集，消除交叉项，对频率的分辨率越高，定位越精确，则更有利于识别特定的地质结构，所以该方法解决了现有的地震信号时频分解方法的不足。本发明解决其技术问题所采用的技术方案是一种基于短时分数阶傅里叶变换的地震信号时频分解方法，包括如下步骤(I)获取一维地震信号；(2)选择高斯窗函数；(3)选择一个初始α值来固定变换核，进行短时分数阶傅里叶变换；(4)调整α值，生成新的变换核，进行下一次短时分数阶傅里叶变换，直至完成11。次分解；(5)对1次分解结果做主成分分析，得到N个时频分解谱；(6)将N个时频分解谱作为最终的时频分解结果，得到大小为N*nt*nf的三维数据体。所述高斯窗函数的表达式为g(x)= —g 4σ,σ>0
2-^πσ
νο所述短时分数阶傅里叶变换是指将信号在时间轴逆时针旋转角度α后在分数阶时-频域的投影。与现有技术相比，本发明的积极效果是克服了传统算法如短时傅里叶和魏格纳威利分布等的缺点，对地震信号频谱分析有重大意义，具体表现如下I)克服了短时傅里叶变换等一类算法对低信噪比地震信号时频分解效果不理想的状况。2)通过旋转局部时频面到最优α值，克服了传统方法的频谱交叉项干扰，使分解后的频谱更精确可信。3)该算法采用PCA对多次时频分解结果进行降维处理，提取出主要信息，简化了对结果进行分析的工作量。4)该算法在提升频率分解效果的同时并未明显增大计算复杂度，具有很高的实用性。
本发明将通过例子并参照附图的方式说明，其中图I是地震信号短时分数阶傅里叶变换流程图；图2是短时傅里叶变换STFT的矩形支撑边界； 图3是STFrFT的平行四边形支撑边界。
具体实施例方式一种基于短时分数阶傅里叶变换的地震信号时频分解方法，如图I所示，包括如下步骤( I)获取一维地震信号地震信号是一种典型的非平稳信号，地震信号一般以地震道为单位进行组织，采用 SEG-Y 文件格式存储。SEG-Y 格式是由 SEG (Society of Exploration Geophysicists)提出的标准磁带数据格式之一，它是石油勘探行业地震数据的最为普遍的格式之一。典型的地震道包含某些低频噪音，例如面波，以及某些高频环境噪音；有用的地震反射能量常常限于10-70HZ，而优势频率在30Hz周围。该频带范围基本属于中低频，所以需要较高的频率分辨率来分析；同时特定的窄频带又反映特定的地质结构，就要尽量避免交叉项的出现；而STFrFT方法满足上述要求。选定一道数据后，截取需要分析的时间段即取得可用的一维地震信号。(2)选择高斯窗函数窗函数是局部时频分析的有力工具，不同的窗函数在频率分辨率，频带展布上有不同表现，本方案采用高斯窗函数，其窗宽度可调且具有较高的频率分辨率，适合地震信号的时频分解。高斯窗表达式为
I -兑g(x) = ~j=e -1. σ > O
2伽(2-1)(3)选择一个初始α值来固定变换核，进行短时分数阶傅里叶变换短时分数阶傅里叶变换(STFrFT)是一种时频分析工具，信号的STFrFT可看成将信号在时间轴逆时针旋转角度α后在分数阶时-频域的投影。(a)对于信号S(U)的分数阶傅里叶变换(FrFT)定义为
权利要求
1.ー种基于短时分数阶傅里叶变换的地震信号时频分解方法，其特征在于包括如下步骤 (1)获取ー维地震信号； (2)选择高斯窗函数； (3)选择ー个初始α值来固定变换核，进行短时分数阶傅里叶变换； (4)调整α值，生成新的变换核，进行下一次短时分数阶傅里叶变换，直至完成ηα次分解； (5)对^次分解结果做主成分分析，得到N个时频分解谱； (6)将N个时频分解谱作为最終的时频分解结果，得到大小为N*nt*nf的三维数据体。
2.根据权利要求I所述的基于短时分数阶傅里叶变换的地震信号时频分解方法，其特征在于所述高斯窗函数的表达式为g(x) = ~^e—
3.根据权利要求I所述的基于短时分数阶傅里叶变换的地震信号时频分解方法，其特征在于所述短时分数阶傅里叶变换是指将信号在时间轴逆时针旋转角度α后在分数阶时-频域的投影。
全文摘要
本发明公开了一种基于短时分数阶傅里叶变换的地震信号时频分解方法，选用高斯窗函数，通过自适应调整窗口的宽度，获得较高的频率分辨率，同时利用分数阶核函数角度旋转的特点，在最优旋转角时得到地震信号的最佳能量聚集，消除交叉项,对频率的分辨率越高，定位越精确，则更有利于识别特定的地质结构，所以该方法解决了现有的地震信号时频分解方法的不足。本发明的积极效果是克服了传统算法如短时傅里叶和魏格纳威利分布等的缺点，对地震信号频谱分析有重大意义。
文档编号G01V1/30GK102798891SQ201210299230
公开日2012年11月28日 申请日期2012年8月22日 优先权日2012年8月22日
发明者钱峰, 黄佳, 胡光岷 申请人:电子科技大学