基于短时分数阶傅里叶变换的地震信号时频分解方法

摘要

本发明公开了一种基于短时分数阶傅里叶变换的地震信号时频分解方法，选用高斯窗函数，通过自适应调整窗口的宽度，获得较高的频率分辨率，同时利用分数阶核函数角度旋转的特点，在最优旋转角时得到地震信号的最佳能量聚集，消除交叉项,对频率的分辨率越高，定位越精确，则更有利于识别特定的地质结构，所以该方法解决了现有的地震信号时频分解方法的不足。本发明的积极效果是：克服了传统算法如短时傅里叶和魏格纳威利分布等的缺点，对地震信号频谱分析有重大意义。

说明书

技术领域

本发明涉及一种基于短时分数阶傅里叶变换(Short Time Fractional Fourier  Transform，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年，Yanghua Wang发表文章，“通过匹配追踪方法的地 震时频谱分解”。地震道可以被分解成一系列子波，这里的子波是通过匹配追 踪算法和时频信号进行匹配计算得到的，每个子波选择的重复过程需要在大的 并且冗余的时频信号字典中进行。

2008年，陈学华等人对S变换改进后得到一种新的广义S变换，用于分 析和补偿地震信号的高频成分，得到了精细的地震层序识别剖面。实际资料处 理表明，它对砂岩油气藏的成层特征的分析具有分辨率高和高信噪比的优点。 地震信号的广义S变换时频谱分解可作为地震层序识别和解释的重要补充。

以上这些方法都有它的局限性和适应性,都受实际地震资料和地下地质条 件的制约,谱分解技术也不例外。

与本发明相关的现有技术包括：

时频特征表示的任务是描述信号的频谱含量在时间上的变化情况，研究并 了解时变频谱在数学和物理方面的概念，最终目的是建立一种分布，以便能在 时间和频率域上同时表示信号的能量或者强度，并对其进行分析和处理。时频 表示方法按照时频联合函数的不同分为线性表示和双线性时频表示两种。基于 在时间和频率均局域化的基本函数（亦称“时频原子”或“原子”）分解的线 性方法，包括短时傅利叶变换、小波变换等。双线性时频表示也称作二次型时 频表示，主要有Cohen类时频分布和仿射类（Affine）双线性时频分布，其中 最著名的是Wigner-Ville分布。此类方法均在分辨率和交叉项干扰之间取得某 种折衷。以下将介绍现有相关的时频分解方法：短时傅里叶变换、魏格纳分布。

1：短时傅里叶变换

尽管傅里叶变换及其离散形式DFT已经成为信号处理，尤其是时频表示 中最常用的工具，但是在信号处理，尤其是非平稳信号处理过程中，特别是地 震信号，人们经常需要对信号的局部频率以及该频率发生的时间段有所了解。 由于标准傅里叶变换只在频域有局部分析的能力，而在时域内不存在局部分析 的能力，故Dennis Gabor于1946年引入短时傅里叶变换（Short-Time Fourier  Transform）。短时傅里叶变换的基本思想是：把信号划分成许多小的时间间隔， 用傅里叶变换分析每个时间间隔，以便确定该时间间隔存在的频率。设信号 f(x),并假定该信号在一个以时间τ为中心且范围有限的窗口函数g(x-τ)内是 稳定的,这样窗口函数的傅里叶变换就定义为短时傅里叶变换：

$\hat{f}\left(\overline{\omega ,\tau }\right)={\int }_{R}f\left(x\right)g_{\tau }\left(x\right)e^{-\mathrm{ix}\bar{\omega }}\mathrm{dx}---(1-1)$

如果g_τ(x)为方波函数:

其中|I_τ|表示I_τ的长度。其中R表示整个实轴。从上述公式很容易看出，为 了分析信号f(x)在时刻τ的局部频域信息，上式实质上是对函数f(x)加上窗函 数g_τ(x)。显然，窗口的长度|I_τ|越小，则越能够反映出信号的局部频域信息。

2：魏格纳―威利分布

魏格纳（Wigner）分布是在定性方面不同于频谱图的一些分布的原型，发 现它的长处和短处已经成为这个领域研究的主要动向。魏格纳分布（Wigner  Distribution，简称WD）是由威利（Ville）引进信号分析的，大约在魏格纳发 表论文15年以后（1959年）。维尔对魏格纳分布给出了一个似乎合理的论证， 并根据特征函数方法推导得出了魏格纳分布。值得注意的是，同样类型的推导 大约在同一时间Moyal也使用了。魏格纳威利变换（WVD）定义如下：

信号s(t)和它的频谱S(ω)的WVD是:

$W(t,\omega )=\frac{1}{2\pi }\int s^{*}(t-\frac{1}{2}\tau )s(t+\frac{1}{2}\tau )e^{-\mathrm{j\tau \omega }}\mathrm{d\tau }---(1-3)$

$=\frac{1}{2\pi }{\int S}^{*}(\omega +\frac{1}{2}\theta )S(\omega -\frac{1}{2}\theta )e^{-\mathrm{jt\theta }}\mathrm{d\theta }---(1-4)$

这两个表达式的等价性通过用频谱表示信号，很容易验证。魏格纳威利分 布是作为信号的双线性表示的，因为信号在其计算中两次出现。

考虑公式（1-3）和（1-4）中不含有任何的窗函数，因此避免短时傅立叶 变换时间分辨率与频率分辨率相互牵制的矛盾，它的时间－带宽达到了测不准 原理给出的下界。但是魏格纳－威利分布本质不是线性的，即两信号和的 WVD分布并不等于每一个信号的WVD分布之和。令x(t)＝x₁(t)+x₂(t)，则：

$W_x(t,\Omega )=\int [x_1(t+\frac{\tau }{2})+x_2(t+\frac{\tau }{2})]*[{x_1}^{*}(t-\frac{\tau }{2})+{x_2}^{*}(t-\frac{\tau }{2})]e^{-\mathrm{j\Omega t}}\mathrm{d\tau }---(1-5)$

$=W_{x_1}(t,\Omega )+W_{x_2}(t,\Omega )+2\mathrm{Re}\left[W_{x_1+x_2}(t,\Omega )\right]---(1-6)$

式中是x₁(t)和x₂(t)的互WVD，称之为交叉项。由公式（1-6） 可以看到：有时魏格纳分布在时间和频率上把这些值置于两个信号的中间；有 时这些值又处在时频平面和所期望的成分争夺位置。因此产生了交叉项。交叉 项极大地干扰了时频分布，同时也抑制了二次型时频分布的推广。

正如上所述：魏格纳威利分布对于单分量信号有很好的能量聚集性，但对 于多分量信号，由于其固有的双线性性质，使得这种时频分解存在严重的交叉 项，对于地震信号中一些较弱的有用分量检测是不利的，短时傅里叶变换 (Short-Time Fourier Transform)能够避免交叉项的出现，但是在低信噪比地震信 号条件下分解信号频谱效果不理想。

发明内容

为了克服现有技术的上述缺点，本发明提供了一种基于短时分数阶傅里叶 变换的地震信号时频分解方法，选用高斯窗函数，通过自适应调整窗口的宽度， 获得较高的频率分辨率，同时利用分数阶核函数角度旋转的特点，在最优旋转 角时得到地震信号的最佳能量聚集，消除交叉项,对频率的分辨率越高，定位 越精确，则更有利于识别特定的地质结构，所以该方法解决了现有的地震信号 时频分解方法的不足。

本发明解决其技术问题所采用的技术方案是：一种基于短时分数阶傅里叶 变换的地震信号时频分解方法，包括如下步骤：

（1）获取一维地震信号；（2）选择高斯窗函数；（3）选择一个初始α值来 固定变换核，进行短时分数阶傅里叶变换；（4）调整α值，生成新的变换核,进 行下一次短时分数阶傅里叶变换,直至完成n_α次分解；（5）对n_α次分解结果做主 成分分析，得到N个时频分解谱；（6）将N个时频分解谱作为最终的时频分解结 果，得到大小为N*n_t*n_f的三维数据体。

所述高斯窗函数的表达式为：

所述短时分数阶傅里叶变换是指将信号在时间轴逆时针旋转角度α后在分 数阶时-频域的投影。

与现有技术相比，本发明的积极效果是：克服了传统算法如短时傅里叶和 魏格纳威利分布等的缺点，对地震信号频谱分析有重大意义，具体表现如下：

1）克服了短时傅里叶变换等一类算法对低信噪比地震信号时频分解效果 不理想的状况。

2）通过旋转局部时频面到最优α值，克服了传统方法的频谱交叉项干扰， 使分解后的频谱更精确可信。

3）该算法采用PCA对多次时频分解结果进行降维处理，提取出主要信息， 简化了对结果进行分析的工作量。

4）该算法在提升频率分解效果的同时并未明显增大计算复杂度，具有很 高的实用性。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是地震信号短时分数阶傅里叶变换流程图；

图2是短时傅里叶变换STFT的矩形支撑边界；

图3是STFrFT的平行四边形支撑边界。

具体实施方式

一种基于短时分数阶傅里叶变换的地震信号时频分解方法，如图1所示，包 括如下步骤：

（1）获取一维地震信号

地震信号是一种典型的非平稳信号，地震信号一般以地震道为单位进行组 织，采用SEG-Y文件格式存储。SEG-Y格式是由SEG（Society of Exploration  Geophysicists）提出的标准磁带数据格式之一，它是石油勘探行业地震数据的 最为普遍的格式之一。典型的地震道包含某些低频噪音,例如面波,以及某些高 频环境噪音;有用的地震反射能量常常限于10-70Hz,而优势频率在30Hz周围。该 频带范围基本属于中低频，所以需要较高的频率分辨率来分析；同时特定的窄 频带又反映特定的地质结构，就要尽量避免交叉项的出现；而STFrFT方法满足 上述要求。选定一道数据后，截取需要分析的时间段即取得可用的一维地震信 号。

（2）选择高斯窗函数

窗函数是局部时频分析的有力工具，不同的窗函数在频率分辨率，频带展 布上有不同表现，本方案采用高斯窗函数，其窗宽度可调且具有较高的频率分 辨率，适合地震信号的时频分解。

高斯窗表达式为：

$g\left(x\right)=\frac{1}{2\sqrt{\mathrm{\pi \sigma }}}{e}^{-\frac{x^2}{4\sigma }},\sigma >0---(2-1)$

（3）选择一个初始α值来固定变换核，进行短时分数阶傅里叶变换

短时分数阶傅里叶变换（STFrFT）是一种时频分析工具，信号的STFrFT可 看成将信号在时间轴逆时针旋转角度α后在分数阶时-频域的投影。

(a)对于信号s(u)的分数阶傅里叶变换（FrFT）定义为：

$S_{\alpha }\left(\upsilon \right)=F^{\alpha }\left[s\left(u\right)\right]=\underset{-\infty }{\overset{+\infty }{\int }}s\left(u\right)K_{\alpha }(u,\upsilon )\mathrm{du}---(2-2)$

(b)其中FrFT变换核为:

$K_{\alpha }(u,\upsilon )=\left(\begin{array}{c}\sqrt{\frac{1-j\cot \alpha }{2\pi }}\exp (j\frac{u^2+{\upsilon }^2}{2}\cot \alpha -\mathrm{ju\upsilon }\csc \alpha ),a\ne \mathrm{n\pi }\\ \delta (u-\upsilon ),a=2\mathrm{n\pi }\\ \delta (u+\upsilon ),a=(2n\pm 1)\pi \end{array}\right)---(2-3)$

由图2与图3可以看出，短时分数阶傅里叶变换(STFrFT)的窗函数的支撑边 界是由短时傅里叶变换STFT的窗函数支撑边界的频率轴逆时针旋转角度a得到 的。

(c)当变化核基函数为窗函数K_α(u,υ)时，该变换即称为短时分数阶傅里叶 变换(STFrFT):

$\mathrm{STFrFT}(t,f)=\sqrt{\frac{1-j\cot \alpha }{2\pi }}\underset{-\infty }{\overset{+\infty }{\int }}s\left(u\right)·g(u-t)·\exp \left[j\frac{(u^2+f^2)\cos \alpha -2f·u}{2\sin \left(\alpha \right)}\right]\mathrm{du},a\ne \mathrm{n\pi }---(2-4)$

其中n∈Z

$g\left(u\right)=\frac{1}{\sqrt{2\pi {{\sigma }_t}^2}}\exp [-{\left(\frac{u}{\sqrt{2}{\sigma }_t}\right)}^2]---(2-5)$

（4）调整α值，生成新的变换核,进行下一次短时分数阶傅里叶变换（即下 一次时频分解）,直至完成n_α次分解。

（5）对n_α次分解结果做主成分分析(PCA)，得到N个时频分解谱：

当做完前四步以后,会生成一个n_α*n_t*n_f的三维数据体（n_α表示取不同α值 的分解次数，n_t表示一维信号时间采样点个数，n_f表示时频分解频率点个数）, 数据量比较大,人为地从n_α次时频分解里面选择最优分解是很难做到的,我们需 要一种自动的方法从三维数据体里面选择比较典型、有代表性的时频分解，主 成分分析(PCA)可以满足这种要求。

主成分分析(Principal Component Analysis,PCA)：是一种提取事物主 要成分的统计分析方法，它可以从多元事物中解析出主要影响因素，揭示事物 的本质，简化复杂的问题。计算主成分的目的是将高维数据投影到较低维空间， 在降低数据维数的同时保留数据携带的主要信息。这样我们仅需人工观察PCA方 法处理过后的少量时频分解结果，就可以得到n_α*n_t*n_f三维数据体的绝大部分 信息，下面给出具体步骤：

（a）将三维数据体表示成关于所有分解次数j，其中j={1,2,…n_α}的n_α*n_α的互相关协方差矩阵C

$C=\left(\begin{array}{ccccc}c(1,1)& ···& c(1,k)& ···& c(1,n_a)\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ c(j,1)& ···& c(j,k)& ···& c(j,n_a)\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ c(n_a,1)& ···& c(n_a,k)& ···& c(n_a,n_a)\end{array}\right)---(2-6)$

其中j行k列的值C_jk为：$C_{\mathrm{jk}}=\underset{n=1}{\overset{n_t}{\Sigma }}\underset{m=1}{\overset{n_f}{\Sigma }}{d}_{\mathrm{mn}}^{\left(j\right)}{d}_{\mathrm{mn}}^{\left(k\right)}---(2-6),$其中与表示第j次与第k次分解时在时间为m，频率为f的位置得到的振幅值。

（b）求解上述协方差矩阵C,使得：

Cv_p＝λ_pv_p        (2-8)

其中λ_p表示该协方差矩阵的特征值，个数为n_α，按照降序排列，v_p表示λ_p对 应的特征向量，每一个v_p为n_α*1大小的数组，个数为n_α。

（c）将n_α次分解的结果投影到前N个主特征向量上，对应生成N个新的主要 的且有代表性的时频分解谱,其中第p个谱时间为m，频率为n的位置对应的系数 为：

$A_{\mathrm{mn}}^{\left(p\right)}=\underset{j=1}{\overset{n_a}{\Sigma }}{v}_p^{\left(j\right)}{d}_{\mathrm{mn}}^j---(2-9)$

N的选取依据是：前N个谱的能量和至少要占到总能量的80%以上，一般情况 下N到4即可。

（6）取第（5）步得到的PCA降维后的N个谱作为最终的时频分解结果（大小 为N*n_t*n_f的三维数据体）。