一种基于GS变换滤波和EMD去噪的地震信号检测算法的制作方法

本发明涉及地震信号检测准确率领域，具体是一种基于gs变换(广义s变换)滤波和emd(经验模态分解)去噪的地震信号检测算法。

背景技术：

传统的地震信号检测中，对采集到的信号只进行单一频带的滤波，且滤波频带是事先设定的，未必是信号的优势频带，所以滤波过程比较盲目。地震信号中有些震相自身波动不明显，而且后来的震相往往会出现在在前一个震相的尾波中，在信噪比不佳的情况下，这些震相就会淹没在噪声中，不能被检测到，极大降低了地震信号的检测概率。

因此，有必要在检测前对信号进行去噪；地震噪声虽然来源众多且频带分布不规律，但是地震信号的特点是在即将到来的某一时刻，信号突然增强，且地震信号的不同震相，有各自的优势频带。

时频分析技术主要用来对信号进行优势频带的分析。常用的s变换比一些基于傅里叶变换的时频分析方法在频率分辨率和时窗灵活性上有一定的改进，但在地震信号处理时，s变换的基本小波固定的特性，依然限制了信噪比和分辨率的有效提升。对于传统的变换去噪方法，例如小波去噪，阈值的合理选择对去噪效果影响显著。

综上，在对地震信号进行处理时，采用何种方法分析出地震信号的优势频带，去噪时尽可能避免由于阈值设置不当造成检测效果不佳，是地震信号滤波去噪过程中亟待解决的问题。

技术实现要素：

本发明针对目前地震信号处理过程中，对信号只进行单一频带滤波且滤波频带比较盲目，造成去除噪声效果不好，地震信号检测率不高的问题，提出一种基于gs(广义s变换)变换滤波和emd(经验模态分解)去噪的地震信号检测算法，有效提高地震检测信号的特征值，从而有利于地震信号的检测率的提高。

具体步骤如下：

步骤一、截取采集的某地震波信号x(i)，并将该地震波信号进行广义s变换，得到该条地震记录的时频分析图。

广义s变换是指从窗函数或者基本小波对s变换进行改进；

具体改进如下：

1)对窗函数的改进是指：窗函数中加入调节因子γgs，用f/γgs将原来的f替换掉，通过调整调节因子γgs改变时窗宽度，进而改变时频分辨率。

改造后的窗函数如下：

f是指时频分析中的频率，t是时频分析中的时间

则广义s变换为：

x(t)是要进行gs变换的时间域上的信号，i是复数中的虚数单位

2)对基本小波的改进；

改造后的基本小波为：

上式中，a为基本小波的幅值；α为能量衰减率；β为能量延迟时间；f0为基本小波视频率。

则广义s变换为：

步骤二、从时频分析图中观察得到该信号x(i)到来时刻对应的优势频带；

时频分析图为二维图，纵坐标轴是频率值，图中信号强弱对应颜色深浅，频率值即为需要的优势频带。

步骤三、选择最开始的两个优势频带，将地震信号x(i)经过滤波器完成滤波并叠加，得到多频带滤波后的地震信号x′(i)。

每个优势频带各对应一个滤波器；将地震信号x(i)分别经过两个滤波器完成滤波，得到滤波后的信号为x1(i)和x2(i)，将滤波后的信号叠加得到多频带滤波后的地震信号x′(i)。

x′(i)＝x1(i)+x2(i)

步骤四、将多频带滤波后的地震信号x′(i)进行emd分解，并将各分量叠加得到去噪后的地震信号x″(t)。

具体过程如下：

步骤401、利用三次样条插值函数构造出信号局部极大值点xmax(t)和局部极小值点xmin(t)的包络线，并求取上下包络线的平均值m1(t)。

步骤402、从原信号中减去平均包络，得到一个去掉低频的高频新序列h1(t)。

h1(t)＝x(t)-m1(t)(6)

原信号x(t)的初始值为多频带滤波后的地震信号x′(i)。

步骤403、判断高频新序列h1(t)是否满足固有模态函数的条件，如果是，将高频新序列h1(t)记为c1(t)，作为第一个分解得到的固有模态函数imf(1)。否则，将h1(t)作为原始信号，返回步骤401，直至满足条件；

c1(t)＝h1(t)(7)

步骤404、用原信号减去第一个固有模态函数imf(1)得到一个去高频分量的序列x1(t)。

x1(t)＝x(t)-c1(t)(8)

原信号x(t)的初始值为多频带滤波后的地震信号x′(i)。

步骤405、重新将去高频分量的序列x1(t)看做原信号，返回步骤401，即可得到固有模态函数的各阶分量imf(k)，直至满足给定的终止条件时结束。

终止条件为：

k取值1,2,3,4,5…代表阶数，ck-1(t)是指上一阶的固有模态函数，也就是imf(k-1)

步骤406、分解过程完成后，原信号被分解成一系列imf分量和一个趋势项r(t)。

原信号x(t)的初始值为多频带滤波后的地震信号x′(i)。

分解得到的各个分量imf(j)通过简单相加还原原始信号，各imf分量都是原始信号的固有模态特征。

步骤407、各imf分量按照频率由高到低依次排列，将包含信息全部为噪声的imf(j)置零后，将各分量叠加得到去噪后的地震信号x″(t)。

步骤五、将去噪后的地震信号x″(t)采用sta/lta(短时窗平均/长时窗平均)算法进行检测。

sta/lta利用信号在相对短和相对长的时间内能量平均值的比值变化，来反映地震信号的能量变化。

sta/lta具体算法如下：

在时间轴上第i个采样点之前取两个时窗，时窗的长度分别为nsta(短时窗长度)和nlta(短时窗长度)。当sta/lta(短时窗平均/长时窗平均)的比值大于设定的阈值时，表明该信号是地震信号，信号被拾取。

本发明的优点在于：

一种基于gs变换滤波和emd去噪的地震信号检测算法，在一定程度上能够提升信号特征值，从而可以进一步提升检测概率，以期降低误检率。

附图说明

图1是现有技术中原有的传统地震信号检测方法流程图；

图2是本发明一种基于gs变换滤波和emd去噪的地震信号检测算法流程图；

图3为本发明将多频带滤波后的地震信号进行emd分解，并将各分量叠加得到去噪后的地震信号的流程图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明一种基于gs变换滤波和emd去噪的地震信号检测算法(seismicsignaldetectionalgorithmbasedongstransformfilteringandemdde-noising)，针对原有的传统地震信号检测方法，如图1所示，原始算法对地震信号进行一次滤波后，直接对信号进行检测算法处理。本发明在原始算法的基础上进行两方面改进。第一，通过对地震信号进行时频分析，得到地震信号的优势频带，对地震信号针对性的分别进行两次滤波，然后叠加得到新的地震信号。第二，对滤波后的地震信号进行经验模态分解检测，完成去噪，从而提高地震信号的检测概率。

如图2所示，具体步骤如下：

步骤一、截取采集的某地震波信号x(i)，并将该地震波信号进行广义s变换，得到该条地震记录的时频分析图。

对信号进行时频分析的方法很多，本发明采用基于s变换的广义s变换对地震信号进行时频分析；广义s变换是指从窗函数或者基本小波对s变换进行改进，以使此类变换比s变换具有更强的局部性、时频分辨性和无损可逆性。

具体改进如下：

1)对窗函数的改进是指：窗函数中加入调节因子γgs，用f/γgs将原来的f替换掉，通过调整调节因子γgs改变时窗宽度，进而改变时频分辨率。

改造后的窗函数如下：

f是指时频分析中的频率，t是时频分析中的时间

则广义s变换为：

x(t)是要进行gs变换的时间域上的信号，i是复数中的虚数单位

2)对基本小波的改进；

改造后的基本小波为：

上式中，a为基本小波的幅值；α为能量衰减率；β为能量延迟时间；f0为基本小波视频率。

则基本小波可调的广义s变换为：

由上式看出，改造后的基本小波不受时窗的影响，在信号分析的抗噪性和稳定性方面产生影响。不论上述哪种改进方式，广义s变换在继承和发展了s变换以后，性能上有一定提升，完全可以作为提升信号分辨率和信噪比的有效手段。

步骤二、从时频分析图中观察得到该信号x(i)到来时刻对应的优势频带；

时频分析图为二维图，横坐标轴是频率值，图中信号强弱对应颜色深浅，频率值即为需要的优势频带；通过对信号进行时频分析，可以在时频图中清楚的看到地震信号到来时刻对应的频率范围，此频带为地震信号的优势频带，据此设定该信号的滤波频带。这样的做法比目前地震信号滤波中采用统一滤波频带对信号进行处理的方法更有针对性，滤波效果也更好。

步骤三、选择最开始的两个优势频带，将地震信号x(i)经过滤波器完成滤波并叠加，得到多频带滤波后的地震信号x′(i)。

一般地震信号的优势频带不止一个，选取最开始的两个优势频带，以此来设定滤波器的滤波范围。每个优势频带各对应一个滤波器；将地震信号x(i)分别经过两个滤波器完成滤波，得到滤波后的信号为x1(i)和x2(i)，将滤波后的信号叠加得到多频带滤波后的地震信号x′(i)。

x′(i)＝x1(i)+x2(i)

步骤四、将多频带滤波后的地震信号x′(i)进行emd分解，并将各分量叠加得到去噪后的地震信号x″(t)。

目前对地震信号的处理过程中，只有滤波没有去噪；考虑到地震信号噪声来源众多，只进行滤波去除的噪声也有限，故本发明在滤波的基础上对信号进行去噪。

传统的变换去噪方法，例如小波去噪，阈值的合理选择对去噪效果影响显著。而经验模态分解(emd)的方法，能够避免由于阈值选择不当引起去噪效果不佳的情况。通过经验模态分解去噪，重构后信号的信噪比更高些，且实现更为简单。

emd分解后得到固有模态函数的各阶分量imf(m)，且各阶分量按照频率由高到低依次排列，将包含信息全部为噪声的imf(m)置零后，将各分量叠加得到去噪后的地震信号x″(t)。如图3所示，具体过程如下：

步骤401、利用三次样条插值函数构造出信号局部极大值点xmax(t)和局部极小值点xmin(t)的包络线，并求取上下包络线的平均值m1(t)。

步骤402、从原信号中减去平均包络，得到一个去掉低频的高频新序列h1(t)。

h1(t)＝x(t)-m1(t)(6)

原信号x(t)的初始值为多频带滤波后的地震信号x′(i)。

步骤403、判断高频新序列h1(t)是否满足固有模态函数的条件，如果是，将高频新序列h1(t)记为c1(t)，作为第一个分解得到的固有模态函数imf(1)。否则，将h1(t)作为原始信号，返回步骤401，直至满足条件；

c1(t)＝h1(t)(7)

步骤404、用原始信号减去第一个固有模态函数imf(1)得到一个去高频分量的序列x1(t)。

x1(t)＝x(t)-c1(t)(8)

原信号x(t)的初始值为多频带滤波后的地震信号x′(i)。

步骤405、重新将去高频分量的序列x1(t)看做原信号，返回步骤401，即可得到固有模态函数的各阶分量imf(k)，直至满足给定的终止条件时结束。

终止条件为：

ck-1(t)中，k取值1,2,3,4,5…代表阶数，然后之前提到过c1(t)＝imf(1)，所以，ck-1(t)是指上一阶的固有模态函数，也就是imf(k)

步骤406、分解过程完成后，原信号被分解成一系列imf分量和一个趋势项r(t)。

原信号x(t)的初始值为多频带滤波后的地震信号x′(i)。

分解得到的各个分量imf(j)通过简单相加还原原始信号，各imf分量都是原始信号的固有模态特征。

步骤407、各imf分量按照频率由高到低依次排列，将包含信息全部为噪声的imf(j)置零后，将各分量叠加得到去噪后的地震信号x″(t)。

步骤五、将去噪后的地震信号x″(t)作为处理过的地震信号，采用sta/lta(短时窗平均/长时窗平均)算法进行检测。

能量比法中长短时窗法sta/lta，利用信号在相对短和相对长的时间内能量平均值的比值变化，来反映地震信号的能量变化，因其原理简单且实现方便，是目前地震信号检测中最常用的方法。

特征函数反映原始信号的变化特征，直接关系拾取精度。此处使用的特征函数的形式为绝对值，sta/lta具体算法如下：

去噪后的地震信号x″(t)在时间上连续，t用i代表时间采样点，取值是1,2,3,…；在时间轴上第i个采样点之前取两个时窗，时窗的长度分别为nsta(短时窗长度)和nlta(短时窗长度)。当sta/lta(短时窗平均/长时窗平均)的比值大于设定的阈值时，表明该信号是地震信号，信号被拾取。

本发明与只经过单一滤波后再使用长短时窗能量比法(sta/lta)进行地震信号检测的方法相比，能够在低信噪比的情况下，尽可能消除噪声对地震信号的影响，有效提高地震信号的检测特征值，从而可以提高信号的检测概率。