CN111709279B - 一种利用svd-emd算法分离微地震噪声混合信号的算法 - Google Patents

Abstract

本发明公开了一种利用SVD‑EMD算法分离微地震噪声混合信号的算法，该算法包括(S1)用奇异值分解法SVD对微地震噪声混合信号做预处理；(S2)得到经SVD预处理后的信号f(x,t)，对该信号进行EMD分解，通过初始化残余函数res(x,t)计算所有极值点的当前步骤信号f₀(x,t)，然后从原始信号中除去平均包络值，让剩余部分信号f₁(x,t)成为下一步的信号；(S3)在当前步骤信号f₁(x,t)和上一步骤信号f₀(x,t)之间将它们置于停止状态，并通过固有模式函数和最终剩余功能重复此筛选过程；(S4)利用希尔伯特变换得到通过EMD分解的各个固有模式函数的频率特性，并计算出相关系数R_i从而完成信号的重构。通过上述方案，本发明达到了对微地震震相识别得到提升，拾取P波效果最佳目的，具有很高的实用价值和推广价值。

Description

一种利用SVD-EMD算法分离微地震噪声混合信号的算法

技术领域

本发明具体涉及一种利用SVD-EMD算法分离微地震噪声混合信号的算法。

背景技术

微地震监测系统采集到的微地震资料都是一些微弱的信号，导致微地震信号的信号特征不易被有效识别出来或者淹没在噪声信号中，容易被监测人员忽略。因为微地震信号具有不易被监测到的特性，使得微地震监测系统中传感器的灵敏度一般都设置较高，这样造成监测得到的数据量特别庞大，使得监测人员很难及时的从这些监测信号中筛选出真正的微地震信号，因此微地震噪声分离是微地震检测技术的基础，为岩爆、矿震等动力灾害预警提高可靠资料。

目前已经有许多种微地震事件辨识方法，Allen等提出了长短时均值比方法(STA/LTA)，该方法的原理是在微地震事件到来时，短时窗能量值迅速增大,导致STA/LTA比值瞬间变大，该方法算法简单，辨识效率高，但该方法需要确定多个参数，不易掌控，对低信噪比信号的辨识效果较差，且经常辨识错误；Maeda等基于地震波形在地震事件到达前后数据统计性质具有显著差别这一理论，提出了赤池信息准则(AIC方法)，同样该方法算法简单，辨识效率高，但该方法计算量大；Zhang Haijing等基于小波变换和AIC方法提出了W-AIC方法该方法首先利用小波变换将信号分解为不同频率段的波形，然后在各个分频波形上应用AIC方法进行到时拾取，如果各个分频波形的拾取结果一致，则认为存在有效微地震事件，虽然该方法克服了单一AIC拾取的缺陷，但是该算法复杂度高；在2002年，Saragiotis等提出了新的拾取方法，称为PAI-S/K方法，该方法的原理依据是地震波形偏斜度和峰度理论，该方法识别效率较高，但该方法受时窗大小影响很大，当时窗长度选取不当时，事件辨识就会出错。可以看出以上方法，均能有效实现微地震事件的识别与分离，但是在算法精度和算法计算量方面都存在不同程度的缺陷，降低了微震数据识别的准确性与可行性。

发明内容

为了克服现有技术中的上述不足，本发明提供一种利用SVD-EMD算法分离微地震噪声混合信号的算法，能有效分离出随机不平稳信号以及不规则的干扰信号,且用EMD分解得到的波形图不易失真。

为了实现上述目的，本发明采用的技术方案如下：

一种利用SVD-EMD算法分离微地震噪声混合信号的算法，包括如下步骤：

(S1)用奇异值分解法SVD对微地震噪声混合信号做预处理；

(S2)得到经SVD预处理后的信号f(x,t)，对该信号进行EMD分解，求得残余函数res(x,t)，计算当前步骤信号f₀(x,t)所有极值点，采用双三样条插值函数计算所有范围函数的上包络up(x,t)和低包络down(x,t)，并获得两个包络的平均包络值avr(x,t)，然后从原始信号中除去平均包络值，让剩余部分信号f₁(x,t)成为下一步的信号；

(S3)在当前步骤信号f₁(x,t)和上一步骤信号f₀(x,t)之间将它们置于停止状态，如果满足停止条件，则从原始信号中移除分解的固有模式函数IMF，并将剩余部分重新初始化为新的残余函数和当前值的阶跃信号，如果不满足停止条件，则以信号f₁(x,t)为初始信号重复步骤(S2)以获取后续步骤的信号直到满足停止条件，并通过固有模式函数和最终残余分量重复此筛选过程；

(S4)利用希尔伯特变换得到通过EMD分解的各个固有模式函数的频率特性，并计算出相关系数R_i从而完成信号的重构，通过相关系数R_i与原信号分解的固有模式函数IMF信号，得到原信号平均值和所有固有模式函数IMF平均值有关。

进一步地，所述算法还包括步骤(S5)：利用步骤(S3)中分解出来的有效的固有模式函数IMF_P和IMF_N重构微地震信号，在进行SVDEMD分解得到新的IMF分量层c’(t)，然后得到信噪比高的有效微地震信号。

进一步地，所述步骤(S1)中对微地震噪声混合信号做预处理是利用微地震信号s(t)中的三个分量X(t)，Y(t)，Z(t)构造出一个m×n的Hankel矩阵A，对该矩阵进行SVD分解，得到互不相关的噪声子空间和信号子空间，然后利用SVD分解求得Hankel矩阵A的奇异值，保留对矩阵内有效信号的特征值，把其余的特征值置为零，初步去除高频信号中噪声子空间的噪声。

进一步地，所述步骤(S1)中SVD分解的计算公式为：

其中，U表示M×N阶正交矩阵；D表示对角元素；D＝diag(σ₁,σ₂,σ₃......σ_n)表示M×N阶的对角矩阵；V表示N×N阶正交矩阵；σ_i表示矩阵A的奇异值；r表示矩阵A的秩，T表示转置；u_i表示矩阵AA^T的第i个特征向量；v_i表示矩阵A^TA的第i个特征向量，P代表有效信号，N代表噪声信号，L代表信号矩阵长度，i就是一个变量；同时，M×N、N×N表示m×n中的一部分。

进一步地，所述步骤(S2)中平均包络和剩余部分信号的计算公式分别为：

avr(x,t)＝(up+down)/2 (3)

f₁(x,t)＝f₀(x,t)-avr(x,t) (4)

进一步地，所述步骤(S3)中新的残余函数和当前值的阶跃信号的计算公式分别为：

res(x,t)＝f₀(x,t)＝f(x,t)-IMF₁ (5)

具体地，所述步骤(S4)中相关系数R_i的计算公式为：

其中，p(i)为原信号，c(i)为IMF信号，为原信号平均值，/>为IMF信号的平均值。

与现有技术相比，本发明具有以下有益效果：

本发明在对岩石微破裂信号进行特征分解时，能有效从信噪比比较低的微地震混合信号中分解出岩石微破裂的特征信号，从而避免由于信噪比问题影响分离精度；本发明采用奇异经验模态分解(SVEMD)算法对信号进行分解；实验表明，EMD分解在对微地震信号进行分解时，能有效保留信号的局部特征，同时通过仿真实验和矿上微地震信号对该方法进行验证，得到结论。并且该方法与AIC方法相结合，可以对微地震震相识别得到提升，拾取P波效果最佳，进一步提升岩石微破裂的位置精度。

附图说明

图1为本发明实施流程框图。

图2为本发明在合成微地震数据SVEMD分解后的前三个IMF分量图。

图3为本发明在合成微地震数据分量前三个IMF的希尔伯特谱。

图4为本发明实践例微震监测检波器及现场设备布设的示意图。

图5为本发明实践例实测微地震数据SVEMD分解后前三个IMF分量图。

图6为本发明实践例实测微地震数据分量IMF的希尔伯特谱。

图7为本发明采用不同种方法对实践例实测微地震数据分解得到IMF1分量图。

图8为本发明实践例实测地震数据经过SVEMD分解后的IMF分量图。

图9为本发明实践例实测地震数据分量IMF的希尔伯特谱。

图10为本发明采用不同种方法对实践例地震数据分解得到IMF1分量图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。

实施例

如图1至图10所示，一种利用SVD-EMD算法分离微地震噪声混合信号的算法，包括以下步骤：

(1)用奇异值分解法(SVD)对微地震噪声混合信号做预处理。首先利用微地震信号s(t)中的三分量X(t)，Y(t)，Z(t)构造出一个m×n的Hankel矩阵A，对该矩阵进行SVD分解，得到互不相关的噪声子空间和信号子空间；接着利用SVD分解求得Hankel矩阵A的奇异值，保留对角阵内有效信号的特征值，把较小的特征值置为零，初步去除高频信号中的噪声。所述SVD分解的计算公式为：

其中，U表示M×N阶正交矩阵；D表示对角元素；D＝diag(σ₁,σ₂,σ₃......σ_n)表示M×N阶的对角矩阵；V表示N×N阶正交矩阵；σ_i表示矩阵A的奇异值；r表示矩阵A的秩，T表示转置；u_i表示矩阵AA^T的第i个特征向量；v_i表示矩阵A^TA的第i个特征向量。

(2)对于经SVD算法预处理的信号f(x,t)，首先初始化残余函数res(x,t)和当前步骤信号f₀(x,t)。并计算所有极值点当前信号f₀(x,t)：采用双三样条插值函数来计算所有局部的上包络up(x,t)和低包络down(x,t)，获得两个包络的平均包络值avr(x,t)。然后从原始信号中除去平均包络值，让剩余部分成为下一步的信号。所述新的残余函数和当前值的阶跃信号为：

avr(x,t)＝(up+down)/2 (3)

f₁(x,t)＝f₀(x,t)-avr(x,t) (4)

(3)在当前步信号f₁和前一步信号f₀之间将它们置于停止状态。如果停止条件不满足，重复此拟合并减去平均过程以获得f₂,f₃,......直到满足停止条件。若满足条件，就从原始信号中移除分解的固有模式函数IMF₁，并将剩余部分重新初始化为新的残余函数和当前值的阶跃信号；对IMF₂，IMF₃.....和最终剩余功能重复此筛选过程。

res(x,t)＝f₀(x,t)＝f(x,t)-IMF₁ (5)

(4)利用希尔伯特变换得到通过EMD分解得到的各个IMF的频率特性，并计算出相关系数R_i，进而完成信号重构。EMD分解加上HAS分析可以解决大部分非平稳随机信号。

其中，p(i)为原信号，c(i)为IMF信号，为原信号平均值，/>为IMF信号的平均值。

(5)通过这几步算法得到相对有效的信号，算法结束，这一步可以分解出微地震和噪声信号。如果信号分解效果不好，接下来还可以利用步骤(S3)中分解出来的有效的固有模式函数IMF_P和IMF_N重构微地震信号，在进行SVDEMD分解得到新的IMF分量层c’(t)，然后得到信噪比高的有效微地震信号。

为了更好地理解本发明，下面采用实验例对本发明做进一步的阐述：

利用matlab软件产生振动数据特点，将本算法进行应用，得到的结果进行验证说明，对比不同算法，证明本算法的优越性。

本实例采用数据合成的三分量微地震数据，主要利用人工模拟得到“3000*3”类型的微地震数据，采样频率6kHz。同时在原始微地震数据加入不同噪声，含有随机噪声的微地震数据和含有高斯白噪声信号的微地震数据。采用几种不同幅度的随机噪声和高斯白噪声数据与微地震数据进行合成，形成多个微地震信号叠加的数据，形成的还有不同信噪比的数据，该数据采样点是6000个点。分别用A、B、C表示三分量合成数据，3组合成数据均由一个按照指数规律衰减的正弦波叠加不同幅值的随机噪声产生。为了保持和实际微震数据类似，合成数据的采样率也设置为6kHz，数据长度定义为3000个采样点，合成的2对数据信噪比分别为11.4dB和7.77dB。合成地震数据的主频率在87Hz。将本发明的方法应用对合成数据的分解过程中，图2显示的利用SVEMD对合成数据进行分解的结果，主要选取前三个IMF分量，图3为对IMF分量进行HAS变换得到的频谱图，同时估计本发明计算时间。表1为合成数据进行SVEMD分离后的IMF的幅度，表2为合成数据进行SVEDM分离后的不同IMF的频率，表3为本发明对比小波分解，EMD分解的计算时间情况(单位：秒)。

表1合成数据采用SVEMD分解的幅度(dB)

表2合成数据采用SVEMD分解的频率(Hz)

表3合成数据采用不同方法的计算时间(单位：秒)

下面采用实践例对本发明做更进一步地阐述：

本实践例采用的数据为我国西南某水电站尾水隧洞开挖岩爆监测的微震监测数据。仪器采用IMS的20通道微震监测设备，现场布置了10支微震检波器，其中，单分量检波器和三分量检波器各5支，采样频率为6kHz，主要用于监测隧道内的岩爆微地震信号，检波器及现场设备布设的示意图如图4所示。沿隧道开挖方向进行布设，本实践例选取了不同信噪比的岩爆事件，分别采用本发明、小波4阶分解方法和EMD方法对数据进行分解。其中图5为采用本方法对微地震数据进行分解的结果，图6为利用本方法对分解得到的前3个IMF进行HAS分析后得到的频谱图，图7为不同种方法对微地震数据分解得到的IMF1分量图。表4为采用不同中方法对微地震数据的计算时间。

表4实际微地震数据采用不同方法的计算时间(单位：秒)

下面采用另外一个实践例对本发明做更进一步地阐述：

本次数据来自于澳大利亚麦格理岛在2018年7月13日9时45分38秒发生地震时采集到的，该次地震产生巴伦尼群岛地区，震源深度10公里，震级5.8级，该数据采样频率为500Hz，采样点为27000个，本事件地震数据分别采用本方法，小波4阶分解方法和EMD方法对数据进行分解，其中图8为采用本方法对地震数据进行分解的结果，图9为利用本方法对分解得到前3个IMF进行HAS分析，得到频谱图，图10为不同种方法对地震数据分解得到的IMF1分量图。表5为采用不同种方法对地震数据的计算时间。

表5地震数据采用不同方法的计算时间(单位：秒)

本发明对微地震信号的分离提出了SVD-EMD算法；用此算法顺利进行实验仿真和实例数据处理，采用的实例数据为白鹤滩微地震监测数据和典型地震数据。要分析此算法对微地震信号分离的效果，可定性地从信号的时间频率特点来分析；除此之外，还可以用能量比或相邻分量的互信息熵等定量的方法来分析。本发明在对岩石微破裂信号进行特征分解时，能有效从信噪比比较低的微地震混合信号中分解出岩石微破裂的特征信号，从而避免了由于信噪比的问题影响到分离精度；本发明针对低信噪比数据的分离具有一定的优点，对信噪比高的混合信号也几乎不受信号复杂度的影响，与其他几种改进算法相比，本发明具有更高的自动化程度和准确度。

上述实施例仅为本发明的优选实施例，并非对本发明保护范围的限制，但凡采用本发明的设计原理，以及在此基础上进行非创造性劳动而做出的变化，均应属于本发明的保护范围之内。

Claims (7)

1.一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，包括如下步骤：

(S1)用奇异值分解法SVD对微地震噪声混合信号做预处理；

(S2)得到经SVD预处理后的信号f(x,t)，对该信号进行EMD分解，求得残余函数res(x,t)，计算当前步骤信号f₀(x,t)所有极值点，采用双三样条插值函数计算所有信号的上包络up(x,t)和低包络down(x,t)，并获得两个包络的平均包络值avr(x,t)，然后从原始信号中除去平均包络值，让剩余部分信号f₁(x,t)成为下一步的信号；

(S3)在当前步骤信号f₁(x,t)和上一步骤信号f₀(x,t)之间将它们置于停止状态，如果满足停止条件，则从原始信号中移除分解的固有模式函数IMF，并将剩余部分重新初始化为新的残余函数和当前值的阶跃信号，如果不满足停止条件，则以信号f₁(x,t)为初始信号重复步骤(S2)以获取后续步骤的信号直到满足停止条件，并通过固有模式函数和最终残余分量R(S)重复此筛选过程；

(S4)利用希尔伯特变换得到通过EMD分解的各个固有模式函数的频率特性，并计算出相关系数R_i从而完成信号的重构，通过相关系数R_i与原信号分解的固有模式函数IMF信号，得到原信号平均值和所有固有模式函数IMF平均值有关。

2.根据权利要求1所述的一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，所述算法还包括步骤(S5)：利用步骤(S3)中分解出来的有效的固有模式函数IMF_P和IMF_N重构微地震信号，在进行SVDEMD分解得到新的IMF分量层c’(t)，然后得到信噪比高的有效微地震信号。

3.根据权利要求1所述的一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，所述步骤(S1)中对微地震噪声混合信号做预处理是利用微地震信号s(t)中的三个分量X(t)，Y(t)，Z(t)构造出一个m×n的Hankel矩阵A，对该矩阵进行SVD分解，得到互不相关的噪声子空间和信号子空间，然后利用SVD分解求得Hankel矩阵A的奇异值，保留对矩阵内有效信号的特征值，把其余的特征值置为零，初步去除高频信号中噪声子空间的噪声。

4.根据权利要求3所述的一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，所述步骤(S1)中SVD分解的计算公式为：

其中，U表示M×N阶正交矩阵；D表示对角元素；D＝diag(σ₁,σ₂,σ₃......σ_n)表示M×N阶的对角矩阵；V表示N×N阶正交矩阵；σ_i表示矩阵A的奇异值；r表示矩阵A的秩，T表示转置；u_i表示矩阵AA^T的第i个特征向量；v_i表示矩阵A^TA的第i个特征向量，P代表有效信号，N代表噪声信号，L代表信号矩阵长度，i就是一个变量；同时，M×N、N×N表示m×n中的一部分。

5.根据权利要求1或2所述的一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，所述步骤(S2)中平均包络和剩余部分信号的计算公式分别为：

avr(x,t)＝(up+down)/2 (3)

f₁(x,t)＝f₀(x,t)-avr(x,t) (4)。

6.根据权利要求1所述的一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，所述步骤(S3)中新的残余函数和当前值的阶跃信号的计算公式分别为：

res(x,t)＝f₀(x,t)＝f(x,t)-IMF₁ (5)

7.根据权利要求1所述的一种利用SVD-EMD算法分离微地震噪声混合信号的算法，其特征在于，所述步骤(S4)中相关系数R_i的计算公式为：

其中，p(i)为原信号，c(i)为IMF信号，为原信号平均值，/>为IMF信号的平均值。