一种消除伪频谱的频谱分析方法

摘要

本发明属于数字信号处理技术领域，尤其涉及地震勘探资料处理中消除伪频谱的分析方法。一种消除伪频谱的频谱分析方法，所述频谱分析方法包括剔除伪频谱步骤，所述伪频谱为数字信号处理中正常频谱以外存在的与正常频谱呈周期关系的频谱；所述剔除伪频谱包括对数字信号进行希尔伯特变换，然后进行傅里叶变换，且将经过希尔伯特变换的结果数据设置为傅里叶变换时输入数据的虚部，经过傅里叶变换后输出频谱，所述伪频谱数字信号被剔除。本发明输出的频谱图上“伪频谱”消失，极大地改善了信号的频谱分析效果。具有突出的效果。

说明书

技术领域

本发明属于数字信号处理技术领域，尤其涉及的地震勘探资料处理中消除 伪频谱的分析方法。

背景技术

频谱分析是现代信号处理与分析中的一个重要手段，它是将时域信号变换 至频域加以分析的方法称为频谱分析，其目的是把复杂的时间历程波形经过傅 里叶变换分解为若干单一的谐波分量来研究，从而获得信号的频率结构以及各 谐波和相位信息。

对信号进行频谱分析可以获得更多有用信息，如求得动态信号中的各个频 率成分和频率分布范围，求出各个频率成分的幅值分布和能量分布，从而得到 主要幅度和能量分布的频率值。然而在对信号进行傅里叶变换后，往往存在“伪 频谱”，类似于数字滤波中的“伪门”，即除了正常频谱外，还存在与正常频谱 呈周期关系的频谱，这对信号的频谱分析非常不利，尤其是当真实频谱与“伪 频谱”有重叠频率时会导致求得的频谱错误。本发明所描述的对常规傅里叶变 换进行优化，将其产生的“伪频谱”进行消除，其剔除和消除方法尚未被公开。

发明内容

本发明介绍了一种消除“伪频谱”的方法，该方法利用优化的傅里叶变换， 将常规傅里叶频谱分析中产生的“伪频谱”进行去除。本发明的目的是为频谱 分析提供了一种消除“伪频谱”的方法。

常规傅里叶变换频谱的特点如下，

众所周知，傅里叶变换具有周期性的特点，如公式（1）

$F\left(m\right)=\mathrm{\Delta t}\underset{n=0}{\overset{N-1}{\Sigma }}f\left(n\right)e^{-i2\mathrm{\pi mn}/N}---\left(1\right)$

N是时间域的抽样点数，也是计算出频谱的频率抽样个数，由连续傅里叶变 换过渡到离散傅里叶变换时，应用了公式（2）

$\mathrm{\Delta t\Delta f}=\frac{1}{N}---\left(2\right)$

公式（2）是完成一对DFT（Discrete Fourier Transform）的条件，否则 就不能进行傅里叶正反变换的对应计算。可以计算出N就是傅里叶变换的频率 抽样点周期，于是，（1）式可以写为：

$F(m+N)=\mathrm{\Delta t}\underset{n=0}{\overset{N-1}{\Sigma }}f\left(n\right)e^{-i2\pi (m+N)n/N}$(3)

$=\mathrm{\Delta t}\underset{n=0}{\overset{N-1}{\Sigma }}f\left(n\right)e^{-i2\mathrm{\pi mn}/N}{e}^{i2\mathrm{\pi n}}$由于e^-i2πn＝1，故（3）式可写为：

$F(m+N)=\mathrm{\Delta t}\underset{n=0}{\overset{N-1}{\Sigma }}f\left(n\right)e^{-i2\mathrm{\pi mn}/N}$(4)

$=F\left(m\right)$

式（4）表示F(m)确实是以N为频率抽样点数的周期，它表示计算F(m)时， 如果给定的f(n)是N个值，那么只要计算N个F(m)就行了，在多计算就重复了。 如图（1）所示，此时N为44。

本发明的技术思路及技术实现

由于傅里叶变换存在频率抽样点数的周期性，这对样点数较少且采样率较 高的数据频谱分析带了干扰，导致数据的真实频谱与它的周期频谱（本专利称 为“伪频谱”）交织重叠在一起，目前尚无对策。

首先对信号进行希尔伯特变换，然后将结果数据放入傅里叶变换时输入数 据的虚部，然后进行傅里叶变换，计算后的频谱图上“伪频谱”消失，极大地 改善了信号的频谱分析效果，尤其是频谱与其周期频谱交织重叠在一起时。具 体见流程图，如图（2）所示。

设实函数道x(t)的希尔伯特变换为h(t)，则希尔伯特变换表达式为：

$h\left(t\right)=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{x\left(\tau \right)}{t-\tau }\mathrm{d\tau }---\left(5\right)$

在常规傅里叶变换的程序设置中将输入时的虚部填零，本发明创新性地将 原数据希尔伯特变换后的数据输入傅里叶变换时的虚部，由于x(t)与h(t)在希尔 伯特域正交，故其频率谱中的“伪频谱”逼近于零。

具体的，本发明的包括，

一种消除伪频谱的频谱分析方法，所述频谱分析方法包括剔除伪频谱步骤， 所述伪频谱为数字信号处理中正常频谱以外存在的与正常频谱呈周期关系的频 谱；所述剔除伪频谱包括对数字信号进行希尔伯特变换，然后进行傅里叶变换， 且将经过希尔伯特变换的结果数据设置为傅里叶变换时输入数据的虚部，经过 傅里叶变换后输出频谱，所述伪频谱数字信号被剔除。

对所述数字信号x(t)进行希尔伯特变换得到h(t)，表达式为：

$h\left(t\right)=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{x\left(\tau \right)}{t-\tau }\mathrm{d\tau }---\left(1\right)$

其中，t为时间，τ为时间延迟，+∞为正无穷，-∞为负无穷。

所述傅里叶变换采用公式：

$F\left(\omega \right)={\int }_{-\infty }^{+\infty }f\left(t\right)e^{\mathrm{i\omega t}}\mathrm{dt}---\left(2\right)$

其中t为时间，ω为圆频率，i为虚数，+∞为正无穷，-∞为负无穷。采用欧拉 公式e^±ix＝cosx±isinx简化（2）式，故（2）式可写为：

$F\left(\omega \right)={\int }_{-\infty }^{+\infty }f\left(t\right)(\cos \left(\mathrm{\omega t}\right)-i\sin \left(\mathrm{\omega t}\right))\mathrm{dt}---\left(3\right)$

（3）式虚部为sin(ωt)，令sin(ωt)＝h(t)。

本发明应用在地震信号频谱分析的方法中，对于含噪的地震数字信号道集进 行剔除为频谱的过程包括：

步骤（1）输入含噪地震数据信号道集；

（2）对所述含噪地震数据信号道集进行希尔伯特变换: 含噪地震数字信号x(t)进行希尔伯特变换得到h(t)，其中， t为时间，τ为时间延迟，+∞为正无穷，-∞为负无穷。

（3）进行傅里叶变换：

完成步骤（2）后，进行傅里叶变换，且设h(t)为傅里叶变换中的虚部；

所述傅里叶变换的公式为：

所述傅里叶变换采用公式：

$F\left(\omega \right)={\int }_{-\infty }^{+\infty }f\left(t\right)e^{\mathrm{i\omega t}}\mathrm{dt}---\left(2\right)$

其中t为时间，ω为圆频率，i为虚数，+∞为正无穷，-∞为负无穷。采用欧拉 公式e^±ix＝cosx±isinx简化（2）式，故（2）式可写为：

$F\left(\omega \right)={\int }_{-\infty }^{+\infty }f\left(t\right)(\cos \left(\mathrm{\omega t}\right)-i\sin \left(\mathrm{\omega t}\right))\mathrm{dt}---\left(3\right)$

（3）式虚部为sin(ωt)，令sin(ωt)＝h(t)。

（4）输出频谱：将经过傅里叶变换后的含噪地震数据信号频谱图输出；

（5）频谱分析：就是傅里叶变换后的数列，然后呈图，是频率-振幅谱， 横轴是频率，纵轴是振幅。

本发明经过对数字信号进行希尔伯特变换和傅里叶变换，且在操作中创造 性的将傅里叶变换中的虚部设置为希尔伯特变换后的结果，由于x(t)与h(t)在希 尔伯特域正交，故其频率谱中的“伪频谱”逼近于零。使计算后的频谱图上“伪 频谱”消失，极大地改善了信号的频谱分析效果。具有突出的效果。

附图说明

图1为常规傅里叶变换频谱图

图2为本发明的主要流程图；

图3为实施例1图，其中图3a是未采用本方法的效果图，图3b是采用本方法 后的消除伪频谱的效果图；

图4为实施例1图，是图3的放大图，其中图4a是未采用本方法的效果图，图 4b是采用本方法后的消除伪频谱的效果图；

图5为实施例2图，其中图5a是未采用本方法的效果图，图5b是采用本方法 后的消除伪频谱的效果图；

将结合具体实施方式对附图进行说明

具体实施方式

图2所示为本发明具体的方法图，本发明应用在地震信号频谱分析的方法中， 对于含噪的地震数字信号道集进行剔除为频谱的过程包括：

步骤（1）输入含噪地震数据信号道集；

（2）对所述含噪地震数据信号道集进行希尔伯特变换:

含噪地震数字信号x(t)进行希尔伯特变换得到h(t)，其中， t为时间，τ为时间延迟，+∞为正无穷，-∞为负无穷。

（3）进行傅里叶变换：

完成步骤（2）后，进行傅里叶变换，且设h(t)为傅里叶变换中的虚部；

所述傅里叶变换的公式为：

所述傅里叶变换采用公式：

$F\left(\omega \right)={\int }_{-\infty }^{+\infty }f\left(t\right)e^{\mathrm{i\omega t}}\mathrm{dt}---\left(2\right)$

其中t为时间，ω为圆频率，i为虚数，+∞为正无穷，-∞为负无穷。采用欧拉 公式e^±ix＝cosx±isinx简化（2）式，故（2）式可写为：

$F\left(\omega \right)={\int }_{-\infty }^{+\infty }f\left(t\right)(\cos \left(\mathrm{\omega t}\right)-i\sin \left(\mathrm{\omega t}\right))\mathrm{dt}---\left(3\right)$

（3）式虚部为sin(ωt)，令sin(ωt)＝h(t)。

（4）输出频谱：将经过傅里叶变换后的含噪地震数据信号频谱图输出；

（5）频谱分析：就是傅里叶变换后的数列，然后呈图，是频率-振幅谱， 横轴是频率，纵轴是振幅。

为了验证本方法的正确性，做如下实验。图3a是某信号道频谱图，含噪， 对其进行常规傅里叶变换后，可以看到在有效时间域的抽样点数内，出现了一 个对称的“伪频谱”。当信号道频率较低且采样点数较多时，有效频谱和其“伪 频谱”一般不会交织重叠在一起。而对频率较高且采样点较少的信号，有效频 谱和其“伪频谱”会交织重叠在一起，尤其某些特殊信号类型。

图3是傅里叶频谱图。a）虚部为零时；b）虚部为希尔伯特变换后的数据 时。图3b是对信号道先进行希尔伯特变换，然后将变换后的输入作为傅里叶变 换时的虚部，得到的频谱如图3b所示，和常规虚部频谱为零的频谱图（图3a） 相比，该频谱图有效压制了“伪频谱”，较好地还原了信号的原频谱。

图4放大后傅里叶频谱图。a）虚部为零时；b）虚部为希尔伯特变换后的 数据时。图4a、图4b是图3a和图3b放大后的频谱图，从图中可以看到常规傅 里叶变换后的频谱和希尔伯特变换后填充虚部变换后的频谱图，在有效频谱区 域几乎完全一致，这说明本专利方法消除“伪频谱”是正确可行的。

图5交叉的傅里叶频谱图。a）虚部为零时；b）虚部为希尔伯特变换后的 数据时。图5a所示的常规频谱图是频率较高且样点数较少的信号道，可以看到 频谱图上信号的真实频谱与“伪频谱”交织重叠在一起。图5b是采用希尔伯特 变换后填充虚部变换后后得到的频谱图，从图中可以看到本方法很好地压制了 “假频谱”，完整地呈现了真实频谱。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言， 在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形， 而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是 优选的，而并不具有限制性的意义。