CN104932012A - 一种地震信号的分数域局部功率谱计算方法 - Google Patents

Abstract

一种地震信号的分数域局部功率谱计算方法，其主旨在于提高功率谱图的时频特性。其方案为：对输入信号做p阶分数傅里叶变换，求该阶次下的时间带宽积，并对变换阶次p遍历；找到信号分数域时间带宽积最小值对应的阶次p₀及时宽、频宽；根据时宽和频宽构造高斯窗函数，再进行p₀阶分数傅里叶变换得最优窗函数；将信号左端和右端进行周期性延拓，用最优窗从信号第一个点为中心开始截取固定长度(即信号本身的长度)得到原信号的局部，对局部进行现代谱估计(如AR模型法、Burg法等)，得到信号的局部功率谱；再将窗中心逐点移动，即可获得有别于传统谱图的局部功率谱时频分布。

Description

一种地震信号的分数域局部功率谱计算方法

技术领域

本发明属于地震信号处理领域。方法基于最优分数阶傅里叶变换，将现代功率谱估计与分数域时频分析技术相结合，提出了新的局部功率谱计算方法。通过该方法可以提高信号功率谱图的时频分辨率，可以在地震信号处理中获得高分辨率的局部功率谱特征，为后续储层流体预测和识别创造有利条件。

技术背景

分数阶傅里叶变换起源于1929年Wiener的研究，兴起于1980年Namias的研究，当时他利用特征值的任意次幂运算明确的提出了阶数为分数的傅里叶变换的概念。20世纪90年代以后，分数傅里叶变换作为一种新型的信号分析工具受到越来越多的关注，并在光学领域最先得到广泛应用。1994年Almeida揭示了分数傅里叶变换和传统时频分析工具的关系，得出了分数傅里叶变换可以解释为信号在时频面上坐标轴绕原点逆时针旋转的重要结论。至此分数傅里叶变换被赋予明确的物理意义，吸引着越来越多的研究者参与相关研究。分数阶傅里叶变换较传统傅里叶变换的一个显著优点就是能为时频分析提供更大的选择余地，以获得某些性能上的改进。

功率谱描述了随机过程中各频率成分平均功率的大小，在工程上广泛应用于随机信号处理，是随机信号的重要表征形式之一。它源于1898年Schuster提出的周期图概念。20世纪50年代后，功率谱估计方法发展迅速，间接法、AR模型法、Burg谱估计法、特征空间分解法等相继出现。到现在信号功率谱密度估计方面已有许多新的算法，这些新方法连同演变出来的各种算法不下几十种。这些方法作为一种功率谱估计手段在生物医学、语音、雷达、地震勘探等领域得到了广泛应用。

然而功率谱估计给出的是一种基于信号全局特征的频谱信息，这只对平稳信号产生较好的分析效果。对于非平稳信号，在不同的时间段包含的频谱信息是不一样的，若对信号直接做功率谱分析，就无法得到信号频谱随时间变化的特征。1946年Gabor提出了一种对信号局部加窗后再做傅里叶变换的方法即Gabor变换。这种变换能够很好的反映出信号局部的频谱信息，于是得到了快速发展。随后短时傅里叶变换、Wigner-Vile分布、Cohen分布、小波变换、S变换等时频分 析技术相继出现。功率谱图与短时傅里叶变换诞生于同一时期，它可由短时傅里叶变换幅值的平方求得，它是一种不同于短时傅里叶变换、Wigner-Vile分布、小波变换等，且能够反映信号在时间和频率上的能量分布特性的一种优秀的时频分析技术。

在地震勘探领域，利用地震资料直接进行储层流体识别，已成为地震储层预测技术的主流发展方向之一。频谱分解技术是其中一种行之有效的技术，它主要利用不同性质的流体对地震波的吸收、衰减不同这一特点，通过地震信号不同层位下的频率属性来反映地下传输介质的物理特性，从而为判断该层位的油气特性提供依据。地震信号作为一种非平稳信号，频谱分解的核心是时频分析技术。现有的时频分析技术因为时频分辨率、交叉项等问题制约着其在地震信号处理中的应用效果，同时也激发了人们寻求新的时频分析方法的热情。

2002年Lufiye等人将分数域傅里叶变换的思想引入到时频分析中，有效地提高了信号的时频分析效果。分数域时频分析在传统时频分析的基础上利用时频旋转特性，通过判定某一旋转角度下信号具有的最小时间带宽积来达到最优的时频分辨率。2013年陈颖频等将基于分数阶的自适应最优阶Gabor变换应用到地震信号频谱分解中，取得了较好的效果。2015年田琳等利用反卷积方法求取了分数阶Gabor谱图并应用到地震储层频谱分解中。该方法仍可以看作是在直接法计算功率谱的概念下做的进一步处理，即现阶段的分数域时频分析技术尚未与功率谱估计技术结合。

发明内容

本发明是针对目前谱图计算提出的新的计算方法，它结合了现代功率谱估计方法并引入分数域最优阶思想提高局部功率谱的分辨率，以提高整体信号功率谱图的时频特性。

为了解决上述技术问题，达到上述目的，本发明采用如下技术方案：

本发明提供了一种地震信号的分数域局部功率谱计算方法，其特征在于包括以下步骤：

步骤(1)、对输入信号做p阶分数傅里叶变换，求取分数域信号的时宽和频宽及时间带宽积；

步骤(2)、对变换阶次p遍历，即求得各个阶次下分数域信号的时宽和频宽 及时间带宽积；

步骤(3)、对各个阶次下的时间带宽积进行遍历搜索，找到时间带宽积最小值以及对应的阶次、对应的时宽、频宽；

步骤(4)、根据求得的时宽和频宽构造高斯窗函数，对高斯窗函数进行相应阶次的分数阶傅里叶变换求得最终的最优窗函数，该窗函数一般是一个复函数；

步骤(5)、将信号左端和右端进行周期性延拓，用最优窗函数从信号第一个点为中心开始截取固定长度，即信号本身的长度，截取的信号称为原信号的局部，其功率谱反映的是原信号的局部功率谱特征，对获得局部应用功率谱计算方法，如间接法、AR模型法、Burg谱估计等，即可获得有别于传统谱图的局部功率谱时频分布。

上述技术方案中，对输入信号x(n)做p阶分数傅里叶变换如下：

$x_p\left(u\right)=\sqrt{\frac{(1-j\cot \mathrm{\α})}{2\mathrm{\π}}}{e}^{j\frac{u^2}{2}\cot \mathrm{\α}}\∫x\left(t\right)e^{j\frac{t^2}{2}\cot \mathrm{\α}}{e}^{-jtu\csc \mathrm{\α}}dt,\mathrm{\α}=\frac{p\mathrm{\π}}{2}---\left(1\right)$

其中p为变换阶次，α为变换角度，u为分数域变量，e为常数，j为虚数单位，t为时间变量，x_p(u)为p阶分数域的变换系数，当初始值p＝0时，x_p(u)即为x(n)，计算x_p(u)的时宽Tx_p，对x_p(u)再做傅里叶变换得X_p+1(v)，可计算出频宽Bx_p：

${\mathrm{Tx}}_p=\frac{\[\∫{(u-{\mathrm{\η}}_u)}^2{\left|x_p\left(u\right)\right|}^{2}du{\]}^{1/2}}{\left|\right|x_p\left(u\right)\left|\right|},{\mathrm{Bx}}_p=\frac{\[\∫{(v-{\mathrm{\η}}_v)}^2{\left|X_{p+1}\left(v\right)\right|}^{2}dv{\]}^{1/2}}{\left|\right|X_{p+1}\left(v\right)\left|\right|}---\left(2\right)$

${\mathrm{\η}}_u=\frac{\∫u{\left|x_p\left(u\right)\right|}^{2}dt}{\left|\right|x_p\left(u\right)|{|}^2},{\mathrm{\η}}_v=\frac{\∫v{\left|X_{p+1}\left(v\right)\right|}^{2}dv}{\left|\right|X_{p+1}\left(v\right)|{|}^2}---\left(3\right)$

其中||·||²为信号二范数的平方即信号能量，η_u为信号p阶分数域的时间中心，η_v为信号p阶分数域的频率中心，即p+1阶分数域的时间中心，求取阶次p下的时间带宽积:

TBP_p＝Tx_p·Bx_p    (4) 

上述技术方案中，鉴于分数阶傅里叶变换关于p＝2的对称性，对p在0～2之间遍历，取步长为0.01，即p＝p+0.01，按步骤(1)计算新阶次下的x_p(u)、时宽Tx_p、X_p+1(v)、频宽Bx_p、时间带宽积TBP_P，遍历结束可获得一个时宽向量Tx， 频宽向量Bx，时间带宽积向量TBP。

上述技术方案中，对时间带宽积进行遍历搜索，找到时间带宽积最小值以及对应的阶次及对应的时宽、频宽，

$p_0\arg \underset{0\&le;p<2}{min}TBP\left\{x_p\left(u\right)\right\}---\left(5\right)$

即阶次p₀满足使信号在该阶次下的时间带宽积是最小的。

上述技术方案中，由时宽Tx₀，频宽Bx₀构造分数阶高斯窗函数gwin(u)(如式6)，再对gwin(u)做p₀阶分数傅里叶变换，可得对应的时域局部功率谱最优窗函数g(n)，n为时间序列：

$gwin\left(u\right)=e^{-\frac{{\mathrm{\&pi;u}}^2{\mathrm{Bx}}_0}{{\mathrm{Tx}}_0}}---\left(6\right)$

$g\left(n\right)=\sqrt{\frac{(1-j{\mathrm{cot\&alpha;}}_0)}{2\mathrm{\&pi;}}}{e}^{j\frac{n^2}{2}{\mathrm{cot\&alpha;}}_0}\&Integral;gwin\left(u\right)e^{j\frac{u^2}{2}{\mathrm{cot\&alpha;}}_0}{e}^{-jnu{\mathrm{csc\&alpha;}}_0}du,{\mathrm{\&alpha;}}_0=\frac{p_0\mathrm{\&pi;}}{2}---\left(7\right)$

上述技术方案中，用窗函数g(n)截取输入信号x(n)，获得信号在m时刻的局部y(n,m)＝x(n)·g(n-m)，m为窗的中心(其初始值为1)，y(n,m)是一个局部观测信号，可以结合功率谱估计的多种方法对y(n,m)求取功率谱，当窗函数的变量m逐渐滑动时即获得了信号x(n)的局部功率谱在不同时间下的时频分布。

上述技术方案中，在m时刻的局部功率谱可用直接法求取如下：

$STPSD(m,k)={\left|\underset{n=1}{\overset{N}{\&Sigma;}}y(n,m)e^{-j\frac{2\mathrm{\&pi;}k}{N}n}\right|}^2---\left(8\right)$

N为信号y(n,m)的长度，k为对应的频率点，STPSD(m,k)表示信号在m时刻在频率分量k处的功率谱密度。

上述技术方案中，在m时刻的局部功率谱可用间接法求取如下：

①计算y(n,m)的自相关函数r(l,m):

$r(l,m)=\underset{n=1}{\overset{N}{\&Sigma;}}y(n,m)y^{*}(n-l,m)---\left(9\right)$

其中l为时延；

②计算局部功率谱：

$STPSD(m,k)=\underset{l=-N}{\overset{N}{\&Sigma;}}r(l,m)e^{-j\frac{2\mathrm{\&pi;}k}{N}n}---\left(10\right)$

上述技术方案中，在m时刻的局部功率谱可用AR模型法求取如下：

①此方法是将局部观测信号y(n,m)看作是一方差为的高斯白噪声信号通过一线性时不变最小相位系统产生的，若用M阶自回归模型来描述系统，系统表述为：

$y(n,m)=-\underset{k=1}{\overset{M}{\&Sigma;}}{a}_k^{m}y(n-k,m)+v(n,m)---\left(11\right)$

其中为m时刻系统待求的第k个模型系数，v(n,m)为m时刻观测信号对应的白噪声信号；

②将式(11)两边同时乘以y^*(n-l,m)(y平移l个单位，并取共轭)并取期望可得y的自相关函数构成的方程组式(12)如下：

$r_y(l,m)=\left\{\begin{array}{cc}-\underset{k=1}{\overset{M}{\&Sigma;}}{a}_k^m{r}_y(l-k,m)& l>0\\ -\underset{k=1}{\overset{M}{\&Sigma;}}{a}_k^m{r}_y(-k,m)& l=0\end{array}\right.---\left(12\right)$

式(12)写成矩阵形式即转化为Yule-Walker方程如式(13)：

其中白噪声方差可由下式(14)求得：

${\mathrm{\&sigma;}}_m^2=r_y(0,m)+a_1^m{r}_y(-1,m)+...a_M^m{r}_y(-M+1,m)---\left(14\right)$

③求解方程获得系数则可求得局部功率谱式(15)如下：

$STPSD(m,k)=\frac{{\mathrm{\&sigma;}}_m^2}{{\left|-\underset{l=1}{\overset{M}{\&Sigma;}}{a}_l^m{e}^{-j\frac{2\mathrm{\&pi;}k}{N}l}+1\right|}^2}---\left(15\right)$

与现有技术相比，本发明具有如下优势：

功率谱估计中现代谱估计较经典谱估计有明显较高的分辨率，本发明将其与分数域时频分析技术相结合，通过引入分数域最优阶思想，获得较传统窗函数性能更优秀的复窗函数，在截断信号时能够提高局部功率谱估计效果，从而最终提高地震信号功率谱图的时频特性，为后续地震信号分析提供便利。

附图说明：

图1为地震信号分数域局部功率谱计算方法原理结构流程图。

图2为获取的单道地震信号。

图3为输入信号在不同阶次下的时间带宽积图。3-a为时间带宽积图，3-b为时宽图，3-c为频宽图。

图4为根据最小时间带宽积求取的最优窗。4-a为窗函数实部，4-b为窗函数虚部，4-c为窗函数幅值。

图5为信号局部功率谱图。5-a为传统谱图，5-b为直接法求取的基于分数傅里叶变换的局部功率谱图，5-c为间接法求取的基于分数傅里叶变换的局部功率谱图，5-d为AR模型法求取的基于分数傅里叶变换的局部功率谱图。

具体实施方式

下面将结合附图对如何求取单道地震信号的分数域局部功率谱时频分布作进一步的描述：

1)从地震资料中读取一道地震信号x(n)(如图2所示)，对地震信号做p阶分数傅里叶变换如式(1)所示：

$x_p\left(u\right)={\mathrm{FrFT}}_p\&lsqb;x\left(t\right)\&rsqb;=\sqrt{\frac{(1-j\cot \mathrm{\&alpha;})}{2\mathrm{\&pi;}}}{e}^{j\frac{u^2}{2}\cot \mathrm{\&alpha;}}\&Integral;x\left(t\right)e^{j\frac{t^2}{2}\cot \mathrm{\&alpha;}}{e}^{-jtu\csc \mathrm{\&alpha;}}dt,\mathrm{\&alpha;}=\frac{p\mathrm{\&pi;}}{2}---\left(1\right)$

初始值p＝0时，x_p(u)即为x(n)。计算x_p(u)的时宽Tx_p，对x_p(u)再做傅里叶变换得X_p+1(v)，可计算出频宽Bx_p，如式(2)(3)所示：

${\mathrm{Tx}}_p=\frac{\&lsqb;\&Integral;{(u-{\mathrm{\&eta;}}_u)}^2{\left|x_p\left(u\right)\right|}^{2}du{\&rsqb;}^{1/2}}{\left|\right|x_p\left(u\right)\left|\right|},{\mathrm{Bx}}_p=\frac{\&lsqb;\&Integral;{(v-{\mathrm{\&eta;}}_v)}^2{\left|X_{p+1}\left(v\right)\right|}^{2}dv{\&rsqb;}^{1/2}}{\left|\right|X_{p+1}\left(v\right)\left|\right|}---\left(2\right)$

${\mathrm{\&eta;}}_u=\frac{\&Integral;u{\left|x_p\left(u\right)\right|}^{2}dt}{\left|\right|x_p\left(u\right)|{|}^2},{\mathrm{\&eta;}}_v=\frac{\&Integral;v{\left|X_{p+1}\left(v\right)\right|}^{2}dv}{\left|\right|X_{p+1}\left(v\right)|{|}^2}---\left(3\right)$

求取阶次p下的时间带宽积:

TBP_p＝Tx_p·Bx_p    (4) 

2)鉴于分数阶傅里叶变换关于p＝2的对称性，对p在0～2之间遍历，一般可取步长为0.01，即p＝p+0.01，按步骤1)计算新阶次下的x_p(u)、时宽Tx_p、X_p+1(v)、频宽Bx_p、时间带宽积TBP_P，遍历结束可获得一个时宽向量Tx(如图3-b)，频宽向量Bx(如图3-c)，时间带宽积向量TBP(如图3-a)。

3)从图3-a遍历搜索可知，最小TBP对应阶次p₀＝1，以及对应时宽Tx₀＝23.51， 频宽Bx₀＝37.41。

$p_0\arg \underset{0\&le;p<2}{min}TBP\left\{x_p\left(u\right)\right\}---\left(5\right)$

4)由时宽Tx₀，频宽Bx₀构造分数阶高斯窗函数gwin(u)，再对gwin(u)做p₀阶分数傅里叶变换，可得对应的时域局部功率谱最优窗函数(如图4)：

$gwin\left(u\right)=e^{-\frac{{\mathrm{\&pi;u}}^2{\mathrm{Bx}}_0}{{\mathrm{Tx}}_0}}---\left(6\right)$

g(n)＝FrFT_p0[gwin(u)]    (7)

5)用窗函数g(n)截取输入信号x(n)，获得信号的局部y(n,m)＝x(n)*g(n-m)。当m固定时，y(n,m)是一个局部观测信号，可以结合功率谱估计的多种方法对y(n,m)求取功率谱，当窗函数的变量m逐渐滑动时即获得了信号x(n)的局部功率谱在不同时间下的时频分布。

方法1：用直接法求取局部功率谱，结果如图5-b：

$STPSD(m,k)={\left|\underset{n=1}{\overset{N}{\&Sigma;}}y(n,m)e^{-j\frac{2\mathrm{\&pi;}k}{N}n}\right|}^2---\left(8\right)$

方法2：用间接法求取局部功率谱，结果如图5-c：

①计算y(n,m)的自相关函数r(l,m):

$r(l,m)=\underset{n=1}{\overset{N}{\&Sigma;}}y(n,m)y^{*}(n-l,m)---\left(9\right)$

②计算局部功率谱：

$STPSD(m,k)=\underset{l=-N}{\overset{N}{\&Sigma;}}r(l,m)e^{-j\frac{2\mathrm{\&pi;}k}{N}n}---\left(10\right)$

方法3：用AR模型法求取局部功率谱：

①此方法是将观测信号y(n,m)看作是一方差为的高斯白噪声信号通过一线性时不变最小相位系统产生的。若用M阶自回归模型来描述系统，系统表述为：

$y(n,m)=-\underset{k=1}{\overset{M}{\&Sigma;}}{a}_k^{m}y(n-k,m)+v(n,m)---\left(11\right)$

②将式(11)两边同时乘以y^*(n-l)并取期望可得式(12)如下：

$r_y(l,m)=\left\{\begin{array}{cc}-\underset{k=1}{\overset{M}{\&Sigma;}}{a}_k^m{r}_y(l-k,m)& l>0\\ -\underset{k=1}{\overset{M}{\&Sigma;}}{a}_k^m{r}_y(-k,m)& l=0\end{array}\right.---\left(12\right)$

式(12)写成矩阵形式即转化为Yule-Walker方程：

其中白噪声方差可由下式求得：

$r_y(0,m)+a_1^m{r}_y(-1,m)+...a_M^m{r}_y(-M+1,m)={\mathrm{\&sigma;}}_m^2---\left(14\right)$

③求解方程获得系数则可求得局部功率谱如下，结果如图5-d：

$STPSD(m,k)=\frac{{\mathrm{\&sigma;}}_m^2}{{\left|-\underset{l=1}{\overset{M}{\&Sigma;}}{a}_l^m{e}^{-j\frac{2\mathrm{\&pi;}k}{N}l}+1\right|}^2}---\left(15\right)$ 。

Claims (9)

1.一种地震信号的分数域局部功率谱计算方法，其特征在于包括以下步骤：

步骤（1）、对输入信号做p阶分数傅里叶变换，求取分数域信号的时宽和频宽及时间带宽积；

步骤（2）、对变换阶次p遍历，即求得各个阶次下分数域信号的时宽和频宽及时间带宽积；

步骤（3）、对各个阶次下的时间带宽积进行遍历搜索，找到时间带宽积最小值以及对应的阶次、对应的时宽、频宽；

步骤（4）、根据求得的时宽和频宽构造高斯窗函数，对高斯窗函数进行相应阶次的分数阶傅里叶变换求得最终的最优窗函数，该窗函数一般是一个复函数；

步骤（5）、将信号左端和右端进行周期性延拓，用最优窗函数从信号第一个点为中心开始截取固定长度，即信号本身的长度，截取的信号称为原信号的局部，其功率谱反映的是原信号的局部功率谱特征，对获得局部应用功率谱计算方法，即可获得有别于传统谱图的局部功率谱时频分布。

2.根据权利要求1所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，对输入信号x(n)做p阶分数傅里叶变换如下：

其中p为变换阶次，α为变换角度，u为分数域变量，e为常数，j为虚数单位，t为时间变量，x_p(u)为p阶分数域的变换系数，当初始值p=0时，x_p(u)即为x(n)，计算x_p(u)的时宽Tx_p，对x_p(u)再做傅里叶变换得X_p+1(v)，可计算出频宽Bx_p：

其中||·||²为信号二范数的平方即信号能量，η_u为信号p阶分数域的时间中心，η_v为信号p阶分数域的频率中心，即p+1阶分数域的时间中心，求取阶次p下的时间带宽积:

TBP_p＝Tx_p·Bx_p       (4)。

3.根据权利要求2所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，鉴于分数阶傅里叶变换关于p=2的对称性，对p在0～2之间遍历，取步长为0.01，即p=p+0.01，按步骤（1）计算新阶次下的x_p(u)、时宽Tx_p、X_p+1(v)、频宽Bx_p、时间带宽积TBP_P，遍历结束可获得一个时宽向量Tx，频宽向量Bx，时间带宽积向量TBP。

4.根据权利要求2所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，对时间带宽积进行遍历搜索，找到时间带宽积最小值以及对应的阶次及对应的时宽、频宽，

即阶次p₀满足使信号在该阶次下的时间带宽积是最小的。

5.根据权利要求2所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，由时宽Tx₀，频宽Bx₀构造分数阶高斯窗函数gwin(u)（如式6），再对gwin(u)做p₀阶分数傅里叶变换，可得对应的时域局部功率谱最优窗函数g(n)，n为时间序列：

。

6.根据权利要求2所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，用窗函数g(n)截取输入信号x(n)，获得信号在m时刻的局部y(n,m)=x(n)·g(n-m)，m为窗的中心,其初始值为1，y(n,m)是一个局部观测信号，可以结合功率谱估计的多种方法对y(n,m)求取功率谱，当窗函数的变量m逐渐滑动时即获得了信号x(n)的局部功率谱在不同时间下的时频分布。

7.根据权利要求6所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，在m时刻的局部功率谱可用直接法求取如下：

N为信号y(n,m)的长度，k为对应的频率点，STPSD(m,k)表示信号在m时刻 在频率分量k处的功率谱密度。

8.根据权利要求6所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，在m时刻的局部功率谱可用间接法求取如下：

①计算y(n,m)的自相关函数r(l,m):

其中l为时延；

②计算局部功率谱：

。

9.根据权利要求6所述的一种地震信号的分数域局部功率谱计算方法，其特征在于，在m时刻的局部功率谱可用AR模型法求取如下：

①此方法是将局部观测信号y(n,m)看作是一方差为的高斯白噪声信号通过一线性时不变最小相位系统产生的，若用M阶自回归模型来描述系统，系统表述为：

其中为m时刻系统待求的第k个模型系数，v(n,m)为m时刻观测信号对应的白噪声信号；

②将式(11)两边同时乘以y^*(n-l,m)并取期望可得y的自相关函数构成的方程组式(12)如下：

式(12)写成矩阵形式即转化为Yule-Walker方程如式（13）：

其中白噪声方差可由下式（14）求得：

③求解方程获得系数则可求得局部功率谱式（15）如下：

。