CN112764108B - 一种基于改进经验小波变换的新型地震资料噪声压制算法 - Google Patents

Abstract

本发明公开了一种基于改进经验小波变换的新型地震资料噪声压制算法,根据S变换谱计算取各频点能量,以能量曲线极大值点频率及ε邻域法确定频谱分割边界，进行改进经验小波变换。对分解后的各IMF计算主频，根据能量曲线确定面波IMF及面波频率峰值，构造带通滤波器对面波IMF进行滤波，实现面波压制。根据有效波频率阈值去除随机噪声IMF，得到单道去噪记录。对二维地震资料中的每道均执行以上操作，得到最终去噪记录。本发明能根据地震信号的频率和能量自适应分解，实现面波和随机噪声的同时压制，并尽可能保护有效波，提高地震资料信噪比。

Description

一种基于改进经验小波变换的新型地震资料噪声压制算法

技术领域

本发明涉及地震勘探资料面波和随机噪声同时压制去噪领域，具体涉及一种基于S变换谱频域能量曲线的改进经验小波变换的新型地震资料噪声压制算法。

背景技术

高质量的地震勘探资料获取为油气田开发提供了可靠基础。但实际采集的地震勘探数据中往往会有各种噪声干扰，其中最普遍存在的是面波和随机噪声。为了提高地震勘探的精度，突出有效波，压制面波和随机噪声干扰成为了地震资料处理工作中极其重要的关键问题。

近年来，为有效压制面波干扰，多种解决方法被提出，如S变换，广义S变换，曲波变换，小波变换等，这些方法都起到了一定效果，但也都有其局限性，容易不同程度地造成有效波的损害。如基于小波变换的压制面波的有效性很大程度上取决于小波函数的选取，S变换和广义S变换将时域地震信号变换到时频域，由于面波的频率难以精确确定，在时频域确定面波区域时往往会导致有效波被损害；曲波变换一是受曲波基提取的影响，且去噪效果与阈值选择密切关联，同样存在损失有效信号的不足。在随机噪声处理方面，f-x域滤波，EMD、小波阈值去噪等方法被引入，但这些方法同样有各自的优势和不足，如f-k域滤波易产生假频，EMD存在模态混叠、伪模态、端点效应等问题，而小波阈值去噪的去噪效果易受阈值、阈值函数选择的影响，难以兼顾增强去噪效果和降低有效信号损伤，难以获得较高信噪比。且这些算法一般针对某一种类型的干扰提出，难以同时处理面波和随机噪声。

2013年，业界推出了经验小波变换(Experiential Wavelet Decomposition，EWT)，通过对信号的傅里叶频谱进行分割，准确提取出信号中固有的模态分量，适合处理非线性非平稳信号。EWT推出后，在地震资料随机噪声去除方面有了一定的应用，但鲜有利用EWT压制面波的应用。如基于EWT将地震信号分解后，再利用时频峰值滤波(time-frequencypeak filtering，TFPF)识别有效信号，实现随机噪声抑制；或利用EWT分解出若干本征模函数(Intrinsic Mode Function，IMF)，再利用聚类算法对含噪分量进行阈值化处理。由于EWT需要根据信号频谱幅值进行频带划分，易受噪声影响，而频谱分割方法决定了分解层数的准确性，进而影响信号分解效果。因此，在分析地震信号时，传统EWT的有效性会受影响，有必要根据地震信号的特点进行频谱分割优化，如在计算地震信号的傅里叶频谱后，基于尺度空间表示提取傅立叶谱的慢变分量，再根据尺度参数和中心频率得到地震信号中包含的频率分量和边界，实现自适应频谱分割，但该方法中尺度参数获取依经验获取，其值的选取会影响频谱分割的准确性。

综上，有必要设计一种能同时压制面波和随机噪声，且最大限度保护有效波的地震资料去噪方法，要能依地震信号特征进行自适应的频谱分割和信号分解，实现面波和随机噪声的准确提取，从而实现二者的同时压制，确保压制噪声的同时不造成对有效波的损失，提高地震资料的信噪比。

发明内容

针对上述不足之处，本发明的目的在于提供一种基于改进经验小波变换的新型地震资料噪声压制算法，能根据地震信号的频率、能量分布自适应进行信号分解，从而准确提取面波和随机噪声分量进行压制，并尽可能减少对有效波的损失，提高地震资料信噪比。

为了实现上述目的，本发明提供如下技术方案：一种基于改进经验小波变换的新型地震资料噪声压制算法，具体步骤如下：

步骤S1、获取单道地震记录s(t)的傅里叶频谱F(w),w∈[0,π]和S变换谱S(f,τ)；对F(w)去除趋势项并进行正则化处理；

步骤S2：根据S变换谱计算频域能量，以能量曲线的N-1个极大值点的频率为频谱划分的初始边界集，根据ε邻域法搜索确定新的边界集w₀'＝0,w_N'＝π；

步骤S3：根据新边界集划分频谱，对各子频带构造Meyer经验小波分解地震信号，自适应得到依频率和能量分解的N个模态分量IMF₁……IMF_N；

具体过程如下：

步骤S301：对子频带[w_n-1',w_n']，确定尺度函数φ_n(w)和小波函数

其中LEN为地震序列长度β(x)＝x⁴(35-84x+70x²-20x³)

步骤S302：分解信号，将信号与小波函数和尺度函数分别做内积，得到细节系数和近似系数/>从而得到原信号的各个模态分量，包括低频分量f₀(t)和高频分量f_k(t)。

步骤S4：根据基于S变换谱的频域能量曲线及面波频率范围提取面波IMF，为保护低频有效波，对面波IMF进行带通滤波，实现精准的面波压制。具体过程如下：

步骤S401：对各IMF进行希尔伯特变换HT，得到其瞬时频率，瞬时振幅，并计算每个IMF的主频；

步骤S402：结合基于S变换谱的频域-能量曲线，找到频率在4-18Hz范围以内，且能量明显强于其他频点的子频带位置，确定面波所在IMF；

步骤S403：为保护面波IMF中有效波，根据面波IMF的频谱对该IMF进行带通滤波，对滤波后的面波IMF剔除，实现面波压制。

步骤S5：根据设定的有效地震勘探波频率阈值(300Hz)及各IMF主频，对阈值范围外的IMF进行剔除，去除随机噪声；对剔除了面波和随机噪声的其余IMF进行重构，得到噪声压制后的地震记录。

步骤S6：对二维地震资料中的每一道均执行步骤S1至步骤S5，得到最终的面波和随机噪声同时压制后的地震记录。

进一步的，所述步骤S2中的基于S变换谱的频域能量曲线定义为：对S变换后的S变换谱的每一个频点f_i，计算其对应的能量谱值LE_i，绘制频率-能量谱曲线。

进一步的，所述步骤S2中ε邻域法搜索确定频谱分割边界方法如下：

根据频域能量曲线的N-1个极大值点对应的频率w_i(i＝1...N-1)作为频谱分割的初始边界集w₀＝0,w_N＝π；

对于第k个边界点w_k(k＝1...N-1)，则其将第k-1和k个子频带分开，设第k-1子频带和第K子频带的长度分别为L_k-1＝w_k-w_k-1,L_k＝w_k+1-w_k；以两子频带长度的一半为搜索半径，即在区间[w_k-L_k-1/2,w_k+L_k/2]内重新寻找该区间内的全局最小值，以该最小值为新的边界点w_k'，按此法依次获取初始边界集{w_n}对应的新边界集w₀'＝0,w_N'＝π

进一步的，所述步骤S4中

1)主频计算如下：

其中，f_m为IMF主频，f为基于HT得到的瞬时频率，s(f)为基于HT得到的瞬时振幅谱。

2)带通滤波器的设计如下：

根据基于S变换谱的能量曲线，根据面波主频范围(4-18Hz以内)和能量确定的面波子频带中对应的峰值频点设计滤波器，该滤波器截止频率为面波IMF频谱峰值的6dB带宽频率。

通过采用上述技术，与现有技术相比，本发明具有如下几个特点：

1.本发明利用基于S变换谱的频域能量曲线进行改进EWT，实现依地震信号频率和能量特征的自适应频谱分割和模态分解，具有更优的频谱分割性能，提高了EWT算法在地震信号处理上的适应性。

2.本发明利用改进EWT分解地震信号，可实现面波和随机噪声的同时压制。

3.本发明在确定面波IMF后，利用带通滤波器对面波IMF滤波，可最大限度地保护低频有效波，提高地震资料信噪比。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲，还可以根据这些附图获得其他的附图。

图1：基于S变换谱频域能量曲线的小波变换的噪声压制算法流程图；

图2：本发明基于Ricker子波的单道合成地震波及含噪地震波；

图3：本发明的基于S变换谱的频域能量曲线；

图4：本发明的频谱分割；

图5：传统经验小波的频谱分割；

图6：本发明的改进EWT所得IMF；

图7：传统EWT所得IMF；

图8:本发明滤波前后的面波分量；

图9：本发明的合成模拟地震剖面；

图10：本发明的去噪效果；

图11：本发明滤波前后的面波IMF记录；

图12：本发明的去噪残差；

图13.本发明与基于广义S变换及二维离散小波变换的面波压制方法的面波压制效果对比；

图14.本发明与基于广义S变换及二维离散小波变换的面波压制方法的面波去噪残差对比；

图15.本发明的随机噪声压制效果(含面波)对比

图16.基于EWT的随机噪声压制残差；

图17.本发明与基于EWT的地震随机噪声压制方法的随机噪声压制效果对比(不含面波)；

图18.本发明与基于EWT的地震随机噪声压制方法的随机噪声压制去噪残差(不含面波)。

具体实施方式

为了使本领域的技术人员更好地理解本发明的技术方案，下面将结合附图对本发明作进一步的详细介绍。

需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括……”或“包含……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的要素。此外，在本文中，“大于”、“小于”、“超过”等理解为不包括本数；“以上”、“以下”、“以内”等理解为包括本数。

实施例一：

如图1-18所示，本实施例中的一种基于改进经验小波变换的地震资料噪声压制算法，其特征在于，包括如下步骤：

S1、对单道地震信号s(t)进行FFT和S变换，分别得到其傅里叶频谱F(w),w∈[0,π]和S变换谱S(f,τ)。对频谱F(w)，分别用多项式和标准差为1.5，窗口大小为10的高斯函数进行趋势项去除及正则化处理，以提高原始地震信号的精度；由于S变换采用带有频率变量的高斯窗函数截取信号，从而实现信号的局部分析，由于其窗函数可随频率的要求进行自适应改变，具有多分辨率的特点，能精确标定地震信号的时频特性。

S2、以基于S变换谱的频域能量曲线的极大值的频率为频谱划分的初始边界，根据ε邻域法搜索确定新的频谱划分边界。

传统EWT依据傅里叶频谱的局部极大值分割频谱，极易受噪声影响，，若直接应用其分解地震信号，其有效性必会受影响，有必要根据各频率成分的能量强弱分解信号。而S变换后，可同时从域和频域分析原地震信号的特性和能量分布，因此利用S变换谱能更好的根据信号的频率和能量进行频谱分割。

具体的频谱划分方法如下：

步骤S201：对步骤一所得的S变换谱中每一个频点fi，计算频域能量谱LE_i，实现从信号时频域到频率-能量域的转换，并根据计算的能量谱绘制频域能量曲线。

其中f_i为S变换后频点i的频率，S(f_i,τ)为该频点对应的S变换谱值。

由于面波与有效波在频率和能量上有区别，为提取地震信号中的面波成分，需要将地震信号按频率和能量分解。由于基于S变换谱的频域能量曲线能清晰描绘信号的不同频率成分的能量强弱，而面波能量强，频率低，可利用能量谱曲线进行自适应频谱分割，从而依频率和能量分解信号，实现面波的提取。，

步骤S202：为进行子频带划分，对步骤201所得的能量曲线进行局部极大值求解，获取其N-1个极大值点对应的频率wi(i＝1...N-1)，以其作为频谱分割的初始边界集w₀＝0,w_N＝π。

步骤203：为更好进行频谱分割，利用ε邻域法重新搜索确定边界。对于第k个边界点w_k(k＝1...N-1)，设其分开的第k-1和k个子频带长度分别为L_k-1＝w_k-w_k-1,L_k＝w_k+1-w_k，以两子频带长度的一半为半径进行搜索，在[w_k-L_k-1/2,w_k+L_k/2]区间内检测该区间的全局最小值，以该最小值为新的边界点w_k'，依此法得到新边界集w₀'＝0,w_N'＝π。

该方法确定的边界，各子频带边界点不重叠，且能较好的反映各频带的频率和能量分布；

S3、根据步骤203所得的边界对傅里叶频谱进行分割，进行改进经验小波变换，将地震信号自适应分解为N个按频率和能量分布的模态分量IMF1……IMFN。

具体过程如下：

步骤S301：根据步骤203所得的边界，对其中的子频带[w_n-1',w_n']，根据Meyer小波构造经验小波，确定经验尺度函数φ_n(w)和小波函数

其中LEN为地震序列长度

β(x)＝x⁴(35-84x+70x²-20x³)

步骤302：根据步骤301的经验小波函数和尺度函数，分解地震信号。将原信号与经验小波函数做内积，得到细节系数将原信号与尺度函数做内积得到近似系数从而得到依地震信号频率和能量分解的各个模态分量，包括低频分量f₀(t)和高频分量f_k(t)。

；

S4、根据频域能量曲线及面波频率范围提取面波IMF；为防止面波IMF中含有面波频率附近的低频有效波，对面波IMF进行带通滤波，剔除滤波后的面波IMF，实现精准的面波压制。

一般面波与有效地震波的频带范围是有差异的，其频率比有效地震波的频率要低，一般为4-18Hz，且能量要明显强于有效波。由于本发明的改进EWT是根据S变换后的频域能量来进行频谱分割的，因此各IMF分量是根据原信号不同频率成分的能量强弱来分解的，从而可以将低频、高能量的面波成分从原信号中分离出来。

具体过程如下：

步骤S401：对步骤302所得的各IMF进行HT，根据基于HT得到的瞬时频率f瞬时振幅谱s(f)计算各IMF的主频f_m：

步骤S402：利用步骤201所得的S变换谱频域能量曲线，找到主频在4-18Hz范围内，且能量明显强于其他频点的子频带位置，根据子频带与分解IMF一一对应的关系，确定对应的面波IMF；

步骤S403：根据步骤201所得能量曲线中面波子频带的频谱峰值设计带通滤波器，该滤波器截止频率为面波IMF频谱功率峰值的6dB带宽频率，即功率下降为峰值25％，幅度下降为峰值50％时的频率为截止频率。

对面波IMF，利用该带通滤波器进行滤波，以精确提取面波，保护低频有效波。

对滤波后的面波IMF进行剔除，实现精准的面波压制而不伤害有效波；

S5、根据设定的有效地震有效波频率阈值，剔除主频在阈值范围外的IMF，实现随机噪声压制，得到单道去噪地震记录。

步骤S501：设定有效地震勘探波频率阈值为300Hz；

由于随机噪声一般以高频率存在于地震资料中，且常规地震勘探中检波器接收的有效地震波频率范围一般在300Hz以内，因此根据有效地震波与随机噪声的频率范围不同，可较好的实现随机噪声衰减。

步骤S502：根据步骤501所得有效波阈值和步骤401所得的各IMF主频，对主频大于300Hz的IMF进行剔除。

步骤S503：对剔除了面波和随机噪声的其余IMF进行重构，得到面波和随机噪声同时压制后的地震记录；

S6、对二维或三维地震资料中的每道均执行步骤1至步骤五，得到最终的二维或三维去噪地震记录；

实施例二：

基于ricker子波合成的单道含面波和随机噪声的地震地震记录为例，说明本发明基于S变换谱频域能量曲线的改进经验小波变换分解地震信号的优势。

利用ricker子波进行单道地震信号模拟记录，时间为1s，采样点1000个。设0.1s处面波主频为10Hz的ricker子波，幅值为5，能量较强，0.5s，0.65s，0.85s处为有效波，频率分别为20Hz，35Hz，50Hz，幅值分别为1，2，3，其合成复合波及加入SNR＝10dB噪声后的含噪复合波如图2所示。

对该地震记录进行S变换后，根据S变换所得的S变换谱，对每个频点计算其能量曲线，如图3所示。从图中可以看出，曲线较好的体现了各频率成分的能量分布，图4为基于该曲线的极大值点及ε邻域法进行子频带划分结果，能更好的体现地震信号中面波和有效波的能量、频率的不同，具有更好的频谱分割性能。图5为传统EWT划分的子频带，不能很好的体现面波，因此不适合进行面波压制。

利用本发明的改进EWT得到的IMF如图6所示，而传统的EWT所得IMF如图7所示.从图中可以看出，本发明改进EWT分解所得IMF能很好的根据信号的频率和能量进行模态分解，图6中IMF主要为0.1s处的面波分量(含有微量的有效波)，IMF2主要为0.5s，0.65s，0.85s处的有效波分量，而IMF3主要为高频随机噪声，因此可以实现面波与随机噪声的有效压制。而传统的EWT分解结果如图7，IMF1和IMF2中均含有面波记录，且有效波与面波叠加在其中，因此无法实现面波提取。可见，本发明的改进EWT能根据地震记录中面波和有效波在频率、能量上的不同，进行自适应模态分解。

为防止面波IMF中的有效波成分被压制，根据S变换的频域-能量曲线中的面波频率峰值，以6dB带宽的带通滤波器进行滤波，保留面波频率范围外的有效波，如图8所示，滤波前面波IMF中含有少量有效波，滤波后，绝大部分有效波被保留，二者与理想纯面波的标准差为0.0202、0.0109；方差为0.1420、0.1042，可见减少了对有效波的损失，提高信噪比。

以合成的含面波和随机噪声的二维地震剖面测试本发明去噪的效果，如图9所示为根据层状速度模型和波场模拟方式生成的模拟地震剖面。面波速度为200m/s,直达波速度为2000m/s，合成地震剖面采样间隔为1ms，共1000个采样点，10道，包含面波、反射波，面波主频为12Hz，反射波频率为40Hz，含噪剖面中添加信噪比为15dB的随机噪声。可以看出，面波与有有效反射波有叠加，导致同相轴模糊。

利用本方法进行面波和随机噪声压制后得到的剖面如图10所示，从图中可以看出，面波得到有效压制，有效地震信号更光滑、更平稳，且有效反射波与面波混叠的同相轴变得清晰，说明了本算法在面波压制方面的性能。

本发明改进EWT分解后的得到的高频IMF记录如图11.a，经滤波后的IMF记录如图11.b，可见滤波后大部分有效波被保留，仅极少部分有用信息被破坏。

最终被去除的噪声残差如图12所示，可见，本发明所得的残差中基本为噪声，仅损失极少量的有效波，且面波压制比较彻底。

为测试本发明的面波压制效果，与基于广义S变换及二维离散小波变换的面波压制方法进行对比，因该算法只适合处理面波，因此，对图9a中的不含随机噪声的地震剖面，用基于广义S变换及二维离散小波变换的面波压制方法去除面波，效果如图13.a所示，从图中可见，面波压制不彻底，在去噪地震剖面中仍残留有少部分面波，从图14所示的去噪残差中可见，基于广义S变换及二维离散小波变换的面波压制方法中损失了部分有效波，导致去噪效果不理想，这是因为该算法在利用广义S变换识别出面波后，利用二维离散小波变换时选择的是高通滤波，这无疑会损失掉一部分低频有效信号，从而影响去噪效果。13.b中为本发明面波压制效果，从图中可以看出，面波去除较彻底，从图14.b中的本发明去噪残差可见，有效波得到充分保护，去噪效果较好。

为对比随机噪声压制效果，对图9.b中的地震剖面用“基于EWT的随机噪声压制算法”去噪，效果如图15.a所示，可以看出，面波未得到有效压制，随机噪声得到压制，从图16所示的基于EWT的去噪残差中可以看出，部分有效波被损害，而从图15.b可见，本发明不仅能压制面波，还能有效压制随机噪声，最大限度保护有效波。

为对比二者性能，采用不含面波，仅含随机噪声的地震剖面进行去噪效果对比，去噪效果如图17所示，去噪残差如图18可见，本发明压制的基本全为噪声，体现了较好的去噪效果，而基于EWT的去噪算法中含有少部分的有效信号。这是因为不同于传统EWT按照傅里叶频谱的极值点分割频谱，本发明在频谱分割时是根据频点的能量来进行的，从而可使面波、有效波、随机噪声根据各自的能量大小自适应分割，从而能更好的实现高频随机噪声的压制。

以上只通过说明的方式描述了本发明的某些示范性实施例，毋庸置疑，对于本领域的普通技术人员，在不偏离本发明的精神和范围的情况下，可以用各种不同的方式对所描述的实施例进行修正。因此，上述附图和描述在本质上是说明性的，不应理解为对本发明权利要求保护范围的限制。

Claims (4)

1.一种基于改进经验小波变换的新型地震资料噪声压制算法，其特征在于，具体步骤如下：

步骤S1、获取单道地震记录s(t)的傅里叶频谱F(w),w∈[0,π]和S变换谱S(f,τ)；对F(w)去除趋势项并进行正则化处理；

步骤S2：根据S变换谱计算频域能量曲线，以能量曲线的N-1个极大值点的频率为频谱划分的初始边界集，根据ε邻域法搜索确定新的边界集

步骤S3：根据新边界集划分频谱，对各子频带构造Meyer经验小波分解地震信号，自适应得到依频率和能量分解的N个模态分量IMF₁……IMF_N；

具体过程如下：

步骤S301：对子频带[w_n-1',w_n']，确定尺度函数φ_n(w)和小波函数

其中LEN为地震序列长度β(x)＝x⁴(35-84x+70x²-20x³)

步骤S302：分解信号，将信号与小波函数和尺度函数分别做内积，得到细节系数和近似系数/>从而得到原信号的各个模态分量，包括低频分量f₀(t)和高频分量f_k(t)；

步骤S4：根据基于S变换谱的频域能量曲线及面波频率范围提取面波IMF，为保护低频有效波，对面波IMF进行带通滤波，实现精准的面波压制；具体过程如下：

步骤S401：对各IMF进行希尔伯特变换HT，得到其瞬时频率，瞬时振幅，并计算每个IMF的主频；

步骤S402：结合基于S变换谱的频域-能量曲线，找到频率在4-18Hz范围以内，且能量明显强于其他频点的子频带位置，确定面波所在IMF；

步骤S403：为保护面波IMF中有效波，根据面波IMF的频谱对该IMF进行带通滤波，对滤波后的面波IMF剔除，实现面波压制；

步骤S5：根据设定的有效地震勘探波频率阈值(300Hz)及各IMF主频，对阈值范围外的IMF进行剔除，去除随机噪声；对剔除了面波和随机噪声的其余IMF进行重构，得到噪声压制后的地震记录；

步骤S6：对二维地震资料中的每一道均执行步骤一至步骤五，得到最终的面波和随机噪声同时压制后的地震记录。

2.如权利要求1所述的一种基于改进经验小波变换的新型地震资料噪声压制算法，其特征在于，所述步骤S2中的基于S变换谱的频域能量曲线定义为：对S变换后的S变换谱的每一个频点f_i，计算其对应的能量谱值LE_i，绘制频率-能量谱曲线

。

3.如权利要求2所述的一种基于改进经验小波变换的新型地震资料噪声压制算法，其特征在于，所述步骤S2中ε邻域法搜索确定频谱分割边界方法如下：

根据频域能量曲线的N-1个极大值点对应的频率w_i(i＝1...N-1)作为频谱分割的初始边界集

对于第k个边界点w_k(k＝1...N-1)，则其将第k-1和k个子频带分开，设第k-1子频带和第K子频带的长度分别为L_k-1＝w_k-w_k-1,L_k＝w_k+1-w_k；以两子频带长度的一半为搜索半径，即在区间[w_k-L_k-1/2,w_k+L_k/2]内重新寻找该区间内的全局最小值，以该最小值为新的边界点w_k'，按此法依次获取初始边界集{w_n}对应的新边界集

4.如权利要求3所述的一种基于改进经验小波变换的新型地震资料噪声压制算法，其特征在于，所述步骤S4中

主频计算如下：

其中，f_m为IMF主频，f为基于HT得到的瞬时频率，s(f)为基于HT得到的瞬时振幅谱；

带通滤波器的设计如下：

根据基于S变换谱的能量曲线，根据面波主频范围和能量确定的面波子频带中对应的峰值频点设计滤波器，该滤波器截止频率为面波IMF频谱峰值的6dB带宽频率。