一种地震信号的分数域局部功率谱计算方法
【技术领域】
[0001] 本发明属于地震信号处理领域。方法基于最优分数阶傅里叶变换，将现代功率谱 估计与分数域时频分析技术相结合，提出了新的局部功率谱计算方法。通过该方法可以提 高信号功率谱图的时频分辨率，可以在地震信号处理中获得高分辨率的局部功率谱特征， 为后续储层流体预测和识别创造有利条件。 技术背景
[0002] 分数阶傅里叶变换起源于1929年Wiener的研宄，兴起于1980年Namias的研宄， 当时他利用特征值的任意次幂运算明确的提出了阶数为分数的傅里叶变换的概念。20世 纪90年代以后，分数傅里叶变换作为一种新型的信号分析工具受到越来越多的关注，并在 光学领域最先得到广泛应用。1994年Almeida揭示了分数傅里叶变换和传统时频分析工具 的关系，得出了分数傅里叶变换可以解释为信号在时频面上坐标轴绕原点逆时针旋转的重 要结论。至此分数傅里叶变换被赋予明确的物理意义，吸引着越来越多的研宄者参与相关 研宄。分数阶傅里叶变换较传统傅里叶变换的一个显著优点就是能为时频分析提供更大的 选择余地，以获得某些性能上的改进。
[0003] 功率谱描述了随机过程中各频率成分平均功率的大小，在工程上广泛应用于随机 信号处理，是随机信号的重要表征形式之一。它源于1898年Schuster提出的周期图概念。 20世纪50年代后，功率谱估计方法发展迅速，间接法、AR模型法、Burg谱估计法、特征空间 分解法等相继出现。到现在信号功率谱密度估计方面已有许多新的算法，这些新方法连同 演变出来的各种算法不下几十种。这些方法作为一种功率谱估计手段在生物医学、语音、雷 达、地震勘探等领域得到了广泛应用。
[0004] 然而功率谱估计给出的是一种基于信号全局特征的频谱信息，这只对平稳信号产 生较好的分析效果。对于非平稳信号，在不同的时间段包含的频谱信息是不一样的，若对 信号直接做功率谱分析，就无法得到信号频谱随时间变化的特征。1946年Gabor提出了 一种对信号局部加窗后再做傅里叶变换的方法即Gabor变换。这种变换能够很好的反映 出信号局部的频谱信息，于是得到了快速发展。随后短时傅里叶变换、Wigner-Vile分布、 Cohen分布、小波变换、S变换等时频分析技术相继出现。功率谱图与短时傅里叶变换诞 生于同一时期，它可由短时傅里叶变换幅值的平方求得，它是一种不同于短时傅里叶变换、 Wigner-Vile分布、小波变换等，且能够反映信号在时间和频率上的能量分布特性的一种优 秀的时频分析技术。
[0005] 在地震勘探领域，利用地震资料直接进行储层流体识别，已成为地震储层预测技 术的主流发展方向之一。频谱分解技术是其中一种行之有效的技术，它主要利用不同性质 的流体对地震波的吸收、衰减不同这一特点，通过地震信号不同层位下的频率属性来反映 地下传输介质的物理特性，从而为判断该层位的油气特性提供依据。地震信号作为一种非 平稳信号，频谱分解的核心是时频分析技术。现有的时频分析技术因为时频分辨率、交叉项 等问题制约着其在地震信号处理中的应用效果，同时也激发了人们寻求新的时频分析方法 的热情。
[0006] 2002年Lufiye等人将分数域傅里叶变换的思想引入到时频分析中，有效地提高 了信号的时频分析效果。分数域时频分析在传统时频分析的基础上利用时频旋转特性，通 过判定某一旋转角度下信号具有的最小时间带宽积来达到最优的时频分辨率。2013年陈颖 频等将基于分数阶的自适应最优阶Gabor变换应用到地震信号频谱分解中，取得了较好的 效果。2015年田琳等利用反卷积方法求取了分数阶Gabor谱图并应用到地震储层频谱分解 中。该方法仍可以看作是在直接法计算功率谱的概念下做的进一步处理，即现阶段的分数 域时频分析技术尚未与功率谱估计技术结合。
【发明内容】
[0007] 本发明是针对目前谱图计算提出的新的计算方法，它结合了现代功率谱估计方法 并引入分数域最优阶思想提高局部功率谱的分辨率，以提高整体信号功率谱图的时频特 性。
[0008] 为了解决上述技术问题，达到上述目的，本发明采用如下技术方案：
[0009] 本发明提供了一种地震信号的分数域局部功率谱计算方法，其特征在于包括以下 步骤：
[0010] 步骤（1)、对输入信号做P阶分数傅里叶变换，求取分数域信号的时宽和频宽及时 间带宽积；
[0011] 步骤（2)、对变换阶次P遍历，即求得各个阶次下分数域信号的时宽和频宽及时 间带宽积；
[0012] 步骤（3)、对各个阶次下的时间带宽积进行遍历搜索，找到时间带宽积最小值以及 对应的阶次、对应的时宽、频宽；
[0013] 步骤（4)、根据求得的时宽和频宽构造高斯窗函数，对高斯窗函数进行相应阶次的 分数阶傅里叶变换求得最终的最优窗函数，该窗函数一般是一个复函数；
[0014] 步骤（5)、将信号左端和右端进行周期性延拓，用最优窗函数从信号第一个点为中 心开始截取固定长度，即信号本身的长度，截取的信号称为原信号的局部，其功率谱反映的 是原信号的局部功率谱特征，对获得局部应用功率谱计算方法，如间接法、AR模型法、Burg 谱估计等，即可获得有别于传统谱图的局部功率谱时频分布。
[0015] 上述技术方案中，对输入信号x(n)做p阶分数傅里叶变换如下：
[0017] 其中p为变换阶次，a为变换角度，u为分数域变量，e为常数，j为虚数单位，t为 时间变量，xp (u)为p阶分数域的变换系数，当初始值p=0时，xp (u)即为x (n)，计算xp (u)的时宽Txp，对Xp(u)再做傅里叶变换得Xp+1 (v)，可计算出频宽Bxp:
[0020] 其中I I *| |2为信号二范数的平方即信号能量，n u为信号P阶分数域的时间中心， nv为信号p阶分数域的频率中心，即p+i阶分数域的时间中心，求取阶次p下的时间带宽 积：
[0021] TBPp=Txp.Bxp (4)
[0022] 上述技术方案中，鉴于分数阶傅里叶变换关于p = 2的对称性，对p在0~2之间 遍历，取步长为〇. 01，即P = P+0. 01，按步骤（1)计算新阶次下的xp(u)、时宽Txp、Xp+1(v)、 频宽Bx p、时间带宽积TBPP，遍历结束可获得一个时宽向量Tx，频宽向量Bx，时间带宽积向 量 TBP。
[0023] 上述技术方案中，对时间带宽积进行遍历搜索，找到时间带宽积最小值以及对应 的阶次及对应的时宽、频宽，
[0025] 即阶次P(l满足使信号在该阶次下的时间带宽积是最小的。
[0026] 上述技术方案中，由时宽TXd，频宽8乂(|构造分数阶高斯窗函数gwin(u)(如式6)， 再对gwin(u)做阶分数傅里叶变换，可得对应的时域局部功率谱最优窗函数g(n)，n为 时间序列：
[0029] 上述技术方案中，用窗函数g(n)截取输入信号x(n)，获得信号在m时刻的局部 y (n, m) = x(n) ? g(n-m)，m为窗的中心（其初始值为1)，y (n, m)是一个局部观测信号，可 以结合功率谱估计的多种方法对y(n，m)求取功率谱，当窗函数的变量m逐渐滑动时即获得 了信号x(n)的局部功率谱在不同时间下的时频分布。
[0030] 上述技术方案中，在m时刻的局部功率谱可用直接法求取如下：
[0032] N为信号y (n，m)的长度，k为对应的频率点，STPSD (m，k)表示信号在m时刻在频 率分量k处的功率谱密度。
[0033] 上述技术方案中，在m时刻的局部功率谱可用间接法求取如下：
[0034] ①计算y(n，m)的自相关函数r(l，m):
[0036] 其中1为时延；
[0037]②计算局部功率谱：
[0039] 上述技术方案中，在m时刻的局部功率谱可用AR模型法求取如下：
[0040] ①此方法是将局部观测信号y (n，m)看作是一方差为的高斯白噪声信号通过一 线性时不变最小相位系统产生的，若用M阶自回归模型来描述系统，系统表述为：
[0042] 其中 < 为m时刻系统待求的第k个模型系数，v (n，m)为m时刻观测信号对应的白 噪声信号；
[0043] ②将式（11)两边同时乘以/(n_l，m)(y平移1个单位，并取共轭）并取期望可得 y的自相关函数构成的方程组式（12)如下：
[0045] 式（12)写成矩阵形式即转化为Yule-Walker方程如式（13):
[0047]其中白噪声方差可由下式（14)求得：
[0049] ③求解方程获得系数<丨，则可求得局部功率谱式（15)如下：
[0051] 与现有技术相比，本发明具有如下优势：
[0052] 功率谱估计中现代谱估计较经典谱估计有明显较高的分辨率，本发明将其与分数 域时频分析技术相结合，通过引入分数域最优阶思想，获得较传统窗函数性能更优秀的复 窗函数，在截断信号时能够提高局部功率谱估计效果，从而最终提高地震信号功率谱图的