CN110514294A - 一种基于emd与vmd的爆破振动信号降噪方法 - Google Patents

Abstract

本发明公开一种基于EMD与VMD的爆破振动信号降噪方法，首先对原始爆破振动信号进行EMD分解，得到本征模态函数及其数量n，然后求得原始爆破振动信号和每个本征模态函数的功率谱密度，找出每个本征模态函数功率谱密度的最大值及对应频率，判断在此频率下本征模态函数功率谱密度与原始爆破振动信号功率谱密度的比值，若比值小于10％，则视为噪声信号，记噪声信号个数为j。最后对原始爆破振动信号进行VND分解，分解个数K为n‑j，将中心频率低于10Hz及高于200Hz的分量滤除，重构剩余信号。

Description

一种基于EMD与VMD的爆破振动信号降噪方法

技术领域

本发明属于信号处理技术领域，涉及一种基于EMD与VMD的爆破振动信号降噪方法。

背景技术

爆破振动信号具有持续时间短、突变快等特点，是一种短时非平稳随机信号。从测量数据中有效的除去噪声干扰，提高信号的分辨率，是正确分析爆破振动信号的前提。对进一步找出控制爆破振动大小的因素，降低爆破振动危害有着重要意义。

经验模态分解(Empirical Mode Decomposition，EMD)可以根据信号本身将非平稳信号分解成一系列本征模态函数，分解中无需信号先验知识，是后验地、完全自适应地分解。但在信号分解的样条插值中有明显的端点效应问题，而且分解后的一系列本征模态函数存在模态混淆问题。

变分模态分解(Variational Mode Decomposition，VMD)将信号分解为一系列有限带宽变分模态分量，将信号转移到变分框架内进行分解，通过寻找模型最优解实现信号分解。VMD方法的分解正确性依赖于分解模态个数K，K是VMD信号分解的关键。

发明内容

为解决上述技术问题，本发明的目的是提供一种基于EMD与VMD的爆破振动信号降噪方法，以克服EMD方法分解后本征模态函数混淆的问题，又解决了VMD方法分解模态数未知的问题，为后续分析提供准确条件。

本发明提供一种基于EMD与VMD的爆破振动信号降噪方法，包括如下步骤：

步骤1：用测振仪对矿山原始爆破振动信号进行采集；

步骤2：对采集的原始爆破振动信号进行EMD分解，得到所有本征模态函数分量IMF_i(i＝1,2,…,n)及数量n；

步骤3：求出各本征模态函数分量IMF_i的功率谱密度及原始爆破振动信号的功率谱密度，根据本征模态函数分量IMF_i的功率谱密度和原始爆破振动信号的功率谱密度的比值判断噪声数量；

步骤4：根据本征模态函数分量的数量和噪声数量计算VMD分解的变分模态分量个数，对采集的原始爆破振动信号进行VMD分解，将中心频率低于10Hz及高于200Hz的变分模态分量滤除，重构剩余分量。

在本发明的基于EMD与VMD的爆破振动信号降噪方法中，所述步骤2中对采集的原始爆破振动信号进行EMD分解具体为：

步骤2.1：找出原始爆破振动信号X(t)的所有极大值和极小值，分别用三次样条函数曲线对所有极大值和极小值进行插值，拟合出原始爆破振动信号X(t)的上下两条包络线，按顺序连接上、下两条包络线的均值得到均值线m₁(t)；

其中，X_max(t)为原始爆破振动信号的上包络线，X_min(t)为原始爆破振动信号的下包络线；

步骤2.2：用原始爆破振动信号X(t)减掉均值线m₁(t)得到h₁(t)，即：

h₁(t)＝X₁(t)-m₁(t) (2)

步骤2.3：将h₁(t)当作原始信号，重复步骤2.1-2.2，获得h₁₁(t)：

h₁₁(t)＝h₁(t)-m₁₁(t) (3)

其中，m₁₁(t)是连接h₁(t)的上、下包络线的均值得到均值线；

步骤2.4：判断h₁₁(t)是否为本征模态函数分量，如果是，记IMF₁＝h₁₁(t)，执行步骤2.5；如果不是将h₁₁(t)当作原始信号继续重复k次步骤2.3，直到满足以下条件，并记IMF₁＝h_1k(t)：

其中，SD取0.2～0.3，IMF₁为第一个本征模态函数分量；

步骤2.5：从原始爆破振动信号X(t)中减去第一本征模态函数分量IMF₁得到剩余信号，即第一残差r₁(t)；

r₁(t)＝X(t)-IMF₁ (5)

步骤2.6：将第一残差r₁(t)当作原始信号重复2.1-2.5操作过程，得到第二个本征模态函数分量IMF₂，循环n次得到所有满足要求的本征模态函数分量IMF_i(i＝1,2，…，n)，并记为：

其中r_n为最终残差，是一个单调函数，代表爆破振动信号的整体趋势。

在本发明的基于EMD与VMD的爆破振动信号降噪方法中，所述步骤3中具体为：

找出每个本征模态函数分量IMF_i功率谱密度的最大值为P_max及其对应频率，判断在此频率下，IMF_i的功率谱密度与原始爆破振动信号的功率谱密度的比值，若比值小于10％，则视为噪声信号，计噪声信号个数为j。

在本发明的基于EMD与VMD的爆破振动信号降噪方法中，所述步骤4中对原始爆破振动信号采用VMD分解为若干变分模态分量，使各个变分模态分量的估计带宽之和最小，各变分模态分量之和等于原始信号X(t)，数学公式为：

式(7)中，{u_k}＝{u₁,u₂…u_K}为VMD分解得到的K个变分模态分量，{ω_k}＝{ω₁,ω₂...ω_K}为各自的中心频率，δ(t)为狄拉克函数，*表示卷积，j²＝-1；式(8)中为变分模态分量之和。

在本发明的基于EMD与VMD的爆破振动信号降噪方法中，求式(7)(8)的最优解，引入增广拉格朗日将变分约束问题变为变分非约束问题，其表达式如下：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子。

在本发明的基于EMD与VMD的爆破振动信号降噪方法中，所述步骤4中对采集的原始爆破振动信号进行VMD分解具体包括：

步骤4.1：初始化和n为0；变分模态分量个数K取EMD方法分解的本征模态函数分量IMF_i个数减去噪声信号个数，即K＝n-j；

步骤4.2：通过式(10)(11)更新和

步骤4.3：通过式(12)更新

式中，τ为保真系数；

步骤4.4：重复步骤4.2-4.3直至(13)满足，最终得到K个变分模态分量：

式中，ε为容差；

步骤4.5：将中心频率低于10Hz及高于200Hz的变分模态分量滤除，重构剩余分量。

本发明的一种基于EMD与VMD的爆破振动信号降噪方法，至少具有以下有益效果：

将一种全新信号分解方法引入到爆破振动信号分析方面，克服了EMD方法分解后，本征模态函数混淆的问题，又解决了VMD方法分解模态数未知的问题，为后续分析提供准确条件。

附图说明

图1是本发明的一种基于EMD与VMD的爆破振动信号降噪方法的流程图。

具体实施方式

如图1所示，首先对原始爆破振动信号进行EMD分解，得到本征模态函数及其数量n，然后求得原始爆破振动信号和每个本征模态函数的功率谱密度，找出每个本征模态函数功率谱密度的最大值及对应频率，判断在此频率下本征模态函数功率谱密度与原始爆破振动信号功率谱密度的比值，若比值小于10％，则视为噪声信号，记噪声信号个数为j。最后对原始爆破振动信号进行VND分解，分解个数K为n-j，将中心频率低于10Hz及高于200Hz的分量滤除，重构剩余信号。

实施例一：针对弓长岭露天矿某次实验，单孔爆破对采空区顶板稳定性研究中，对采空区顶板上部地表进行检测。本次检测放置6台测振仪，分别对采空区顶板厚度最小处、采空区顶板正中心及采空区四周进行检测，其中最大振速出现在距离炮孔最近的1号设备，对其Z方向振动图像进行研究。测试设备为TC-4850爆破测振仪、三矢量传感器，参数设置电平为0.1cm/s，采集频率为2000Hz，采集时间为2s，延时为-100ms。

具体工艺实施流程如下：

本发明的一种基于EMD与VMD的爆破振动信号降噪方法，包括如下步骤：

步骤1：用测振仪对矿山原始爆破振动信号进行采集；

步骤2：对采集的原始爆破振动信号进行EMD分解，得到所有本征模态函数分量IMF_i(i＝1,2,…,n)及数量n，EMD分解具体过程为：

步骤2.1：找出原始爆破振动信号X(t)的所有极大值和极小值，分别用三次样条函数曲线对所有极大值和极小值进行插值，拟合出原始爆破振动信号X(t)的上下两条包络线，按顺序连接上、下两条包络线的均值得到均值线m₁(t)；

其中，X_max(t)为原始爆破振动信号的上包络线，X_min(t)为原始爆破振动信号的下包络线；

步骤2.2：用原始爆破振动信号X(t)减掉均值线m₁(t)得到h₁(t)，即：

h₁(t)＝X₁(t)-m₁(t) (2)

步骤2.3：将h₁(t)当作原始信号，重复步骤2.1-2.2，获得h₁₁(t)：

h₁₁(t)＝h₁(t)-m₁₁(t) (3)

其中，m₁₁(t)是连接h₁(t)的上、下包络线的均值得到均值线；

步骤2.4：判断h₁₁(t)是否为本征模态函数分量，如果是，记IMF₁＝h₁₁(t)，执行步骤2.5；如果不是将h₁₁(t)当作原始信号继续重复k次步骤2.3，直到满足以下条件，并记IMF₁＝h_1k(t)：

其中，SD取0.2～0.3，IMF₁为第一个本征模态函数分量；

步骤2.5：从原始爆破振动信号X(t)中减去第一本征模态函数分量IMF₁得到剩余信号，即第一残差r₁(t)；

r₁(t)＝X(t)-IMF₁ (5)

步骤2.6：将第一残差r₁(t)当作原始信号重复2.1-2.5操作过程，得到第二个本征模态函数分量IMF₂，循环n次得到所有满足要求的本征模态函数分量IMF_i(i＝1,2，…，n)，并记为：

其中r_n为最终残差，是一个单调函数，代表爆破振动信号的整体趋势。

步骤3：求出各本征模态函数分量IMF_i的功率谱密度(power spectral density，PSD)及原始爆破振动信号的功率谱密度，根据本征模态函数分量IMF_i的功率谱密度和原始爆破振动信号的功率谱密度的比值判断噪声数量；

具体实施时，找出每个本征模态函数分量IMF_i功率谱密度的最大值为P_max及其对应频率，判断在此频率下，IMF_i的功率谱密度与原始爆破振动信号的功率谱密度的比值，若比值小于10％，则视为噪声信号，计噪声信号个数为j。

步骤4：根据本征模态函数分量的数量和噪声数量计算VMD分解的变分模态分量个数，对采集的原始爆破振动信号进行VMD分解，将中心频率低于10Hz及高于200Hz的变分模态分量滤除，重构剩余分量。

所述步骤4中对原始爆破振动信号采用VMD分解为若干变分模态分量，使各个变分模态分量的估计带宽之和最小，各变分模态分量之和等于原始信号X(t)，数学公式为：

式(7)中，{u_k}＝{u₁,u₂…u_K}为VMD分解得到的K个变分模态分量，{ω_k}＝{ω₁,ω₂...ω_K}为各自的中心频率，δ(t)为狄拉克函数，*表示卷积，j²＝-1；式(8)中为变分模态分量之和。

求式(7)(8)的最优解，引入增广拉格朗日将变分约束问题变为变分非约束问题，其表达式如下：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子。

所述步骤4中对采集的原始爆破振动信号进行VMD分解具体过程为：

步骤4.1：初始化和n为0；变分模态分量个数K取EMD方法分解的本征模态函数分量IMF_i个数减去噪声信号个数，即K＝n-j；

步骤4.2：通过式(10)(11)更新和

步骤4.3：通过式(12)更新

式中，τ为保真系数；

步骤4.4：重复步骤4.2-4.3直至(13)满足，最终得到K个变分模态分量：

式中，ε为容差；

步骤4.5：将中心频率低于10Hz及高于200Hz的变分模态分量滤除，重构剩余分量。

实施例二：针对关宝山台阶爆破时，附近居民反应强烈，因此对距离爆区最近一户居民的房屋质点振速进行检测。本次检测选择了3个测点，分别为距离爆区最近的围墙处，居民主房屋门前及物品存放仓库内。测试过程中采用的仪器设备有TC-4850爆破测振仪、三矢量传感器。参数设置电平为0.05cm/s，采集频率为4000Hz，采集时间为3s，延时为-100ms。

具体工艺实施步骤按上述步骤进行操作。

实施例三：针对多宝山铜矿底部掘沟爆破过程中，考虑爆破振动高程放大作用，防止边坡失稳发生滑坡，对爆区上部四层台阶布置测点进行监测，其中在每层台阶的西北和西南方向处放置测振仪，测试设备选用四川拓普数字设备有限公司开发研制的Topbox振动自记仪。

具体工艺实施步骤按上述步骤进行操作。

实施例四：针对大孤山铁矿爆破施工过程中，其爆区下部有一巷道，为探究爆破施工是否对其稳定性有影响而进行振动监测，由于巷道内环境复杂，只能将测振仪安置在下底板。测振仪分别布置在巷道内标号68、69、70、71点处。测试设备为TC-4850爆破测振仪、三矢量传感器，参数设置电平为0.1cm/s，采集频率为8000Hz，采集时间为2s，延时为-100ms。

以上所述仅为本发明的较佳实施例，并不用以限制本发明的思想，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于EMD与VMD的爆破振动信号降噪方法，其特征在于，包括如下步骤：

步骤1：用测振仪对矿山原始爆破振动信号进行采集；

步骤2：对采集的原始爆破振动信号进行EMD分解，得到所有本征模态函数分量IMF_i(i＝1,2,…,n)及数量n；

步骤3：求出各本征模态函数分量IMF_i的功率谱密度及原始爆破振动信号的功率谱密度，根据本征模态函数分量IMF_i的功率谱密度和原始爆破振动信号的功率谱密度的比值判断噪声数量；

步骤4：根据本征模态函数分量的数量和噪声数量计算VMD分解的变分模态分量个数，对采集的原始爆破振动信号进行VMD分解，将中心频率低于10Hz及高于200Hz的变分模态分量滤除，重构剩余分量。

2.如权利要求1所述的基于EMD与VMD的爆破振动信号降噪方法，其特征在于，所述步骤2中对采集的原始爆破振动信号进行EMD分解具体为：

步骤2.1：找出原始爆破振动信号X(t)的所有极大值和极小值，分别用三次样条函数曲线对所有极大值和极小值进行插值，拟合出原始爆破振动信号X(t)的上下两条包络线，按顺序连接上、下两条包络线的均值得到均值线m₁(t)；

其中，X_max(t)为原始爆破振动信号的上包络线，X_min(t)为原始爆破振动信号的下包络线；

步骤2.2：用原始爆破振动信号X(t)减掉均值线m₁(t)得到h₁(t)，即：

h₁(t)＝X₁(t)-m₁(t) (2)

步骤2.3：将h₁(t)当作原始信号，重复步骤2.1-2.2，获得h₁₁(t)：

h₁₁(t)＝h₁(t)-m₁₁(t) (3)

其中，m₁₁(t)是连接h₁(t)的上、下包络线的均值得到均值线；

步骤2.4：判断h₁₁(t)是否为本征模态函数分量，如果是，记IMF₁＝h₁₁(t)，执行步骤2.5；如果不是将h₁₁(t)当作原始信号继续重复k次步骤2.3，直到满足以下条件，并记IMF₁＝h_1k(t)：

其中，SD取0.2～0.3，IMF₁为第一个本征模态函数分量；

步骤2.5：从原始爆破振动信号X(t)中减去第一本征模态函数分量IMF₁得到剩余信号，即第一残差r₁(t)；

r₁(t)＝X(t)-IMF₁ (5)

步骤2.6：将第一残差r₁(t)当作原始信号重复2.1-2.5操作过程，得到第二个本征模态函数分量IMF₂，循环n次得到所有满足要求的本征模态函数分量IMF_i(i＝1,2，…，n)，并记为：

其中r_n为最终残差，是一个单调函数，代表爆破振动信号的整体趋势。

3.如权利要求1所述的基于EMD与VMD的爆破振动信号降噪方法，其特征在于，所述步骤3中具体为：

找出每个本征模态函数分量IMF_i功率谱密度的最大值为P_max及其对应频率，判断在此频率下，IMF_i的功率谱密度与原始爆破振动信号的功率谱密度的比值，若比值小于10％，则视为噪声信号，计噪声信号个数为j。

4.如权利要求3所述的基于EMD与VMD的爆破振动信号降噪方法，其特征在于，所述步骤4中对原始爆破振动信号采用VMD分解为若干变分模态分量，使各个变分模态分量的估计带宽之和最小，各变分模态分量之和等于原始信号X(t)，数学公式为：

式(7)中，{u_k}＝{u₁,u₂...u_K}为VMD分解得到的K个变分模态分量，{ω_k}＝{ω₁,ω₂...ω_K}为各自的中心频率，δ(t)为狄拉克函数，*表示卷积，j²＝-1；式(8)中为变分模态分量之和。

5.如权利要求4所述的基于EMD与VMD的爆破振动信号降噪方法，其特征在于，求式(7)(8)的最优解，引入增广拉格朗日将变分约束问题变为变分非约束问题，其表达式如下：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子。

6.如权利要求5所述的基于EMD与VMD的爆破振动信号降噪方法，其特征在于，所述步骤4中对采集的原始爆破振动信号进行VMD分解具体包括：

步骤4.1：初始化和n为0；变分模态分量个数K取EMD方法分解的本征模态函数分量IMF_i个数减去噪声信号个数，即K＝n-j；

步骤4.2：通过式(10)(11)更新和

步骤4.3：通过式(12)更新

式中，τ为保真系数；

步骤4.4：重复步骤4.2-4.3直至(13)满足，最终得到K个变分模态分量：

式中，ε为容差；

步骤4.5：将中心频率低于10Hz及高于200Hz的变分模态分量滤除，重构剩余分量。