一种基于时频域谱减法的隧道地震波数据去噪方法及系统与流程
本公开属于地球物理勘探领域，具体涉及一种基于时频域谱减法的隧道地震波数据去噪方法及系统。
背景技术：
本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。
在隧道环境下，由于受到隧道空间限制和施工干扰影响，现场采集到的地震信号中往往包含有较多的干扰波，影响地震记录质量，多种类型的噪声相互混叠，严重影响了隧道地震波法超前预报有效信号的提取，从而导致隧道地震探测成像结果中存在假异常，解译结果出现误差。对于隧道地震波法超前探测而言，高信噪比的地震数据是准确定位成像的前提。如果采用合理的去噪方法去除地震数据中的干扰波，则能够提高地震解释结果的精度。
在隧道地震超前探测中，地震信号的信噪比受隧道环境噪声影响较大，开展针对性的去噪方法研究十分必要。目前专注于隧道地震超前探测数据去噪的研究还相对较少，除带通滤波等常用基本去噪方法之外，主要还有：tst超前探测系统中，用到f-k滤波的随机噪声去噪方法，根据围岩波速信息提取同相轴信息，能够更有效的保留反射波同相轴，并在国内一些工程中进行了效果验证；sap超前探测系统中，同时将带通滤波、τ-p滤波、f-k滤波等方法应用到隧道地震超前探测数据去噪处理中，进一步提升了去噪效果。在各种滤波方法中，时空域滤波方法发展相对较早，方法较为成熟；而目前，变换域滤波也逐渐成为地震数据去噪研究的热点和发展方向。其中谱减法是变换域滤波中效果较好的一种方法。
但是据发明人了解，采用谱减法进行隧道地震波数据去噪存在以下难题：采用谱减法对数据进行去噪处理，需满足两个假设条件：其一，噪声需与有效信号相互独立，并满足统计平稳条件；其二，有效信号不能同时分布在数据的所有时段中，数据中至少有约为2-3个有效信号周期的连续纯噪声。由于基本假设与实际情况之间具有差异，在进行谱减法去噪过程中会产生“音乐噪声”和“相位噪声”，其中“音乐噪声”对处理结果的影响更大。但由于隧道中噪声种类复杂、难以确保在时域上保持平稳，易导致噪声谱估计出现误差，引入额外噪声。
技术实现要素：
本公开为了解决上述问题，提出了一种基于时频域谱减法的隧道地震波数据去噪方法及系统，本公开将谱减法转换到时频域进行处理，通过将数据在时频域进一步区分，提升谱减法对噪声谱误差的修正效果，实现对隧道环境下噪声的消减和去除，提高了地震数据低信噪比噪声去除能力。
根据一些实施例，本公开采用如下技术方案：
一种基于时频域谱减法的隧道地震波数据去噪方法，包括以下步骤：
获取时域地震记录数据的采样频率信息，将数据通过傅里叶变换转换到频率域中，进行频谱分析，确定采集信号主频，将时域数据转换至时频域；
基于采样频率信息，建立滤波窗长与信号主频之间的关系，并求得最优滤波窗长；
利用小波系数矩阵中模值最大值的位置，确定时频谱中对应主频段位置，自适应选取高斯滤波系数的峰值位置；
利用高斯概率分布密度函数和对地震数据频谱的分析，自适应选取高斯滤波系数的方差；
在时频域数据各频段分别选取最后若干时窗的小波系数矩阵求平均值，作为各个频段噪声谱的估计值，利用添加系数的高斯滤波窗分频段对噪声谱估计矩阵进行消减，得到去噪后小波系数矩阵，将其转换到时间域中，完成去噪。
作为可选择的实施方式，利用采样频率信息，建立滤波窗长与信号主频之间的关系，并求得最优滤波窗长的具体过程为：
从子波的到达时、波形、波峰、波谷和相位参数的角度进行分析，以选取对雷克子波影响最小的滤波窗为目的，确定最优滤波窗长。
作为可选择的实施方式，自适应选取高斯滤波系数的峰值位置的具体过程包括：
基于时频域变换后的小波系数矩阵中数值大小与原始信号强度相对应，在有效信号强度大于噪声强度的情况下，通过小波系数矩阵中的最大值选取加权系数峰值点；
进行同步挤压小波变换后，原始时域数据可转换为时频域的小波系数矩阵，通过计算矩阵中各元素值的大小，选出小波系数最大的元素，并确定该元素在矩阵中的位置，位置所对应的点为所求的峰值点。
作为可选择的实施方式，自适应选取高斯滤波系数的方差的具体过程包括：根据实际数据特征总结经验公式确定最优方差，最优方差与小波系数矩阵频段数之间的联系，通过进行正演模拟拟合出最优方差与频段数之间的关系曲线。
作为可选择的实施方式，利用添加系数的高斯滤波窗分频段对噪声谱估计矩阵进行消减的具体过程包括：利用添加系数的高斯滤波窗分频段对噪声谱估计矩阵相减，若相减后矩阵的模增大，则改为对噪声谱估计矩阵的一半进行相减。
作为可选择的实施方式，所述添加系数的高斯滤波窗函数为：在传统二维高斯滤波窗函数的幅值并不随时间或频率的改变而改变的基础上，引入与频率相关的加权系数，以增强信号主频附近数据的权重，得到改进后的地震超前探测滤波窗函数。
作为可选择的实施方式，利用小波变换的逆变换，将小波系数矩阵转换到时间域中。
一种基于时频域谱减法的隧道地震波数据去噪系统，包括：
用于获取时域地震记录数据的采样频率信息，将数据通过傅里叶变换转换到频率域中，进行频谱分析，确定采集信号主频，将时域数据转换至时频域的模块；
用于基于采样频率信息，建立滤波窗长与信号主频之间的关系，并求得最优滤波窗长的模块；
用于利用小波系数矩阵中模值最大值的位置，确定时频谱中对应主频段位置，自适应选取高斯滤波系数的峰值位置的模块；
用于利用高斯概率分布密度函数和对地震数据频谱的分析，自适应选取高斯滤波系数的方差的模块；
用于在时频域数据各频段分别选取最后若干时窗的小波系数矩阵求平均值，作为各个频段噪声谱的估计值的模块；
用于利用添加系数的高斯滤波窗分频段对噪声谱估计矩阵进行消减，得到去噪后小波系数矩阵，将转换到时间域中的模块。
一种计算机可读存储介质，其中存储有多条指令，所述指令适于由终端设备的处理器加载并执行所述的一种基于时频域谱减法的隧道地震波数据去噪方法。
一种终端设备，包括处理器和计算机可读存储介质，处理器用于实现各指令；计算机可读存储介质用于存储多条指令，所述指令适于由处理器加载并执行所述的一种基于时频域谱减法的隧道地震波数据去噪方法。
与现有技术相比，本公开的有益效果为：
针对传统谱减法对全频段数据统一处理，噪声谱估计误差较大的问题，提出了利用时频域阈值滤波方法将谱减法转换到时频域中进行处理的时频域谱减法，优化了传统谱减法在频率方向的区分能力。
针对隧道地震超前探测中大量噪声影响数据质量，易诱发假异常和解译偏差的关键问题，将时频域谱减法引入地震超前探测数据处理中，实现了隧道环境下噪声的消减和去除。
利用隧道地震数据主频集中的特征，改进了滤波窗函数，并通过大量模拟算例总结出改进后窗函数系数与原始数据之间的关系，实现了窗函数系数的自适应选取。
附图说明
构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。
图1为一种基于时频域谱减法的隧道地震波数据去噪方法流程图；
图2(a)、图2(b)为本实施例的含噪地震记录去噪效果对比图。
具体实施方式：
下面结合附图与实施例对本公开作进一步说明。
应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
本实施例公开了一种基于时频域谱减法的隧道地震波数据去噪方法，如图1所示，包括以下步骤：
步骤s1，获取时域地震记录数据的采样频率信息，将数据通过傅里叶变换转换到频率域中，并通过频谱分析确定采集信号主频，通过同步挤压小波变换，将时域数据转换至时频域；
本实施例采用sigmoid模型。sigmoid模型由多个倾斜地层、一个断层、一个不整合面和一组正弦状地层组成，其地震记录包含210道数据、采样总时间长度为2s、信号主频为30hz，并施加信噪比为0db的高斯白噪声，以模拟隧道探测环境下噪声强度较弱的情况。
采用同步挤压小波变换对含噪数据x(k)进行时频域转换，同步挤压小波变换以小波变换为基础，首先对信号x(t)进行小波变换，得到小波系数矩阵wx(a,b)：
然后根据帕瓦尔定理，将小波变换转换成频率域形式：
式中，
步骤s2，利用采样频率等信息，建立滤波窗长与信号主频之间的关系，并求得最优滤波窗长；
最优滤波窗长选取如下：
为使信号在滤波窗范围内近似保持平稳，从子波的到达时、波形、波峰、波谷和相位等各种参数的角度进行分析，以选取对雷克子波影响最小的滤波窗为目的，试验并总结出适合地震信号去噪处理的滤波窗长经验公式：
其中：wl表示所选用的窗长；fs表示采样频率；fd表示有效信号的主频，若计算得出的wl非整数时，则一般取整为最接近的奇数。
步骤s3，利用小波系数矩阵中模值最大值的位置，确定时频谱中对应主频段位置，自适应选取高斯滤波系数的峰值位置；
高斯滤波系数的峰值位置的自适应选取包括：
时频域变换后的小波系数矩阵中数值大小与原始信号强度相对应。在有效信号强度大于噪声强度的情况下，可通过小波系数矩阵中的最大值选取加权系数峰值点μ。
进行同步挤压小波变换后，原始时域数据可转换为时频域的小波系数矩阵amn，共分为m个频率段和n个时间段；通过计算矩阵amn中各元素值的大小，可选出小波系数最大的元素max(aij)，并确定该元素在矩阵中的位置(i,j)。此时，i即为所求的峰值点μ。
步骤s4，利用高斯概率分布密度函数和对地震数据频谱的分析，自适应选取高斯滤波系数的方差；
高斯滤波系数的方差的自适应选取包括：
加权系数的最优方差σ的选取目前尚无成熟的经验公式，难以直接建立最优方差σ与小波系数矩阵之间的理论关系。本公开根据实际数据特征总结经验公式确定最优方差。由于引入的加权系数λω与原始数据频率相关，建立最优方差σ与小波系数矩阵频段数f之间的联系，通过进行大量正演模拟拟合出最优方差σ与频段数f之间的关系曲线并构建经验公式：
σ＝0.063×f-0.169
其中，σ为选取的最优加权系数方差，f为原始数据转换到时频域后频率方向的频段数。
步骤s5，各频段分别选取最后若干个时窗的小波系数矩阵求平均值，作为噪声谱估计矩阵；
在获得时频域信号wx(a,b)之后，利用改进后的二维高斯窗函数进行分窗选取：
式中：m、n分别表示频率方向和时间方向上的步数，wl为窗长，wx(i+wl,j+wl)为每个窗口内小波系数矩阵。
在本实施例，在时频域数据wx(i+wl,j+wl)中各频段均选取时间方向最后五列滤波窗内的小波系数矩阵相加取平均，作为各个频段噪声谱的估计值：
式中：
步骤s6，分别用同频段上各时窗的小波系数矩阵与对应频段的噪声的小波系数矩阵估计值相减，并对相减前后的矩阵模的大小进行判断，若相减后所得的小波系数矩阵的模较相减前更大，则减去噪声谱的一半：
式中
步骤s7，通过小波变换的逆变换，将小波系数矩阵转换到时间域中，完成滤波流程。传统谱减法滤波效果和本公开所述方法滤波效果分别如图2(a)、(b)所示。
还提供以下产品实施例：
一种基于时频域谱减法的隧道地震波数据去噪系统，包括：
用于获取时域地震记录数据的采样频率信息，将数据通过傅里叶变换转换到频率域中，进行频谱分析，确定采集信号主频，将时域数据转换至时频域的模块；
用于基于采样频率信息，建立滤波窗长与信号主频之间的关系，并求得最优滤波窗长的模块；
用于利用小波系数矩阵中模值最大值的位置，确定时频谱中对应主频段位置，自适应选取高斯滤波系数的峰值位置的模块；
用于利用高斯概率分布密度函数和对地震数据频谱的分析，自适应选取高斯滤波系数的方差的模块；
用于在时频域数据各频段分别选取最后若干时窗的小波系数矩阵求平均值，作为各个频段噪声谱的估计值的模块；
用于利用添加系数的高斯滤波窗分频段对噪声谱估计矩阵进行消减，得到去噪后小波系数矩阵，将转换到时间域中的模块。
一种计算机可读存储介质，其中存储有多条指令，所述指令适于由终端设备的处理器加载并执行所述的一种基于时频域谱减法的隧道地震波数据去噪方法。
一种终端设备，包括处理器和计算机可读存储介质，处理器用于实现各指令；计算机可读存储介质用于存储多条指令，所述指令适于由处理器加载并执行所述的一种基于时频域谱减法的隧道地震波数据去噪方法。
本领域内的技术人员应明白，本公开的实施例可提供为方法、系统、或计算机程序产品。因此，本公开可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
本公开是参照根据本公开实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。
上述虽然结合附图对本公开的具体实施方式进行了描述，但并非对本公开保护范围的限制，所属领域技术人员应该明白，在本公开的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本公开的保护范围以内。
- 还没有人留言评论。精彩留言会获得点赞！