时频域谱减法和经验模态分解联合的隧道滤波方法及系统与流程

本公开属于地球物理勘探技术领域，涉及一种时频域谱减法和经验模态分解联合的隧道滤波方法及系统。

背景技术：

本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。

在隧道超前地质预报中，地震波法超前探测由于对不良地质体边界具有较好的分辨效果，且有效探测距离相对较长，得到了较广泛的应用。在地震数据处理中，提高信噪比有助于更为精确地定位不良地质体和准确成像。在隧道环境中，由于隧道有限的空间和施工干扰影响，在采集的地震数据中存在着“频带范围较宽、振幅强度不同、视速度和传播方向难以确定”的来源广泛的噪声，使得从隧道地震数据中提取有效信号受到严重的制约，成像结果与真实地质情况存在较大误差，增加了解释的难度。因此开展针对隧道环境的地震波法超前探测数据去噪研究具有重要理论和实际意义。

据发明人了解，目前隧道地震超前探测数据去噪的研究还相对较少，除带通滤波等常用基本去噪方法之外，主要有利用视速度差异采用τ-p变换进行波场分离，提取隧道前方的有效反射波，提高预报准确性；在隧道前方反射层倾向与隧道轴线交角变化较大范围内(90°～45°)，利用f-k变换有效分离多方向和多波反射事件，提取主要来自掌子面前方的有效反射波剖面。此外，针对隧道地震超前探测中大量平稳噪声影响数据质量，采用时频域谱减法进行地震超前探测数据处理，可以实现隧道环境下平稳噪声的消减和去除。针对隧道超前探测数据中噪声非平稳、非连续的问题，经验模态分解可将非平稳时间域信号分解为一系列按频率降序排列的固有模态函数分量，降低了原始信号的非平稳特性，进而采用平稳信号去噪方法进行数据去噪。

据发明人了解，实现隧道地震超前探测的数据去噪存在以下两个难题：

第一，隧道地震超前探测中噪声种类较多、振幅较强，易导致有效波波形畸变，造成地震数据解释产生偏差。现有方法往往从地表去噪技术出发，直接将地表方法应用与隧道环境中。然而，隧道内环境相对狭小，炮点和检波器数均远少于地表地震探测，一些地表地震探测数据中常用的去噪方法在隧道环境下往往难以直接应用。

第二，单一数据去噪方法存在一定的局限性，采用单一数据去噪方法进行隧道地震超前探测数据去噪不能达到较好的效果。如采用时频域谱减法对数据进行去噪处理，需满足平稳噪声假设条件，但实际隧道超前探测数据中噪声非平稳、非连续，易导致时频域谱减法中噪声谱估计误差较大，影响滤波效果。传统的经验模态分解认为噪声主要集中于高频段的模态函数中，同时也在其余模态分量中存在，在处理中往往直接将最高频模态分量直接去除，这样做虽然提升了信号的信噪比，但同时这种处理方式会造成高频有效信号的缺失，对后续处理成像等处理带来不利影响。

技术实现要素：

本公开为了解决上述问题，提出了一种时频域谱减法和经验模态分解联合的隧道滤波方法及系统，本公开针对隧道超前探测数据中噪声非平稳、非连续的问题，改进了传统的经验模态分解，采用经验模态分解将非平稳时间域信号分解为一系列按频率降序排列的固有模态函数分量，降低了原始信号的非平稳特性，进而采用时频域谱减法进行数据去噪。

根据一些实施例，本公开采用如下技术方案：

一种时频域谱减法和经验模态分解联合的隧道滤波方法，包括以下步骤：

基于经验模态分解方法，将时域隧道含噪地震记录分解为多个固有模态函数分量；

对分解得到的各个固有模态函数分量进行加权；

对加权后的各个时间域imf分量分别转到时频域进行时频域谱减法滤波；

将各个imf分量进行时频域谱减法滤波后得到的小波系数利用同步挤压小波变换的逆变换分别转换为时域信号；

将获得的时域信号直接相加，实现数据重构和滤波。

作为可选择的实施方式，基于经验模态分解方法，将时域隧道含噪地震记录分解为多个固有模态函数分量的具体过程包括：对时域信号中所有极值点进行识别和选取，利用插值分别绘制上包络线和下包络线，根据包络线计算均值作为迭代的目标函数。

作为进一步的，从原始信号中将其均值减掉，得到余量，用其代替原始信号进行迭代，计算极值包络和余量，直至余量满足预先设定的筛选准则。

作为进一步的，在迭代过程中，对余量进行筛选，筛选条件包括：利用对连续两个筛分标准差sd进行限制：

其中，t是时间域数据的总采样点数，ri,k-1(t)和ri,k(t)分别是在计算第i个imf分量时，相邻位置的两个余量，设定当0.2<sd<0.3时结束筛选。

作为进一步的，经过经验模态分解后获得的各个模态分量再分别通过希尔伯特变换求得各个分量的瞬时频率，并将其在同一个时频谱中进行重构，得到完整数据的时频谱，从而实现对原始数据的时频分析。

作为可选择实施方式，对加权后的各个时间域imf分量分别转到时频域进行时频域谱减法滤波的具体过程包括：对第一个imf分量采用频率方向不加权的窗函数滤波，后续的各个imf分量采用在频率方向加权的窗函数滤波。

作为可选择的实施方式，将获得的时域信号直接相加的具体过程包括：将经过时频域谱减法滤波后获得的各个时域数据对应相加获得最终滤波结果。

一种时频域谱减法和经验模态分解联合的隧道滤波系统，包括：

用于将时域隧道含噪地震记录分解为多个固有模态函数分量的模块；

用于对分解得到的各个固有模态函数分量进行加权的模块；

用于对加权后的各个时间域imf分量分别转到时频域进行时频域谱减法滤波的模块；

用于将各个imf分量进行时频域谱减法滤波后得到的小波系数利用同步挤压小波变换的逆变换分别转换为时域信号的模块；

用于将获得的时域信号直接相加，实现数据重构和滤波的模块。

一种计算机可读存储介质，其中存储有多条指令，所述指令适于由终端设备的处理器加载并执行所述的一种时频域谱减法和经验模态分解联合的隧道滤波方法。

一种终端设备，包括处理器和计算机可读存储介质，处理器用于实现各指令；计算机可读存储介质用于存储多条指令，所述指令适于由处理器加载并执行所述的一种时频域谱减法和经验模态分解联合的隧道滤波方法。

与现有技术相比，本公开的有益效果为：

本公开针对隧道地震探测方法中的地震数据，进行了滤波处理。在实际隧道地震探测数据中存在着非平稳、非连续的噪声，易导致时频域谱减法中噪声谱估计误差较大，影响滤波效果。首先采用经验模态分解降低信号的非平稳特征，然后应用时频域谱减法处理，进一步减小了时频域谱减法处理非平稳信号时的噪声谱估计误差，拓展了时频域谱减法的适用性。

传统的经验模态分解认为噪声主要集中于高频段的模态函数中，同时也在其余模态分量中存在，在处理中往往直接将最高频模态分量直接去除，这样做虽然提升了信号的信噪比，但会造成高频有效信号的缺失，使数据分辨率降低；并且将后续imf分量中，比重较高的低频数据影响引入重构数据之中。相比于地表地震探测数据，隧道地震探测数据的主频相对较高，imf1中有效数据成分相对较多，不宜直接舍弃，本公开通过对imf分量进行加权，在保留部分高频信号的同时，削减低频信号对有效信号的影响。

附图说明

构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。

图1是本实施例的方法流程图；

图2(a)-(b)是本实施例的含噪数据去噪地震记录对比图。

具体实施方式：

下面结合附图与实施例对本公开作进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

本实施例公开了一种时频域谱减法和经验模态分解联合的隧道滤波方法，能够实现对隧道环境下非平稳的地震数据的滤波处理，首先采用经验模态分解将含噪地震数据分解为多个固有模态函数分量，并提出了各固有模态函数分量加权重构方法，降低数据的非平稳特性，然后采用时频域谱减法对各个固有模态函数分量进行滤波处理，最后对处理后的固有模态分量进行求和以完成数据重构，对非平稳的地震数据实现了较好的滤波效果。通过将时频域谱减法和经验模态分解进行联合滤波，拓展了时频域谱减法的适用性，提高了滤波效果。

具体过程如图1所示，包括以下步骤：

步骤s1，基于经验模态分解方法将时域隧道含噪地震记录分解为多个固有模态函数分量；

本实施例主要针对隧道地震探测中噪声较强的情况进行模拟，这里隧道断层正演模型所采用的测线布置形式为从隧道掌子面后方5米处向后沿左右边墙各布置一排检波器，每排30个，间隔1m；在检波器测线后方10m处左右边墙各布置一个炮点。

本实施例所采用的隧道断层正演模型尺寸为150m×80m，其中隧道尺寸为50m×10m，网格宽度δx＝δy＝0.5m，模型四周各额外设置宽度为20的吸收边界，在隧道掌子面正前方50m处设置与隧道轴线呈60°倾角的断层，断层厚度为7m，隧道围岩波速为4000m/s，断层内波速为2000m/s，信号主频300hz，采样间隔为dt＝5×10^-4s，采样点数3000，采样时间共0.15s。

在选取位于右边墙的震源并利用左边墙所有检波器接收的共炮点道集地震记录上加噪，构成信噪比为-10db的含噪信号n1，该含噪记录中前20道数据中反射波同相轴较为清晰，后10道数据中反射波同相轴难以分辨，二次反射波的同相轴在图中完全无法分辨，利用该噪声模拟隧道地震探测中噪声较强的情况。

步骤s2，采用所述经验模态分解方法将上述时域隧道含噪地震记录分解为多个固有模态函数分量；

基于经验模态分解方法将时域隧道含噪地震记录分解为多个固有模态函数分量步骤包括：

对时域信号x(t)中所有极值点进行识别和选取，利用插值分别绘制上包络线emax和下包络线emax，并计算均值m1(t)准备作为迭代的目标函数：

从原始信号x(t)中将其均值m1(t)减掉，得到余量h1(t)，用其代替原始信号x(t)进行迭代，计算极值包络和余量，直至余量满足预先设定的筛选准则。

在整个迭代过程中，对余量hi(t)进行筛选的准则十分重要。目前常见的筛选条件主要为：利用对连续两个筛分标准差sd进行限制：

其中，t是时间域数据的总采样点数，ri,k-1(t)和ri,k(t)分别是在计算第i个imf分量时，相邻位置的两个余量，通常情况下，设定当0.2<sd<0.3时结束筛选。

经过经验模态分解后获得的各个模态分量再分别通过希尔伯特变换求得各个分量的瞬时频率，并将其在同一个时频谱中进行重构，最终获得完整数据的时频谱，从而实现对原始数据的时频分析。

步骤s3，对经验模态分解得到的各个固有模态函数分量，给定相应的加权系数进行加权；

步骤s4，对加权后的各个时间域imf分量分别转到时频域进行时频域谱减法滤波；

步骤s5，将各个imf分量进行时频域谱减法滤波后得到的小波系数利用同步挤压小波变换的逆变换分别转换为时域信号；

对各个imf分量进行时频域谱减法进行滤波处理时，对第一个imf分量采用频率方向不加权的窗函数滤波，后续的各个imf分量采用在频率方向加权的窗函数滤波。

步骤s6，将获得的时域信号直接相加，完成数据重构和全部滤波过程。经过时频域谱减法滤波后获得的各个时域数据对应相加获得最终滤波结果x(t)。

将联合滤波后的地震记录与利用去噪中常用的预测滤波方法进行滤波处理得到的地震记录进行对比，从图2(a)中可以看出，通过预测滤波处理后的地震记录中，噪声幅度整体较高，图中前5道和后5道地震道中噪声振幅明显高于中间道噪声振幅，对于倾斜界面所产生的绕射波，并没有能够得到有效的区分，而是与反射波同相轴混为一团，二次波同相轴与噪声混合，难以分辨；从图2(b)中可以看出，通过联合滤波处理后的地震记录中，噪声振幅明显小于预测滤波处理后噪声的振幅，反射波同相轴清晰可见，反射波同相轴与绕射波同相轴有一个较好的区分，二次反射波同相轴由于本身振幅较小，且存在噪声干扰，在前5道数据中能够分辨出来。

还提供以下产品实施例：

一种时频域谱减法和经验模态分解联合的隧道滤波系统，包括：

用于将时域隧道含噪地震记录分解为多个固有模态函数分量的模块；

用于对分解得到的各个固有模态函数分量进行加权的模块；

用于对加权后的各个时间域imf分量分别转到时频域进行时频域谱减法滤波的模块；

用于将各个imf分量进行时频域谱减法滤波后得到的小波系数利用同步挤压小波变换的逆变换分别转换为时域信号的模块；

用于将获得的时域信号直接相加，实现数据重构和滤波的模块。

一种计算机可读存储介质，其中存储有多条指令，所述指令适于由终端设备的处理器加载并执行所述的一种时频域谱减法和经验模态分解联合的隧道滤波方法。

一种终端设备，包括处理器和计算机可读存储介质，处理器用于实现各指令；计算机可读存储介质用于存储多条指令，所述指令适于由处理器加载并执行所述的一种时频域谱减法和经验模态分解联合的隧道滤波方法。

本领域内的技术人员应明白，本公开的实施例可提供为方法、系统、或计算机程序产品。因此，本公开可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本公开是参照根据本公开实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

上述虽然结合附图对本公开的具体实施方式进行了描述，但并非对本公开保护范围的限制，所属领域技术人员应该明白，在本公开的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本公开的保护范围以内。