一种微地震信号时频域初至检测方法

摘要

本发明提供一种微地震信号时频域初至检测方法，利用ICEEMDAN进行数据处理时，将微地震信号作为初始数据，在分解的每一阶段添加一个特定白噪声，并计算一个唯一残差以得到每个IMF，ICEEMDAN能自适应地将一个复杂信号分解为一系列IMF分量，且IMF分量满足从高频到低频系列分布；将IMF分量作为输入，对噪声主导的模态直接去除，对其他的模态进行DFA去噪，通过间隔硬阈值去除残余噪声；将不同尺度去噪后的结果进行融合重构获得去噪后的地震记录，即为有效信号；通过高精度时频分析检测得到的有效信号的初至信息。本发明与经验模式分解结果相比，模态混叠现象有了较为明显的降低，能够提供原始信号的精确重构，具有更好的收敛性。

说明书

技术领域

本发明涉及勘探地球物理领域，特别是涉及到地震资料处理中的一种基于总体平均经验模态分解ICEEMDAN与去趋势波动分析法(DFA)联合的微地震信号检测方法。

背景技术

微地震事件信号通常能量较弱，加之传播过程中存在能量损失，这就导致地面接收的地震数据有效信号存在能量弱、信噪比低的缺点，因此提高微地震资料的信噪比是微地震数据处理和解释的首要任务。随着提高定位精度和震源机制全矩张量反演需求的增加，对去噪技术的要求也逐步提高。传统压制随机噪音的方法有很多，可分为空间域及变换域的方法，前者主要包括均值滤波、中值滤波及各向异性扩散滤波等，后者主要包括傅里叶变换域滤波方法、基于小波变换、曲波变换等的阈值去噪方法等。地面微地震资料具有强噪声、弱有效信号的特点，常规的去噪方法往往难以获得较好的去噪效果，因此研究专门针对地面微地震资料的去噪方法是十分重要的。针对经验模态分解(EMD)方法不稳定和产生模态混叠现象(Huang，1998)，总体经验模态分解(EEMD)利用高斯白噪声频谱均匀分布的统计特性，向原始信号中加入不同的白噪声，使得信号在不同尺度上具有连续性，但是该方法计算效率不高(Wu，2009)；完备总体经验模态分解方法(CEEMD)通过加入正、负成对的辅助噪声形式，能够有效消除重构信号中的残余辅助噪声，而且计算效率也能够得到提高(Yeh，2009)，但是重构时精度会有所欠缺；Torres(2011，2014)采用改进的完备总体经验模态分解方法(ICEEMDAN)能够精确重构原始信号，有效减少虚假模态和模态中的噪音，同时计算成本也有所降低。

但是，上述方法尚未能够有效去除噪声，达到精确检测微地震信号的目的。

发明内容

本发明的目的是提供一种基于ICEEMDAN与DFA联合的微地震信号时频域初至检测方法，用以有效去除噪声，以弥补现有技术的不足。

由于ICEEMDAN分解的IMF分量中含有噪音，通过消除波动趋势分析方法(DFA)分离出有效信号和噪音，从而达到消除噪音的目的。因此本发明的目的可通过如下技术措施来实现：

一种基于ICEEMDAN与DFA联合的微地震信号时频域初至检测方法，包括以下步骤：

(1)在利用ICEEMDAN进行数据处理时，将微地震信号作为初始数据，在分解的每一阶段添加一个特定白噪声，并计算一个唯一残差以得到每个IMF，ICEEMDAN能自适应地将一个复杂信号分解为一系列本征模态函数(IMF)分量，且IMF分量满足从高频到低频系列分布；

(2)将(1)处理获得的一系列本征模态函数(IMF)分量作为输入，对噪声主导的模态直接去除，对其他的模态进行DFA去噪；

(3)将去噪后的数据重构获得联合去噪后的结果：将步骤(2)中不同尺度去噪后的结果进行融合重构获得去噪后的地震记录，即为有效信号；

(4)通过高精度时频分析检测步骤(3)得到的有效信号的初至信息。

进一步的，所述步骤(1)中，ICEEDAN的步骤如下描述：

101：通过EMD实现加入不同噪音后信号的计算xⁱ＝x+ε₀E₁(ωⁱ)，之后获得一级残差

102：计算第一个IMF分量IMF₁＝x-r₁；

103：EMD实现r₁+ε₁E₂(ωⁱ)，计算二级残

104：对k＝3,…K，计算k阶残差

105：计算第k个IMF分量IMF_k＝r_k-1-r_k；

106：重复步骤104、105直至残差不能被分解；

其中：定义算子E_j(.)为对给定信号通过EMD求得第j个模态；ωⁱ为单位方差的零均值高斯白噪声，i＝1,2….I，xⁱ＝x+ωⁱ为加入不同噪声后的信号；ε_k允许在每个阶段选择信噪比；M(.)表示局部均值算子，EMD中E₁＝x-M(x)。

进一步的，所述步骤(2)中，DFA去噪方法具体为：

对于一个给定的微地震信号的第i个IMF分量信号IMF_i(t)，计算其累积离差y_i(t)，首先滤去了该序列的平均值，然后进行序列重构，对y_i(t)分别进行等长分割，以长度k将长度为n序列分割成m个不重叠的区间，其中m＝[n/k](取整数)；由于序列长度并不总是增量k的整数倍,因此,序列尾端有时会出现小部分的数据信息未能被利用；为了充分利用数据,对序列的颠倒顺序进行同样的操作，共得到个等长度的区间；然后对每个区间，用最小二乘法分别对每个区间所包含的k个数据进行一阶线性拟合；计算每个区间滤去趋势后的均方差，此处将顺序和逆序分别公式进行计算；对所有等长度区间求均值并开方,计算得到DFA波动函数。

进一步的，所述步骤(2)中，所述DFA波动函数通过阈值范围去除IMF分量来降低白噪音，阈值范围定义为α＝φ±0.5。

本发明的优点和有益效果：

本发明中的总体平均经验模态分解与DFA联合的微地震信号时频域初至检测方法，利用了总体平均经验模态分解具有很好的多尺度的特点以及DFA的有效性,可以压制微地震信号中较强的随机噪声。本发明中的总体平均经验模态分解与DFA联合的微地震信号检测方法，在总体平均经验模态分解的基础上，将其与改进的DFA方法相结合，进一步提高了微地震信号的检测效果。本发明利用总体平均经验模态分解与DFA，得到了高信噪比的微地震处理结果，有利于后续的微地震资料定位与机制分析等，该方法为微地震的反演定位等奠定了基础。

本发明与经验模式分解结果相比，模态混叠现象有了较为明显的降低，能够提供原始信号的精确重构，具有更好的收敛性。去趋势波动分析法(DFA)能够检测含有噪声的微地震信号中哪些内蕴模态是噪声，哪些内蕴模态是纯净信号，进而在保留微地震信号物理意义的前提下，去除微地震中所含有的噪声，提高地震信号的信噪比。

附图说明

图1为本发明的具体流程图；

图2为实施例中微地震信号的波形记录图；

图3是图2中微地震信号经过ICEEMADAN高精度时频分析的结果图；

图4是图2中微地震信号经过ICEEMADAN分解出的本征模态图；

图5是实施例中去噪后的微地震记录和经过ICEEMADAN高精度时频分析的结果图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合所附图式，作详细说明如下。

实施例1：选用某地区的微地震数据。

如图1所示，一种基于ICEEMDAN与DFA联合的微地震信号时频域初至检测方法，包括以下步骤：

(1)在利用ICEEMDAN进行数据处理时，将微地震信号作为初始数据，在分解的每一阶段添加一个特定白噪声，并计算一个唯一残差以得到每个IMF，ICEEMDAN能自适应地将一个复杂信号分解为一系列本征模态函数(IMF)分量，且IMF分量满足从高频到低频系列分布；

(2)将(1)处理获得的一系列本征模态函数(IMF)分量作为输入，对噪声主导的模态直接去除，对其他的模态进行DFA去噪；

(3)将去噪后的数据重构获得联合去噪后的结果，将步骤(2)中不同尺度去噪后的结果进行融合重构获得去噪后的地震记录，即为有效信号；

(4)通过高精度时频分析检测步骤(3)得到的有效信号的初至信息。

上述方法具体为：

(1)选定需要进行噪声压制的微地震波形记录，如图2所示。在利用ICEEMDAN进行数据处理时，将微地震信号看做初始数据，对原始信号进行ICEEMDAN分解，分解出的IMF分量满足从高频到低频的系列分布，图3给出了去噪前所有分量的时频分布。定义算子E_j(.)为对给定信号通过EMD求得第j个模态；ωⁱ为单位方差的零均值高斯白噪声，i＝1,2….I，xⁱ＝x+ωⁱ为加入不同噪声后的信号；ε_k允许在每个阶段选择信噪比；M(.)表示局部均值算子，EMD中E₁＝x-M(x)，ICEEDAN的步骤如下描述：

101：通过EMD实现加入不同噪音后信号的计算xⁱ＝x+ε₀E₁(ωⁱ)，之后获得一级残差

102：计算第一个IMF分量IMF₁＝x-r₁；

103：EMD实现r₁+ε₁E₂(ωⁱ)，计算二级残差

104：对k＝3,…K，计算k阶残差

105：计算第k个IMF分量IMF_k＝r_k-1-r_k；

106：重复步骤104、105直至残差不能被分解；

流程进入到步骤(2)；

(2)将步骤(1)处理获得的分解后的微地震记录，如图4所示，作为输入，对噪声主导的模态可以直接去除，对其他的模态进行DFA去噪，消除波动趋势分析方法是识别含噪音IMF分量的一种新方法：

对于一个给定的微地震信号的第i个IMF分量信号IMFi(t)，计算其累积离差yi(t)，首先滤去了该序列的平均值，然后进行序列重构，对yi(t)分别进行等长分割，以长度k将长度为n序列分割成m个不重叠的区间，其中m＝[n/k](取整数)；由于序列长度并不总是增量k的整数倍,因此,序列尾端有时会出现小部分的数据信息未能被利用。为了充分利用数据,对序列的颠倒顺序进行同样的操作，共得到个等长度的区间。然后对每个区间，用最小二乘法分别对每个区间所包含的k个数据进行一阶线性拟合。计算每个区间滤去趋势后的均方差(此处将顺序和逆序分别公式进行计算)。对所有等长度区间求均值并开方,计算得到DFA波动函数；

改变窗口大小，最小为5个采样周期，最大不超过时间序列样本数量的四分之一。画出对应时窗长度的均方根波动的重对数坐标图，直线斜率为标度指数；标度指数可作为粗糙度的指标；值越大，信号越平稳；简而言之，值越小表示信号更迅速的波动。利用DFA斜率实现基于ICEEMD的去噪算法，阈值能够区别出含噪音的IMF分量；总的来说阈值范围为α＝φ±0.25；该方法第一阶段ICCEMDAN的目的是通过阈值范围去除IMF分量来降低白噪音，因此，阈值范围定义为α＝φ±0.5；处理后如果该IMF分量还存在残余噪声，可以采用间隔硬阈值方法进行后处理去除残余噪声；之后流程进入到步骤(3)；

(3)将步骤(2)中不同尺度去噪后的结果进行融合重构获得去噪后的地震记录，即为联合去噪方法得到的结果，如图5所示；

(4)将步骤(3)中将去噪后的微地震记录进行基于ICEEMDAN的高分辨率时频分析，分析信号的时频分布，为检测有效信号的初至信息。

其中，图3是微地震信号经过ICEEMADAN高精度时频分析的结果；由于微地震数据中包含大量的噪声，有效的初至信息难以有效拾取；图4是图2中微地震信号经过ICEEMADAN分解出的本征模态；图5是去噪后的微地震记录和经过ICEEMADAN高精度时频分析的结果，从图5中可以方便的拾取微地震数据的初至，如图5中虚线所示；从图2到图5可以看出，基于ICEEMADAN分解和DFA的联合去噪方法对于理论模型的应用效果好，噪声得到了有效的压制。本发明提供的联合去噪方法可以有效地压制剖面中的随机噪音，突显有效的微地震信号，由基于CEEMADAN分解和DFA的联合去噪方法对实际微地震记录进行噪声压制效果良好，信噪比得到显著提高。该方法为微地震的反演定位等奠定了基础。