挤压短时傅里叶变换的地震勘探信号随机噪声压制方法与流程

本发明属于勘探地球物理技术领域，具体涉及一种利用挤压短时傅里叶变换的地震勘探信号随机噪声压制方法。

背景技术：

随机噪声主要来源于地面微震、仪器以及由于介质不均匀性造成的散射波。它在地震记录上表现为杂乱无章的振动，并且频带很宽。同时由于复杂的地表和地形条件，如山地、黄土塬、崎岖海底、沙漠、戈壁和多期叠合构造等因素，对深层地震记录产生各种随机噪声干扰。随机噪声干扰会降低地震记录的信噪比及质量，对后续的地震信号处理及解释环节造成很大的影响。因此为了更加充分和准确的应用采集到的地震数据，需要对勘探地震数据中的随机噪声进行压制以提高地震信号的信噪比，进而保证后续处理解释环节的可靠性。

目前已有的地震勘探信号随机噪声压制方法包含傅里叶域滤波类方法，频率-波数域滤波类方法，时频域变换滤波类方法和其它二维小波变换类滤波方法，针对傅里叶域滤波类方法，宽频带地震信号与随机噪声在频率域相互交叠，并不能将两者完全分离；同时随机噪声在时域的分布并不完全均匀，直接在频率域压制噪声并不完全且对信号损伤较大；针对频率-波数域滤波类方法，宽频带地震信号与随机噪声在频率-波数域相互交叠，相较于傅里叶域滤波方法交叠减少，但也不能完全将两者完全分离；同时随机噪声在时间-空间域的分布并不完全均匀，直接在频率-波数域压制噪声并不完全且对信号损伤较大。针对时频域变换滤波类方法，计算效率较低且压制噪声的效果依赖于时频变换工具，时频变换的分辨率及聚集度对滤波影响很大。针对其它二维小波变换类滤波方法，计算量大，且各种变换的分辨率受不确定性原理的约束都不高，因此信号的系数在变换域并不聚焦。

技术实现要素：

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，对高维地震数据体中的每一道信号进行短时傅里叶变换，对大于某个阈值的短时傅里叶变换系数进行挤压操作得到挤压短时傅里叶变换系数，然后利用相邻道挤压短时傅里叶变换系数的位置和大小为约束，同时利用阈值对当前道挤压短时傅里叶变换的系数进行阈值等处理，将处理完的短时傅里叶变换系数进行重构得到压制随机噪声后的信号，为后续地震数据处理也解释提供高信噪比的基础数据。

本发明采用以下技术方案：

挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，先对高维地震数据体中的每一道信号进行短时傅里叶变换；然后对大于设定阈值1的短时傅里叶变换系数进行挤压操作得到挤压短时傅里叶变换系数；再利用相邻道挤压短时傅里叶变换系数的位置和大小为约束，对当前道挤压短时傅里叶变换的系数进行阈值处理；最后将处理完的短时傅里叶变换系数进行重构得到压制随机噪声后的信号，为后续地震数据处理也解释提供高信噪比的基础数据，重复以上步骤直到二维数据体中所有的nx道地震信号全部被压制随机噪声。

具体的，对高维地震数据体中的每一道信号进行短时傅里叶变换具体为：

s1、采集二维含噪地震数据体，选取第m道信号xm(t)作为当前道；

s2、对选取的第m道地震信号xm(t)做短时傅里叶变换，得到xm(t)的短时傅里叶变换为stft_g(xm；t,ω)。

进一步的，步骤s2中，取实数窗函数为g(t)，对第m道地震信号xm(t)做短时傅里叶变换，得到窗函数为g(t)时xm(t)的短时傅里叶变换为stft_g(xm；t,ω)为：

stft_g(xm；t,ω)＝∫rxm(τ)g(τ-t)e^-jω(τ-t)dτ

其中，τ为临时积分变量。

具体的，对大于设定阈值1的短时傅里叶变换系数进行挤压操作得到挤压短时傅里叶变换系数具体为：

s3、对第m道地震信号xm(t)的短时傅里叶变换stft_g(xm；t,ω)进行挤压操作，得到窗函数为g(t)时xm(t)的挤压短时傅里叶变换stft_g(xm；t,ω)；

s4、利用设定的阈值e2获取第m道信号挤压短时傅里叶变换的掩模函数mask(xm；t，w)。

进一步的，步骤s3中，xm(t)的挤压短时傅里叶变换sstft_g(xm；t,ω)为：

其中，ωt为临时积分变量，δ(ω)为脉冲函数，当ω为0时该函数取值为1，否则为0。

进一步的，步骤s4中，如果第m道信号的挤压短时傅里叶变换sstft_g(xm；t，w)系数模值大于等于阈值e2，对应的掩模函数mask(xm；t，w)设定为1；如果第m道信号的挤压短时傅里叶变换sstft_g(xm；t，w)系数模值小于阈值e2，对应的掩模函数mask(xm；t，w)设定为0。

具体的，利用相邻道挤压短时傅里叶变换系数的位置和大小为约束，对当前道挤压短时傅里叶变换的系数进行阈值处理具体为：

s5、选取第m道左右各q道，利用从第m-q道到第m+q道挤压短时傅里叶变换的掩模函数共同确定第m道挤压短时傅里叶变换的最终掩模函数final_mask(xm；t，w)；

s6、选取第m道左右各p道，利用从第m-p道到第m+p道挤压短时傅里叶变换系数共同确定第m道挤压短时傅里叶变换的平均系数mean_sstft_g(xm；t，w)。

进一步的，步骤s5中，若mask(xm-q；t，w)、mask(xm-q+1；t，w)，……，mask(xm；t，w)，mask(xm+1；t，w)，……，mask(xm+q，t，w)中有一个取值为1，则第m道该点对应的掩模函数final_mask(xm；t，w)设定为1；若mask(xm-q；t，w)、mask(xm-q+1；t，w)，……，mask(xm；t，w)，mask(xm+1；t，w)，……，mask(xm+q，t，w)全部为0，则第m道该点对应的掩模函数final_mask(xm；t，w)设定为0。

进一步的，步骤s6中，第m道挤压短时傅里叶变换的平均系数mean_sstft_g(xm；t，w)为：

其中，p为选取的左右相邻地震道的道数，n为临时求和变量，sstft_g(xm+n；t,ω)为二维数据中第m+n道的挤压短时傅里叶变换。

具体的，利用第m道挤压短时傅里叶变换的最终掩模函数final_mask(xm；t，w)将第m道的平均挤压短时傅里叶变换系数mean_sstft_g(xm；t，w)中的某些系数置零，重构出经过随机噪声压制的第m道地震信号如下：

其中，g(0)为窗函数g(t)在0时刻的取值，ω为积分变量。

与现有技术相比，本发明至少具有以下有益效果：

本发明一种利用挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，对高维地震数据体中的每一道信号进行短时傅里叶变换，对大于某个阈值的短时傅里叶变换系数进行挤压操作得到挤压短时傅里叶变换系数，然后利用相邻道挤压短时傅里叶变换系数的位置和大小为约束，同时利用阈值对当前道挤压短时傅里叶变换的系数进行阈值等处理，将处理完的短时傅里叶变换系数进行重构得到压制随机噪声后的信号，为后续地震数据处理也解释提供高信噪比的基础数据。充分利用具有高时频分辨率的挤压短时傅里叶变换及地震信号的空间相关性，能够更加充分地压制地震随机噪声，且对信号的损伤较少。

进一步的，对选取的第m道地震信号做短时傅里叶变换，将地震信号拓展到更高的维度，有利于提高后续去噪的灵活性。

进一步的，对第m道地震信号的短时傅里叶变换进行挤压操作，得到其挤压短时傅里叶变换sstft_g(xm；t,ω)，有利于提高有效信号在时频域的时频聚集度。

进一步的，利用设定的阈值e2获取第m道信号挤压短时傅里叶变换的掩模函数，有利于确定噪声和信号在时频域的分布范围。

进一步的，利用从第m-q道到第m+q道挤压短时傅里叶变换的掩模函数(从mask(xm-q；t,ω)到mask(xm+q；t,ω))共同确定第m道挤压短时傅里叶变换的最终掩模函数，有利于利用信号的空间相关性进行去噪。

进一步的，利用从第m-p道到第m+p道挤压短时傅里叶变换系数共同确定第m道挤压短时傅里叶变换的平均系数，有利于利用信号的空间相关性进行去噪。

综上所述，本发明能够充分利用挤压短时傅里叶变换的高时频聚集度特性以及有效信号和随机噪声在时频域聚集度不同等特征，同时可利用高维地震信号的空间相关性，具有效率高、可靠性高、简便易行等特点，可为地震勘探信号随机噪声压制提供了一个全新的方法。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明流程图；

图2为一个含有随机噪声的二维地震数据；

图3为第2道含有随机噪声的地震信号；

图4为第2道地震信号的挤压短时傅里叶变换谱；

图5为第1道地震信号的挤压短时傅里叶变换谱；

图6为第3道地震信号的挤压短时傅里叶变换谱；

图7为经过阈值处理的第2道地震信号的挤压短时傅里叶变换系数绝对值；

图8为经过重构得到的经去噪后的第2道地震信号；

图9为经过去噪后的二维地震数据。先对高维地震数据体中的每一道信号进行短时傅里叶变换，然后对大于设定阈值1的短时傅里叶变换系数进行挤压操作得到挤压短时傅里叶变换系数，接着利用相邻道挤压短时傅里叶变换系数的位置和大小为约束，对当前道挤压短时傅里叶变换的系数进行阈值处理，最后将处理完的短时傅里叶变换系数进行重构得到压制随机噪声后的信号，为后续地震数据处理也解释提供高信噪比的基础数据

具体实施方式

本发明提供了一种挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，对高维地震数据体中的每一道信号进行短时傅里叶变换，然后对大于设定阈值1的短时傅里叶变换系数进行挤压操作得到挤压短时傅里叶变换系数，接着利用相邻道挤压短时傅里叶变换系数的位置和大小为约束，对当前道挤压短时傅里叶变换的系数进行阈值处理，最后将处理完的短时傅里叶变换系数进行重构得到压制随机噪声后的信号，为后续地震数据处理提供高信噪比的基础数据。相比较于常规去噪方法，本发明仅在利用具有高时频分辨率的挤压短时傅里叶变换的同时，能够充分地利用相邻道的空间相关性，能够充分压地震信号中的随机噪声。

请参阅图1，本发明一种挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，包括以下步骤：

s1、在采集到的二维含噪地震数据体中选取第m道信号xm(t)；

若二维含噪地震数据体大小为nt*nx，其中，nt表示时间采样点数，nx表示二维含噪地震数据体中的地震道数，同时获取二维含噪地震数据体的采样间隔为δt；选取第m道地震信号xm(t)作为当前道。

s2、对选取的第m道地震信号xm(t)做短时傅里叶变换，得到xm(t)的短时傅里叶变换为stft_g(xm；t,ω)；

取实数窗函数为g(t)，对第m道地震信号xm(t)做短时傅里叶变换，得到窗函数为g(t)时xm(t)的短时傅里叶变换为stft_g(xm；t,ω)为：

stft_g(xm；t,ω)＝∫rxm(τ)g(τ-t)e^-jω(τ-t)dτ(1)

其中，τ为临时积分变量。

s3、对第m道地震信号xm(t)的短时傅里叶变换stft_g(xm；t,ω)进行挤压操作，得到窗函数为g(t)时xm(t)的挤压短时傅里叶变换stft_g(xm；t,ω)；

首先利用窗函数g(t)的导数g’(t)作为窗函数计算xm(t)的另外一个短时傅里叶变换stft_g'(xm；t,ω)：

stft_g'(xm；t,ω)＝∫rxm(τ)g'(τ-t)e^-jω(τ-t)dτ(2)

其中，τ为临时积分变量。

然后在时频域计算每一点的瞬时频率ωl(xm；t，ω)：

其中，im(z)表示取复数z的虚部。接着将模值大于阈值e1的stft_g(xm；t,ω)累加到其对应的瞬时频率ωl(xm；t，ω)位置上，得到xm(t)的挤压短时傅里叶变换sstft_g(xm；t,ω)：

其中，ωt为临时积分变量，δ(ω)为脉冲函数，当ω为0时该函数取值为1，否则为0。

s4、利用设定的阈值e2获取第m道信号挤压短时傅里叶变换的掩模函数mask(xm；t，w)；

第m道信号挤压短时傅里叶变换的掩模函数mask(xm；t，w)通过如下规则获得：

(a)如果第m道信号的挤压短时傅里叶变换sstft_g(xm；t，w)系数模值大于等于阈值e2，则该点对应的掩模函数mask(xm；t，w)设定为1；

(b)如果第m道信号的挤压短时傅里叶变换sstft_g(xm；t，w)系数模值小于阈值e2，则该点对应的掩模函数mask(xm；t，w)设定为0。

s5、选取第m道左右各q道，利用从第m-q道到第m+q道挤压短时傅里叶变换的掩模函数(从mask(xm-q；t，w)到mask(xm+q；t，w))共同确定第m道挤压短时傅里叶变换的最终掩模函数final_mask(xm；t，w)；

第m道挤压短时傅里叶变换的最终掩模函数final_mask(xm；t，w)可通过如下方法获得：

(a)若mask(xm-q；t，w)、mask(xm-q+1；t，w)，……，mask(xm；t，w)，mask(xm+1；t，w)，……，mask(xm+q，t，w)中有一个取值为1，则第m道该点对应的掩模函数final_mask(xm；t，w)设定为1；

(b)若mask(xm-q；t，w)、mask(xm-q+1；t，w)，……，mask(xm；t，w)，mask(xm+1；t，w)，……，mask(xm+q，t，w)全部为0，则第m道该点对应的掩模函数final_mask(xm；t，w)设定为0。

s6、选取第m道左右各p道，利用从第m-p道到第m+p道挤压短时傅里叶变换系数(从sstft_g(xm-p；t，w)到sstft_g(xm+p；t，w))共同确定第m道挤压短时傅里叶变换的平均系数mean_sstft_g(xm；t，w)；

第m道挤压短时傅里叶变换的平均系数mean_sstft_g(xm；t，w)可通过如下方法获得：

s7、利用第m道挤压短时傅里叶变换的最终掩模函数final_mask(xm；t，w)将第m道的平均挤压短时傅里叶变换系数mean_sstft_g(xm；t，w)中的某些系数置零，然后重构出经过随机噪声压制的第m道地震信号；

将final_mask(xm；t，w)与mean_sstft_g(xm；t；w)逐点相乘，实现利用第m道挤压短时傅里叶变换的最终掩模函数final_mask(xm；t，w)将sstft_g(xm；t，w)中的某些系数置零，在此基础之上重构经过噪声压制的第m道信号

s8、重复步骤s1～s7直到二维数据体中所有的nx道地震信号全部被压制随机噪声。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明一种利用挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，以图2所示的二维合成含噪地震勘探信号为例，其受随机噪声的影响很大。图3为第2道合成含噪地震勘探信号，受随机噪声干扰，并不能清晰地看出有六个反射子波。图4、图5及图6分别为第2道合成含噪地震勘探信号、第1道合成含噪地震勘探信号及第1道合成含噪地震勘探信号的挤压加窗傅里叶变换谱，信号在时频域分布清晰，时频聚集度较高。图7为第2道合成含噪地震勘探信号经过时频域压制噪声后的时频域结果。图8为利用第2道合成含噪地震勘探信号去噪后的挤压加窗傅里叶变换系数重构出的去噪信号，可见六个明显的反射子波。图9为去噪的二维地震勘探信号，相比于原始的含噪信号，信号的有效结构更加清晰。

综上所述，本发明一种挤压短时傅里叶变换的地震勘探信号随机噪声压制方法，采用具有高时频分辨率及重构性能的挤压短时傅里叶变换对含噪信号进行时频分解，充分利用信号与随机噪声在时频域的分布不同，同时考虑信号的空间相关性在时频域压制随机噪声，最后利用挤压时频变换的可重构特性重构出含噪信号。相比较于常规时频域去噪方法，本发明仅在利用具有高时频分辨率的挤压短时傅里叶变换的同时，能够充分地利用相邻道的空间相关性，能够充分压地震信号中的随机噪声。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。