CN111007569A - 一种集成经验模态分解的样本熵阈值微地震信号降噪方法 - Google Patents

Abstract

本公开公开了一种集成经验模态分解的样本熵阈值微地震信号降噪方法，所述方法包括：步骤S1，获取含有噪声的微地震信号，并对所述微地震信号进行集成经验模态分解，得到一系列不同频率成分的本征模态函数IMF；步骤S2，基于样本熵从所述步骤S1得到的本征模态函数中提取符合微地震信号特征的目标本征模态函数；步骤S3，对于所述目标本征模态函数进行重构，得到去噪后的微地震信号。本公开利用EEMD方法有效解决微地震信号中存在的模态混叠现象，进而能够有效去除微地震信号中的噪声，提高微震数据的信噪比。

Description

一种集成经验模态分解的样本熵阈值微地震信号降噪方法

技术领域

本公开属于信息处理技术领域，尤其涉及一种集成经验模态分解的样 本熵阈值微地震信号降噪方法。

背景技术

地下工程岩爆的发生由于其不仅可能会造成人员伤亡、设备损坏、工 期延误等情况，严重时甚至还可能会诱发地震，因此已经成为地下工程在 施工和运行期间必须解决的问题之一。据美国国家安全与健康研究所 (NIOSH)的Iannacchione等人对Moonee矿微地震数据的整理和研究发现： 90％的顶板塌落事件可以提前1分钟被预测，73％的顶板塌落事件可以提 前54分钟被预测(Iannacchione，2005)。2017年成都理工大学隧道与地 下工程团队利用微地震监测技术成功预测了米仓山隧道岩爆，避免了施工 人员伤亡(顾华宁，2017)，冯夏庭教授团队在锦屏二级水电站深埋隧道 群建立了基于微地震监测技术的岩爆预警方法，在监测期间，共发生275 次岩爆，预报了其中的243次，未预测到的以轻微岩爆为主(丰光亮等， 2016)。

微地震监测技术已经被成功的应用于大坝矿山的安全监测、页岩气开 采中水力压裂监测、深埋地下工程和隧洞稳定性监测、矿山盗采等领域 (Obert et al,1941；Xu etal,2016；Song et al,2010；撒利明等,2012；宋广 东等,2016；沈统等,2017)。微地震监测技术相比于其他监测技术而言，其 优势是可以描述监测区域内微地震事件的时空分布，然而，由于微地震检 波器自身的技术局限性以及外部监测环境等因素的影响，实际采集的微地 震资料不可避免地与噪声混合在一起，如电磁干扰、背景噪声、机械振动、 爆破冲击等。这些噪声通常具有非线性、非平稳性和随机性的特点，甚至 会导致微地震信号产生畸变，掩盖岩石破裂扩展的真实过程，进而影响震 源位置定位、震源机理解释以及其他后续工作。微地震信号与外界随机噪 声的频带存在部分重叠，传统的线性滤波和频谱分析方法的去噪效果并不 理想[何柯等，2016]。集成经验模态分解在非平稳信号处理领域显示出独特的优势，该方法在不需要信号先验知识的前提下，能够自适应地将原信 号分解为一系列本征模态函数(IMF)，理论上，在信号去噪处理领域，只 要找出信号和噪声之间的界限后再进行重构，就可以得到去噪信号。因此 需要寻找合适的方法区分IMF分量中的信号和噪声。

目前针对微地震信号中非平稳随机噪声的去噪方法主要有：小波阈值 法[Beenamol et al,2012；Gaci et al,2014]，经验模态分解法(EMD)[李月等， 2013；陈凯，2009]等等。其中，小波阈值法具有计算量小、计算简单等优 点，但由于目前还没有确定小波基和分解水平的方法，因此限制了小波阈 值法的使用。经验模态分解(EMD)方法是傅立叶变换、小波变换等传统时 频分析方法的重大突破，此方法在非平稳信号处理领域显示出独特的优 势，已被广泛应用于机械故障诊断、图像处理、地震信号处理等领域[Li etal.2015；Linderhed et al.2014；Wang et al.2012]，该方法在不需要信号先验 知识的前提下，能够自适应地将原信号分解为一系列本征模态函数(IMF)， 并根据信号的特性，按频率进行降序排列。理论上，在信号去噪处理领域， 只要找出信号和噪声之间的界限后再进行重构，就可以得到去噪信号。然 而，IMF在EMD分解后，存在模态混叠问题，由此也导致EMD在信号 去噪领域的推广和应用受到了很大的限制。

发明内容

本公开的目的在于解决上述微地震信号采集过程中不可避免地混合 有非平稳随机噪声，而传统的线性滤波和频谱分析方法对这种混合信号的 去噪效果并不理想，进而导致微震数据信噪比较低的问题。

为了解决上述技术问题，本公开提出一种集成经验模态分解的样本熵 阈值微地震信号降噪方法。

具体地，所述一种集成经验模态分解的样本熵阈值微地震信号降噪方 法包括：

步骤S1，获取含有噪声的微地震信号，并对所述微地震信号进行集成 经验模态分解，得到一系列不同频率成分的本征模态函数IMF；

步骤S2，基于样本熵从所述步骤S1得到的本征模态函数中提取符合 微地震信号特征的目标本征模态函数；

步骤S3，对于所述目标本征模态函数进行重构，得到去噪后的微地震 信号。

在本公开一实施例中，所述步骤S1包括以下步骤：

步骤S11，为所述微地震信号x(n)增加高斯白噪声w(n)，得到目标信号 X(n)；

步骤S12，使用经验模态分解方法对于所述目标信号X(n)进行分解， 得到所述目标信号X(n)的分解表示；

步骤S13，重复上述步骤S11和S12 N次，每次将不同的高斯白噪 声添加到所述微地震信号中，得到N个目标信号序列的分解表示X_i(n)；

步骤S14，对于所述N个目标信号序列的分解表示中的本征模态函数 IMF序列进行平均运算，得到平均本征模态函数IMF序列imf_j′；

步骤S15，基于所述平均本征模态函数IMF序列imf_j′，对于所述微地 震信号x(n)重新进行表示，得到一系列不同频率成分的本征模态函数IMF。

在本公开一实施例中，所述步骤S11中，目标信号X(n)表示为：

X(n)＝x(n)+w(n)，

其中，n表示采样点数。

在本公开一实施例中，所述步骤S12中，所述目标信号X(n)分解表示 为：

其中，k表示IMF序列的数量，imf_j表示一系列本征模态函数IMF序列，r表示残余量。

在本公开一实施例中，所述步骤S13中，第i个目标信号序列的分解 表示为：

在本公开一实施例中，所述步骤S14中，平均本征模态函数IMF序列 imf_j′表示为：

在本公开一实施例中，所述步骤S15中，所述微地震信号x(n)重新表 示为：

其中，r_k(n)表示残余序列。

在本公开一实施例中，所述步骤S2包括以下步骤：

步骤S21，基于所述本征模态函数构造得到一组m维向量：

imf_jm′(1),imf_jm′(2),…,imf_jm′(N-m+1)，

其中，imf_j′(i)为分解分量IMF，j＝1,2,...,k表示模态分量的序号，i＝1,2,...,N表 示数据的长度，imf_jm′(i)＝{imf_j′(i),imf_j′(i+1),...,imf_j′(i+m-1)}；

步骤S22，对于给定的imf_jm′(i)，统计imf_jm′(i)与imf_jm′(l)之间距离小于d的l的数目，并记作B_i，同时计算B_i与总距离的比值

其中，向量imf_jm′(i)与imf_jm′(l)之间的距离d[imf_jm′(i),imf_jm′(l)]为两者对应元素中最 大差值的绝对值，即d[imf_jm′(i),imf_jm′(l)]＝max_{k＝0,...,m-1}(|imf_j′(i+k)-imf_j′(l+k)|)；

步骤S23，计算比值

对于所有i的均值，得到第一比值均值 B^(m)(d)：

步骤S24，将维度m增加到m+1，重复所述步骤S21-S23得到比值

和第二比值均值A^m+1(d)：

其中，A_i表示当维度m增加到m+1时，对于给定的imf_jm′(i)，统计得到的imf_jm′(i)与imf_jm′(l)之间距离小于d的l的数目，

表示A_i与总距离的比值，A^m+1(d) 表示比值

对于所有i的均值；

步骤S25，计算第一比值均值B^(m)(d)和第二比值均值A^m+1(d)的样本熵 SampEn(N,m,r)：

步骤S26，设置样本熵阈值，将样本熵小于所述样本熵阈值的本征模 态函数作为符合微地震信号特征的目标本征模态函数。

本公开能够利用EEMD方法有效解决微地震信号中存在的模态混叠 现象，进而有效去除微地震信号中的噪声，提高微震数据的信噪比。

附图说明

图1是根据本公开一实施例的一种集成经验模态分解的样本熵阈值微 地震信号降噪方法流程图；

图2为根据本公开一实施例的微地震信号的EEMD分解结果示意图；

图3为根据本公开一实施例的IMF样本熵值示意图；

图4(a)为根据本公开一实施例的原始信号分量示意图，图4(b) 为根据本公开一实施例的使用所述去噪方法处理后得到的有效信号分量 示意图；

图5为根据本公开一实施例的实测微地震数据频谱分析与时频分析示 意图；

图6(a)是根据本公开一实施例的利用所述去噪方法对于实测微地震数 据去噪后的效果示意图；

图6(b)是根据本公开一实施例的利用小波多尺度分析方法对于实测微 地震数据去噪后的效果示意图。

具体实施方式

下文中，将参考附图详细描述本公开的示例性实施方式，以使本领域 技术人员可容易地实现它们。此外，为了清楚起见，在附图中省略了与描 述示例性实施方式无关的部分。

在本公开中，应理解，诸如“包括”或“具有”等的术语旨在指示本说明 书中所公开的特征、数字、步骤、行为、部件、部分或其组合的存在，并 且不欲排除一个或多个其他特征、数字、步骤、行为、部件、部分或其组 合存在或被添加的可能性。

另外还需要说明的是，在不冲突的情况下，本公开中的实施例及实施 例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本公 开。

针对上述微地震信号采集过程中不可避免地混合有非平稳随机噪声， 传统的线性滤波和频谱分析方法对这种混合信号的去噪效果并不理想，进 而导致微震数据信噪比较低的问题，本公开提出一种新的集成经验模态分 解的样本熵阈值微地震信号降噪方法。

图1是根据本公开一实施例的一种集成经验模态分解的样本熵阈值微 地震信号降噪方法流程图，如图1所示，所述集成经验模态分解的样本熵 阈值微地震信号降噪方法包括以下步骤：

步骤S1，获取含有噪声的微地震信号，并对所述微地震信号进行集成 经验模态分解(EEMD)，得到一系列不同频率成分的本征模态函数(IMF)；

进一步地，所述步骤S1包括以下步骤：

步骤S11，为所述微地震信号x(n)增加高斯白噪声w(n)，得到目标信号 X(n)：

X(n)＝x(n)+w(n) (1)

其中，n表示采样点数。

步骤S12，使用经验模态分解EMD方法对于所述目标信号X(n)进行分 解，得到所述目标信号X(n)的分解表示：

其中，k表示IMF序列的数量，imf_j表示一系列本征模态函数IMF序列， r表示残余量。

步骤S13，重复上述步骤S11和S12N次，每次将不同的高斯白噪声 添加到所述微地震信号中，得到N个目标信号序列的分解表示X_i(n)：

步骤S14，根据高斯白噪声频谱零均值理论，对于所述N个目标信号 序列的分解表示中的本征模态函数IMF序列进行平均运算，以显著减小多 次加入白噪声对真实IMF的影响，得到平均本征模态函数IMF序列imf_j′：

步骤S15，基于所述平均本征模态函数IMF序列imf_j′，对于所述微地 震信号x(n)重新进行表示，得到一系列不同频率成分的本征模态函数 (IMF)：

其中，r_k(n)表示残余序列。

EEMD是一种噪声辅助信号分析方法，其是在EMD方法的基础上发 展而来的，它通过在待测信号中加入高斯白噪声来平衡信号，并经过多次 EMD分解得到IMF的平均分量，该方法能够有效地解决传统EMD方法 存在的模态混叠问题。上文对于EEMD方法的使用，体现在噪声实现次数 和增加白噪声幅值这两个参数设置上，在本公开一实施例中，将增加的白噪声幅度设置为所述微地震信号的0.2倍，噪声实现次数N设置为100次， 单个IMF分量的最大筛选迭代次数为100次，在此基础上，可将信噪比为 -5dB主频为120Hz的微地震信号自适应地分解为11个分量，如图2所示， 图2中，所述微地震信号在时间尺度下，由高频到低频依次分解，其中， IMF1频率最高，IMF2次之，……，Residual表示频率最低的残余序列。

步骤S2，基于样本熵从所述步骤S1得到的本征模态函数中提取符合 微地震信号特征的目标本征模态函数；所述微地震信号通过上述过程自适 应地分解为了若干个IMF，为了达到信号降噪的目的，该步骤通过计算各 个IMF的样本熵，根据预先设置的样本熵阈值来提取符合微地震信号特征 的目标IMF，后续对于提取得到的目标IMF进行信号重构，即可得到去除 噪声后的微地震信号。

样本熵(SampEn)是一种基于近似熵(ApEn)改进的时间序列复杂度度 量方法，其已被用于评估生理时间序列的复杂性和诊断病理状态。在本公 开实施例中，利用样本熵作为区分噪声信号与微地震信号的判断标准，即 根据样本熵的大小来判断信号的随机性，进而得出噪声与地震信号的分界 点。

进一步地，所述步骤S2包括以下步骤：

步骤S21，基于所述本征模态函数构造得到一组m维向量：

imf_jm′(1),imf_jm′(2),…,imf_jm′(N-m+1) (6)

其中，imf_j′(i)为EEMD的分解分量IMF，j＝1,2,...,k表示模态分量的序号， i＝1,2,...,N表示数据的长度，imf_jm′(i)＝{imf_j′(i),imf_j′(i+1),...,imf_j′(i+m-1)}；

步骤S22，对于给定的imf_jm′(i)，统计imf_jm′(i)与imf_jm′(l)之间距离小于d的 l(1≤l≤N-m+1,l≠i)的数目，并记作B_i，同时计算B_i与总距离的比值

其中，向量imf_jm′(i)与imf_jm′(l)之间的距离d[imf_jm′(i),imf_jm′(l)]为两者对应元素中最 大差值的绝对值，即d[imf_jm′(i),imf_jm′(l)]＝max_{k＝0,...,m-1}(|imf_j′(i+k)-imf′_j(l+k)|)。

步骤S23，计算比值

对于所有i的均值，得到第一比值均值B^(m)(d)：

步骤S24，将维度m增加到m+1，重复所述步骤S21-S23得到比值

和第二比值均值A^m+1(d)：

其中，类似地，A_i表示当维度m增加到m+1时，对于给定的imf_jm′(i)，统计得 到的imf_jm′(i)与imf_jm′(l)之间距离小于d的l(1≤l≤N-m+2,l≠i)的数目，

表示 A_i与总距离的比值，A^m+1(d)表示比值

对于所有i的均值。

步骤S25，计算第一比值均值B^(m)(d)和第二比值均值A^m+1(d)的样本熵 SampEn(N,m,r)；

理论上，第一比值均值B^(m)(d)和第二比值均值A^m+1(d)的样本熵 SampEn(N,m,r)可表示为：

实际上，N不可能是∞，因此，当N为一个有限值时，样本熵可以用下 式估计：

从定义中可以看出，A^m+1(d)不大于B^(m)(d)，所以SampEn(N,m,r)的值非零 即正。由于样本熵主要用于衡量信号复杂度，对于时间序列数据集， SampEn(N,m,r)值越小，表示自相似性越高，噪声越小。在本公开一实施例中， 设置维数m＝1，相似容限d＝0.15×std(x(n))，其中，std(x(n))表示信号x(n)的 标准差，对上述IMF分量进行样本熵计算，得到的计算结果如图3所示， 从图3中可以看出，IMF分量的样本熵随着IMF分量阶数的增大而减小，这说明IMF分量的阶数越高其复杂度越低，所含的噪声也依次减少。因此， 为了能准确地提取有效的微地震信号，可计算不同IMF分量的样本熵，借 助样本熵对于含有大量噪声的信号分量进行筛选。

步骤S26，设置样本熵阈值，将样本熵小于所述样本熵阈值的本征模 态函数作为符合微地震信号特征的目标本征模态函数。

鉴于集成经验模态分解方法在不需要信号先验知识的前提下，能够自 适应地将原信号分解为一系列频率由高到低依次排列的本征模态函数 (IMF)的特点，对低信噪比的微地震数据进行EEMD处理，得到IMF分量， 后续即可通过对于上述IMF分量进行用于衡量信号复杂度的样本熵计算， 来进行噪声信号的筛选，实现正常微震信号与噪声的区分。在本公开一实 施例中，所述样本熵阈值可设置为：δ＝0.2，即后续可对于所有样本熵值小 于0.2的IMF分量进行重组。

步骤S3，对于所述目标本征模态函数进行重构，得到去噪后的微地震 信号。

通过对于符合微地震信号特征的目标本征模态函数进行信号重构，可 得到更贴近于原始微地震信号的信号，即去噪后的微地震信号，由此达到 抑制随机噪声的目的。

为了测试上述去噪方法的可行性，分别对上述原始信号和重构信号进 行了频谱分析和时频分析。图4(a)为根据本公开一实施例的原始信号分 量示意图，对其进行频谱分析和时频分析发现，原始信号分量的主频是 120Hz，其中包括部分低频噪声和大量高频噪声；图4(b)为根据本公开 一实施例的使用所述去噪方法处理后得到的有效信号分量示意图，可以看 出重构后的信号主频与原始数据保持着一定的一致性，从时频谱中可以看 出重构信号与原始信号的低频能量也具有一致性，且重构后的信号有效地 避免了相位畸变(模拟信号的起跳点预设为第1000个采样点)，另一个重 要特性是能够有效保存到达波的振幅信息，振幅信息对震源机制等的研究 非常重要。

接下来对实测的微地震信号进行去噪处理，所选实测数据为某在建水 电站岩爆监测数据中的一次岩爆数据，由于当地特殊的地表条件和复杂多 变的表层结构，采集的微地震资料中存在多种特殊干扰波。如图5所示， 根据对实测的原始监测数据的频谱分析以及时频分析发现，有效信号分量 具有一定的频率范围，而且噪声一部分分布在低频端和高频端，另一部分 与微地震有效信号的频率范围重叠，150Hz附近存在一个强能量的伪单频干扰噪声。分别采用上述去噪方法、小波多尺度分析法分别对此岩爆数据 进行去噪处理。如图6所示，实验中的两种方法都有效地去除了噪声，提 高了信噪比，其中，图6(b)是根据本公开一实施例的利用小波多尺度分析 方法对于实测微地震数据去噪后的效果示意图，如图6(b)所示，该方法虽 能有效保持有效信号的振幅信息，但对伪单频这种常见的干扰波没有很好 的压制；图6(a)是根据本公开一实施例的利用所述去噪方法对于实测微地 震数据去噪后的效果示意图，如图6(a)所示，该方法在保持信号的振幅信 息的基础上能有效地剔除高频噪声和削弱伪单频噪声干扰，且信号的起跳 更加明显，更利于信号的初至拾取和震源机制的研究。

以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行 了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而 已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修 改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (8)

1.一种集成经验模态分解的样本熵阈值微地震信号降噪方法，其特征在于，所述方法包括：

步骤S1，获取含有噪声的微地震信号，并对所述微地震信号进行集成经验模态分解，得到一系列不同频率成分的本征模态函数IMF；

步骤S2，基于样本熵从所述步骤S1得到的本征模态函数中提取符合微地震信号特征的目标本征模态函数；

步骤S3，对于所述目标本征模态函数进行重构，得到去噪后的微地震信号。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1包括以下步骤：

步骤S11，为所述微地震信号x(n)增加高斯白噪声w(n)，得到目标信号X(n)；

步骤S12，使用经验模态分解方法对于所述目标信号X(n)进行分解，得到所述目标信号X(n)的分解表示；

步骤S13，重复上述步骤S11和S12 N次，每次将不同的高斯白噪声添加到所述微地震信号中，得到N个目标信号序列的分解表示X_i(n)；

步骤S14，对于所述N个目标信号序列的分解表示中的本征模态函数IMF序列进行平均运算，得到平均本征模态函数IMF序列imf_j′；

步骤S15，基于所述平均本征模态函数IMF序列imf_j′，对于所述微地震信号x(n)重新进行表示，得到一系列不同频率成分的本征模态函数IMF。

3.根据权利要求2所述的方法，其特征在于，所述步骤S11中，目标信号X(n)表示为：

X(n)＝x(n)+w(n)，

其中，n表示采样点数。

4.根据权利要求3所述的方法，其特征在于，所述步骤S12中，所述目标信号X(n)分解表示为：

其中，k表示IMF序列的数量，imf_j表示一系列本征模态函数IMF序列，r表示残余量。

5.根据权利要求4所述的方法，其特征在于，所述步骤S13中，第i个目标信号序列的分解表示为：

6.根据权利要求5所述的方法，其特征在于，所述步骤S14中，平均本征模态函数IMF序列imf_j′表示为：

7.根据权利要求6所述的方法，其特征在于，所述步骤S15中，所述微地震信号x(n)重新表示为：

其中，r_k(n)表示残余序列。

8.根据权利要求1所述的方法，其特征在于，所述步骤S2包括以下步骤：

步骤S21，基于所述本征模态函数构造得到一组m维向量：

imf_jm′(1),imf_jm′(2),…,imf_jm′(N-m+1)，

其中，imf_j′(i)为分解分量IMF，j＝1,2,...,k表示模态分量的序号，i＝1,2,...,N表示数据的长度，imf_jm′(i)＝{imf_j′(i),imf_j′(i+1),...,imf_j′(i+m-1)}；

步骤S22，对于给定的imf_jm′(i)，统计imf_jm′(i)与imf_jm′(l)之间距离小于d的l的数目，并记作B_i，同时计算B_i与总距离的比值

其中，向量imf_jm′(i)与imf_jm′(l)之间的距离d[imf_jm′(i),imf_jm′(l)]为两者对应元素中最大差值的绝对值，即d[imf_jm′(i),imf_jm′(l)]＝max_{k＝0,...,m-1}(|imf_j′(i+k)-imf_j′(l+k)|)；

步骤S23，计算比值

对于所有i的均值，得到第一比值均值B^(m)(d)：

步骤S24，将维度m增加到m+1，重复所述步骤S21-S23得到比值

和第二比值均值A^m+1(d)：

其中，A_i表示当维度m增加到m+1时，对于给定的imf_jm′(i)，统计得到的imf_jm′(i)与imf_jm′(l)之间距离小于d的l的数目，

表示A_i与总距离的比值，A^m+1(d)表示比值

对于所有i的均值；

步骤S25，计算第一比值均值B^(m)(d)和第二比值均值A^m+1(d)的样本熵SampEn(N,m,r)：

步骤S26，设置样本熵阈值，将样本熵小于所述样本熵阈值的本征模态函数作为符合微地震信号特征的目标本征模态函数。

Images