CN111881858B - 一种微震信号多尺度去噪方法、装置及可读存储介质 - Google Patents

一种微震信号多尺度去噪方法、装置及可读存储介质 Download PDF

Info

Publication number
CN111881858B
CN111881858B CN202010760856.4A CN202010760856A CN111881858B CN 111881858 B CN111881858 B CN 111881858B CN 202010760856 A CN202010760856 A CN 202010760856A CN 111881858 B CN111881858 B CN 111881858B
Authority
CN
China
Prior art keywords
component
signal
denoising
matrix
imf
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010760856.4A
Other languages
English (en)
Other versions
CN111881858A (zh
Inventor
彭康
尚雪义
郭宏扬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing University
Central South University
Original Assignee
Chongqing University
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing University, Central South University filed Critical Chongqing University
Priority to CN202010760856.4A priority Critical patent/CN111881858B/zh
Publication of CN111881858A publication Critical patent/CN111881858A/zh
Application granted granted Critical
Publication of CN111881858B publication Critical patent/CN111881858B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Abstract

本发明公开了一种微震信号多尺度去噪方法、装置及可读存储介质，所述方法包括：步骤1：获取微震信号并进行EMD或EEMD分解，再滤除分解信号中的高频噪音；步骤2：分别构建每个IMF分量的Hankel矩阵；步骤3：基于每个IMF分量的Hankel矩阵分别进行奇异值分解得到主成分分析的得分向量并进行一次去噪，以及对一次去噪后的各个分量信号和残余分量分别进行软阈值二次去噪；步骤4：将二次去噪后的分量信号和残余分量进行叠加得到去噪后的微震信号。本发明将奇异值分解与主成分分析相关联，以奇异值分解的信息作为主成分分析的得分向量，简化PCA计算过程，克服了单一列向量无法利用奇异值分解的得分向量进行去噪的缺陷。

Description

一种微震信号多尺度去噪方法、装置及可读存储介质
技术领域
本发明属于微震信号处理技术领域，具体涉及一种微震信号多尺度去噪方法、装置及可读存储介质。
背景技术
微震监测技术作为一种先进高效的地压活动监测手段，在矿山、水电、隧道等领域已得到广泛的应用。微震监测技术通过传感器采集岩爆、断层滑移等产生的地震波信号，再对地震波信号进行处理分析，从而得到事件发生的位置、震级大小、能量等信息，为矿山评估风险提供依据。然而，矿山环境复杂，传感器监测到的信号通常带有大量背景噪音，这会严重影响信号信噪比、能量计算以及降低P波初至拾取精度，进而影响震源定位和震源机制反演等后续工作。因此，开展微震信号降噪方法的研究具有重要意义。
现有技术中，利用傅里叶变换可得到信号主频信息，从而用滤波器筛选得到所需频段，实现信号去噪。傅里叶变换常用低通、带通、高通进行滤波，但其存在一个明显的缺点：滤掉的频率可能含有有用成分。因此，研究者提出了将小波分析应用于微震信号处理，其中，Tangna等结合小波分析和自适应阈值对微震信号进行降噪，提出了用于微震信号去噪的自适应阈值剪切波变换(ATST)方法。Mousavi等将同步压缩小波变换和自定义单通道阈值相结合来进行微震信号降噪，取得了较好的降噪效果。Wei等提出了一种经验小波变换和基于聚类的阈值自动噪声衰减算法，计算简便，降噪效果较好。此外，Zhang等基于小波变换和Akaike信息标准(AIC)选取器开发了一种自动P波到达检测和选取算法，能够降低噪音对微震震相拾取的影响。小波分析方法是一种窗口大小固定但其形状可改变，时间窗和频率窗都可改变的时频局部化分析方法。其可将信号分解成不同的频道和频率成分，而且可以通过伸缩、平移聚焦到信号的任一细节加以分析。但是小波变换在信号高频段频率分辨率较差，而在低频段其时间分辨率又较差，并且其需选取合适的小波基才能达到较好的降噪效果。
进一步地，Huang等提出了经验模态分解方法(EMD)，可将时间信号分解为具有不同特征尺度的数据序列(IMF)，在微震信号降噪中的应用非常广泛，在一定程度上克服了小波变换的不足。而针对EMD分解过程中会出现模态混叠现象，为了降低模态混叠的影响，Wu和Huang提出了基于噪声辅助的集合经验模态分解(Ensemble empirical modedecomposition，EEMD)。基于EMD和EEMD的优势，如何充分利用该算法提高微震去噪效果或者提供另一种技术思路来实现微震信号去噪是可以进一步研究与探讨的。
PCA(Principal component analysis主成分分析)是一种分析、简化数据集的技术，其能保持数据集中对方差贡献最大的特征。现有技术中已存在将PCA应用于信号去噪中，然而利用PCA处理单一分量信号时，常规是选择协方差矩阵特征值作为信息筛选，然而计算繁琐，针对单一分量信号如何实现另一种信息筛选是值得探讨的。
发明内容
本发明的目的是将提供一种全新的技术手段来实现微震信号去噪，其中，先对微震信号进行EMD或EEMD分解，然后针对单一分量信号，将奇异值分解与主成分分析关联，以奇异值分解的信息作为主成分分析的得分向量，省去了常规主成分分析手段中均值、协方差计算等步骤，简化了计算过程；同时，为了实现单一分量信号的奇异值分解，本发明将单一分量信号构建成一个多维的Hankel矩阵，克服了单一分量信号(单一列向量)无法利用奇异值分解的得分向量去噪的缺陷。
一方面，本发明提供一种微震信号多尺度去噪方法，包括如下步骤：
步骤1：获取微震信号并进行EMD或EEMD分解，再滤除分解信号中的高频噪音，其中，EMD或EEMD分解后得到多个IMF分量和1个残余分量；
步骤2：分别构建每个IMF分量的Hankel矩阵，所述Hankel矩阵为多维矩阵；
步骤3：基于每个IMF分量的Hankel矩阵分别进行奇异值分解得到主成分分析的得分向量并进行一次去噪，以及对一次去噪后的各个分量和残余分量信号分别进行软阈值二次去噪；
步骤4：将二次去噪后的分量信号和残余分量信号进行叠加得到去噪后的微震信号。进一步优选，一个IMF分量的Hankel矩阵如下所示：
式中，Hi表示第i个IMF分量ci的Hankel矩阵，ci(N)表示IMF分量ci的第N个采样点，p、q为Hankel矩阵Hi的行数、列数；N为IMF分量ci的信号长度。
进一步优选，步骤3中基于第i个IMF分量ci的Hankel矩阵Hi进行一次去噪的过程如下：
首先，对Hankel矩阵Hi进行奇异值分解得到得分向量和负荷向量，分解如下：
式中，为第j个得分向量，/>为第j个负荷向量，/>Ti为得分矩阵，/>Pi为负荷矩阵，T为矩阵转置符号；R为实数，p、q为Hankel矩阵Hi的行数、列数；
然后，将得分向量按照长度从大到小排列，再选择得分向量长度与总得到向量长度的累积百分比超过阈值时确定的前k个负荷向量和前k个得分向量；
接着，利用所述前k个负荷向量和前k个得分向量进行重构得到Hankel重构矩阵；
最后，利用Hankel重构矩阵中第一行和最后一列确定IMF分量ci对应的一次去噪后的分量。
其中，优选阈值设定为85％，IMF分量ci对应的一次去噪后的分量c′i＝[c′i(1)c′i(2) … c′i(q) c′i(q+1) c′i(N)]。
进一步优选，所述Hankel重构矩阵如下所示：
式中，表示Hankel矩阵Hi对应的Hankel重构矩阵，/>表示前k个得分向量构成的得分矩阵，/>表示前k个负荷向量构成的负荷矩阵。
进一步优选，对一次去噪后的分量信号和残余分量进行软阈值二次去噪的过程如下：
式中，Thi表示IMF分量ci对应的软阈值二次去噪后的分量信号，c′i表示IMF分量ci对应的一次去噪后的分量信号，残余分量作为cn+1，sgn(·)为符号函数；T为阀值，并满足：σ为分量信号c′i方差，N′为分量信号c′i长度。
进一步优选，步骤1中基于方差贡献率滤除分解信号中的高频噪音，过程为：计算出每个IMF分量的方差贡献率，再滤除方差贡献率低于预设值的IMF分量；
其中，被滤除的IMF分量为所述高频噪音，所述方差贡献率的计算公式如下：
其中，将将EMD或EEMD分解后得到的残余分量作为cn+1，v(i)表示第i个IMF分量ci的方差贡献率，N为IMF分量ci的信号长度，ci(j)表示IMF分量ci的第j个采样点信号，n为将EMD或EEMD分解后得到IMF分量的总个数。
进一步优选，步骤1中的所述微震信号是从监测信号中截取得到，过程如下：
首先，利用如下公式计算出监测信号的长短时窗均值比时间序列；
式中，λ(k)表示长短时窗均值比时间序列中第k点的值，s(t)表示监测信号中第t个采样点信号，W1为短时窗的长度，W2为长时窗的长度；
然后，基于长短时窗均值比时间序列获取微震信号P波初至，其中，将长短时窗均值比大于Th时为微震信号P波初至，2≤Th≤3；
最后，截取微震信号P波初至前m1个点和后m2个点范围对传感器监测到的监测信号进行截取得到微震信号。
进一步优选，短时窗的长度W1的取值范围为[40,60]，长时窗长度W2的取值范围为[150,250]，m1为1000，m2为3000。
第二方面，本发明还提供一种基于上述微震信号多尺度去噪方法的去噪装置，其包括：
微震信号获取单元，用于获取微震信号；
分解单元，用于对微震信号进行将EMD或EEMD分解；
高频噪音去噪单元，用于滤除分解信号中的高频噪音；
Hankel矩阵构建单元，用于构建每个IMF分量的Hankel矩阵；
主成分和软阈值去噪单元，用于基于每个IMF分量的Hankel矩阵分别进行奇异值分解得到主成分分析的得分向量并进行一次去噪，以及用于对一次去噪后的各个分量信号和残余分量分别进行软阈值二次去噪；
重构单元，用于将二次去噪后的分量信号和残余分量信号进行叠加得到去噪后的微震信号。
第三方面，本发明还提供一种终端，其包括处理器和存储器，存储器内存储了计算机程序，处理器调用所述存储器存储的计算机程序以执行上述微震信号多尺度去噪方法的步骤。
第四方面，本发明还提供一种可读存储介质，存储了计算机程序指令，所述计算机程序指令被处理器调用执行上述微震信号多尺度去噪方法的步骤。
有益效果
本发明提供的一种微震信号多尺度去噪方法首先利用EEMD或EMD对微震信号进行分解，再针对每个单一分量信号结合奇异值分解与主成分分析，将奇异值分解得到的U∑作为主成分分析的得分向量，省去PCA分析的均值、协方差计算等繁琐步骤，且与常规的协方差矩阵特征值是完全不同的一种主要信息筛选方式。尤其是，本发明是分别针对单一分量信号进行奇异值分解，而单一的列向量进行奇异值分解时，分解后的奇异值仅有一个非零值，即仅有一个得分向量为非零向量，由此不能使用得分向量长度进行去噪。为此，本发明构建出多维Hankel矩阵来增加数据维度，继而得到多个非零奇异值，再以得分向量的长度筛选主要信息，实现微震信号去噪。
附图说明
图1是本发明实施例提供的所述方法的流程图；
图2是“Blocks”原始信号(a)和含噪信号(b)；
图3是图2中图(b)含噪信号的EEMD分解结果图；
图4是图3中c3构成的Hankel矩阵的得分向量与累积得分向量百分比图；
图5是第3个内禀模态函数c3去噪前后对比图；
图6是降噪后的“Blocks”信号图。
具体实施方式
如图1所示，本发明提供的一种微震信号多尺度去噪方法，其对于EEMD或EMD分解得到的高频模态分量直接去除，剩余模态分量则分别构建Hankel矩阵，再借助Hankel矩阵进行奇异值分解得到主成分分析中的得分向量，进而实现一次去噪，以及再利用软阀值进行二次去噪，在保持模态分量特征的同时进行去噪，最后将去噪后的模态分量相加得到去噪后的微震信号。下面将结合实施例对本发明做进一步的说明。
本发明实施例提供的一种微震信号多尺度去噪方法，其使用的是集合经验模态分解，具体包括如下步骤：
Step1：借助长短时窗均值比法从监测信号中截取出微震信号。其中，传感器监测到的监测信号s(t)(连续信号)大部分为噪音序列，本实施例中利用下述公式(1)计算出监测信号s(t)的长短时窗均值比(STA/LTA)时间序列：
式中，λ(k)表示长短时窗均值比时间序列中第k点的值，s(t)表示监测信号中第t个信号，W1为短时窗的长度，W2为长时窗的长度。
然后，设置长短时窗均值比大于Th(2≤Th≤3)时为微震信号P波初至，进而得到微震信号P波初至；
最后，取P波初至前1000个点和后3000个点范围对传感器监测到的信号时间序列进行截取，得到微震信号x(t)。
其他可行的实施例中，截取的范围在满足需求的前提下可以进行适当调整，本发明对此不进行具体的限定。
Step2：利用集合经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)分解微震信号，再利用方差贡献率去除高频噪音。
其中，借助集合经验模态分解具有高自适应的特点，将微震信号x(t)自适应分解为n个IMF分量ci和1个残余分量rn，i＝1,2,…,n。
然后，利用公式(2)计算出每个IMF分量的方差贡献率，并滤除方差贡献率低于0.01的IMF分量。
其中，残余分量rn记为了cn+1，v(i)表示第i个IMF分量ci的方差贡献率，N为IMF分量ci的信号长度(信号总采样点数)，ci(j)表示IMF分量ci的第j个采样点信号。其他可行的实施例中，在满足需求的前提下不限定于0.01。
Step3：构建剩余IMF分量中第一个IMF分量的Hankel矩阵。其中，普通技术而言，单个微震信号通常是一个行向量，难以开展矩阵分解、主成分分子等操作，本发明针对单个微震信号构建多维矩阵，即Hankel矩阵。
第i个IMF分量ci的Hankel矩阵如下所示：
式中，Hi表示第i个IMF分量ci的Hankel矩阵，ci(N)表示IMF分量ci的中第N个采样点，p、q为Hankel矩阵Hi的行数、列数。
Step4：基于Hankel矩阵进行奇异值分解得到得分向量，并进行一次去噪；其中，IMF分量ci通常包含噪音，及对应的Hankel矩阵Hi包含冗余信息，为此，本发明借助主成分分析筛选出主要特征，去除冗余信息。
原理如下：
Hi由奇异值可分解为Ui∈Rp×p，∑∈Rp×q，Vi∈Rq×q，记Ti＝Ui∑，Pi＝Vi，则Hi可由q个向量的外积之和，即：
式中，为得分向量，/>为负荷向量。/>为得分矩阵；为负荷矩阵，T为矩阵转置。
将得分向量按照其长度由大到小排列，如：/>选用得分向量长度与总得分向量长度的累计百分比85％为阀值，确定队列中对应的前k个负荷向量进行重构，得到重构后的Hankel重构矩阵为：
再以Hankel重构矩阵中的第一行和最后一列得到IMF分量ci对应地主成分分析后的去噪分量c′i。
Step5：对一次去噪后的分量信号进行软阈值二次去噪，软阈值计算公式如下：
式中，Thi表示IMF分量ci对应的软阈值二次去噪后的分量信号，残余分量作为cn+1，c′i表示IMF分量ci对应的一次去噪后的分量信号，sgn(·)为符号函数；T为阀值，并满足：σ为分量信号c′i方差，N′为分量信号c′i长度。
Step6：依次对剩余IMF分量分别重复上述Step3-Step5，对残余分量也要进行软阈值去噪。
Step7：将去噪后的分量信号和残余分量信号进行叠加得到去噪的微震信号。
综上所述，本实施例中使用EEMD对微震信号进行分解，再结合奇异值分解与主成分分析的关联性，以奇异值分解得到的U∑作为主成分分析的得分向量，可省去PCA分析的均值、协方差计算等步骤。此外，单一的列向量奇异值分解时，奇异值仅有一个非零值，即仅有一个得分向量为非零向量，由此不能使用得分向量长度进行去噪。为此，本发明构建了Hankel矩阵来增加数据维度，进而得到多个非零奇异值，再以得分向量的长度筛选主要信息。由上可知，本发明以以奇异值分解得到的U∑作为主成分分析的得分向量筛选主要信息，其与常规的协方差矩阵特征值是完全不同的主要信息筛选方式，且本专利简化了PCA计算过程；同时，构建Hankel矩阵能够克服单一列向量无法使用得分向量去噪的缺陷，更好地适用于本发明中单一分量信号的处理。
在一些实施例中，本发明还提供一种微震信号多尺度去噪装置，包括通信连接的微震信号获取单元、分解单元、高频噪音去噪单元、Hankel矩阵构建单元、主成分和软阈值去噪单元以及重构单元。
其中，微震信号获取单元，用于获取微震信号；
分解单元，用于对微震信号进行集合经验模态分解或经验模态分解；
高频噪音去噪单元，用于滤除分解信号中的高频噪音；
Hankel矩阵构建单元，用于构建每个IMF分量的Hankel矩阵；
主成分和软阈值去噪单元，用于基于每个IMF分量的Hankel矩阵分别进行奇异值分解得到主成分分析的得分向量并进行一次去噪，以及用于对一次去噪后的各个分量信号和残余分量分别进行软阈值二次去噪；
重构单元，用于将二次去噪后的分量信号和残余分量信号进行叠加得到去噪后的微震信号。
应当理解，上述单元模块的具体实现过程参照方法内容，本发明在此不进行具体的赘述，且上述功能模块单元的划分仅仅是一种逻辑功能的划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个装置，或一些特征可以忽略，或不执行。同时，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
在一些实施例中，本发明还提供一种终端，其包括处理器和存储器，存储器内存储了计算机程序，处理器调用所述存储器存储的计算机程序以执行上述微震信号多尺度去噪方法的步骤。
在一些实施例中，本发明还提供一种可读存储介质，其存储了计算机程序指令，所述计算机程序指令被处理器调用执行上述微震信号多尺度去噪方法的步骤。
应当理解，在本发明实施例中，所称处理器可以是中央处理单元(CentralProcessing Unit，CPU)，该处理器还可以是其他通用处理器、数字信号处理器(DigitalSignal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field-Programmable GateArray，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据。存储器的一部分还可以包括非易失性随机存取存储器。例如，存储器还可以存储设备类型的信息。
所述可读存储介质为计算机可读存储介质，其可以是前述任一实施例所述的控制器的内部存储单元，例如控制器的硬盘或内存。所述可读存储介质也可以是所述控制器的外部存储设备，例如所述控制器上配备的插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(Secure Digital,SD)卡，闪存卡(Flash Card)等。进一步地，所述可读存储介质还可以既包括所述控制器的内部存储单元也包括外部存储设备。所述可读存储介质用于存储所述计算机程序以及所述控制器所需的其他程序和数据。所述可读存储介质还可以用于暂时地存储已经输出或者将要输出的数据。
基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的可读存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，RandomAccess Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
实施例：为了验证本发明所述方法的可行性，本发明进行了如图2-图6的实验。
图2中的(a)为常用的“Blocks”测试信号，取信号采样长度为1024个点，对其施加0.2*randn(1,1024)的噪音，即噪音满足0.2倍标准正态分布，加噪信号如图2中的(b)所示，由此可见，加噪后“Blocks”信号的信噪比较低。
图3是图2中(b)图含噪信号的EEMD分解结果图。该含噪信号EEMD分解得到了9个IMF分量c1(t)～c9(t)和一个残余分量r9(t)。由图可知，随着EEMD的分解，IMF分量的整体频率降低，前2个IMF分量的高频噪音非常大。根据方差贡献率公式(2)计算各分量的方差贡献率，去除前几个小于0.01的分量，本测试去除了前2个IMF分量。
图4是图3中c3构成的Hankel矩阵的得分向量与累积得分向量百分比图。对c3构建Hankel矩阵H，对矩阵H进行奇异值分解，以奇异值分解得到的U∑作为主成分分析的得分向量，将得分向量的长度由大到小进行排序，并计算累计得分向量百分比(图5中(a))；设置累积得分向量百分比阈值为85％，将累积得分向量百分比大于85％后的得分向量置为0，由剩余载荷重构该分量。再使用软阀值，对重构分量进行去噪，得到图5中(b)。
图5是第3个内禀模态函数c3去噪前后的对比图。由图可知，c3分量在某些部位存在较大的噪音，而降噪后的c′3分量增强了这些部位的特征，且对其他段的振幅进行了压制，验证了该去噪方法的有效性。
图6是降噪后的“Blocks”信号图。对比图6与图2可知，该方法具有较好的降噪效果，去噪后的信号较好的保留了原信号特征。为定量评价降噪效果，计算得到降噪前和降噪后的信噪比分别为7.00和11.30，降噪效果明显。
需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。

Claims (10)

1.一种微震信号多尺度去噪方法，其特征在于：包括如下步骤：
步骤1：获取微震信号并进行EMD或EEMD分解，再滤除分解信号中的高频噪音，其中，EMD或EEMD分解后得到多个IMF分量和1个残余分量；
步骤2：分别构建每个IMF分量的Hankel矩阵，所述Hankel矩阵为多维矩阵；
步骤3：基于每个IMF分量的Hankel矩阵分别进行奇异值分解得到主成分分析的得分向量并进行一次去噪，以及对一次去噪后的各个分量和残余分量信号分别进行软阈值二次去噪；
步骤4：将二次去噪后的分量信号和残余分量信号进行叠加得到去噪后的微震信号。
2.根据权利要求1所述的方法，其特征在于：一个IMF分量的Hankel矩阵如下所示：
式中，Hi表示第i个IMF分量ci的Hankel矩阵，ci(N)表示IMF分量ci的第N个采样点，p、q为Hankel矩阵Hi的行数、列数；N为IMF分量ci的信号长度。
3.根据权利要求1所述的方法，其特征在于：步骤3中基于第i个IMF分量ci的Hankel矩阵Hi进行一次去噪的过程如下：
首先，对Hankel矩阵Hi进行奇异值分解得到得分向量和负荷向量，分解如下：
式中，为第j个得分向量，/>为第j个负荷向量，/>Ti为得分矩阵，/>Pi为负荷矩阵，T为矩阵转置符号；R为实数，p、q为Hankel矩阵Hi的行数、列数；
然后，将得分向量按照长度从大到小排列，再选择得分向量长度与总得到向量长度的累积百分比超过阈值时确定的前k个负荷向量和前k个得分向量；
接着，利用所述前k个负荷向量和前k个得分向量进行重构得到Hankel重构矩阵；
最后，利用Hankel重构矩阵中第一行和最后一列确定IMF分量ci对应的一次去噪后的分量c′i。
4.根据权利要求3所述的方法，其特征在于：所述Hankel重构矩阵如下所示：
式中，表示Hankel矩阵Hi对应的Hankel重构矩阵，/>表示前k个得分向量构成的得分矩阵，/>表示前k个负荷向量构成的负荷矩阵。
5.根据权利要求1所述的方法，其特征在于：对一次去噪后的分量信号和残余分量进行软阈值二次去噪的过程如下：
式中，Thi表示IMF分量ci对应的软阈值二次去噪后的分量信号，c′i表示IMF分量ci对应的一次去噪后的分量信号，残余分量作为cn+1，n为EMD或EEMD分解后得到IMF分量的总个数，sgn(·)为符号函数；T为阀值，并满足：σ为分量信号c′i方差，N'为分量信号c′i长度。
6.根据权利要求1所述的方法，其特征在于：步骤1中基于方差贡献率滤除分解信号中的高频噪音，过程为：计算出每个IMF分量的方差贡献率，再滤除方差贡献率低于预设值的IMF分量；
其中，被滤除的IMF分量为所述高频噪音，所述方差贡献率的计算公式如下：
其中，将EMD或EEMD分解后得到的残余分量作为cn+1，v(i)表示第i个IMF分量ci的方差贡献率，N为IMF分量ci的信号长度，ci(j)表示IMF分量ci的第j个采样点信号，n为EMD或EEMD分解后得到IMF分量的总个数。
7.根据权利要求1所述的方法，其特征在于：步骤1中的所述微震信号是从监测信号中截取得到，过程如下：
首先，利用如下公式计算出监测信号的长短时窗均值比时间序列；
式中，λ(k)表示长短时窗均值比时间序列中第k点的值，s(t)表示监测信号中第t个采样点信号，W1为短时窗的长度，W2为长时窗的长度；
然后，基于长短时窗均值比时间序列获取微震信号P波初至，其中，将长短时窗均值比大于Th时为微震信号P波初至，2≤Th≤3；
最后，截取微震信号P波初至前m1个点和后m2个点范围对传感器监测到的监测信号进行截取得到微震信号。
8.根据权利要求7所述的方法，其特征在于：短时窗的长度W1的取值范围为[40,60]，长时窗长度W2的取值范围为[150,250]，m1为1000，m2为3000。
9.一种基于权利要求1-8任一项所述方法的去噪装置，其特征在于：包括：
微震信号获取单元，用于获取微震信号；
分解单元，用于对微震信号进行EMD或EEMD分解；
高频噪音去噪单元，用于滤除分解信号中的高频噪音；
Hankel矩阵构建单元，用于构建每个IMF分量的Hankel矩阵；
主成分和软阈值去噪单元，用于基于每个IMF分量的Hankel矩阵分别进行奇异值分解得到主成分分析的得分向量并进行一次去噪，以及用于对一次去噪后的各个分量信号和残余分量分别进行软阈值二次去噪；
重构单元，用于将二次去噪后的分量信号和残余分量信号进行叠加得到去噪后的微震信号。
10.一种可读存储介质，其特征在于：存储了计算机程序指令，所述计算机程序指令被处理器调用执行权利要求1-8任一项所述方法的步骤。
CN202010760856.4A 2020-07-31 2020-07-31 一种微震信号多尺度去噪方法、装置及可读存储介质 Active CN111881858B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010760856.4A CN111881858B (zh) 2020-07-31 2020-07-31 一种微震信号多尺度去噪方法、装置及可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010760856.4A CN111881858B (zh) 2020-07-31 2020-07-31 一种微震信号多尺度去噪方法、装置及可读存储介质

Publications (2)

Publication Number Publication Date
CN111881858A CN111881858A (zh) 2020-11-03
CN111881858B true CN111881858B (zh) 2024-02-13

Family

ID=73204941

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010760856.4A Active CN111881858B (zh) 2020-07-31 2020-07-31 一种微震信号多尺度去噪方法、装置及可读存储介质

Country Status (1)

Country Link
CN (1) CN111881858B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112633225B (zh) * 2020-12-31 2023-07-18 矿冶科技集团有限公司 矿用微震信号滤波方法
CN112818830B (zh) * 2021-01-29 2022-11-25 上海跃磁生物科技有限公司 心磁信号降噪方法、系统、介质及装置
CN113392732B (zh) * 2021-05-31 2023-05-26 国网山东省电力公司电力科学研究院 一种局部放电超声信号抗干扰方法及系统
CN114200525B (zh) * 2021-12-10 2022-12-02 河北地质大学 一种自适应的多道奇异谱分析地震数据去噪方法
CN116451029B (zh) * 2023-06-15 2023-09-01 深圳瑞福来智能科技股份有限公司 一种除湿器工作状态预警方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011013610A1 (ja) * 2009-07-31 2011-02-03 富士フイルム株式会社 画像処理装置及び方法、データ処理装置及び方法、並びにプログラム及び記録媒体
EP2530607A1 (en) * 2011-06-01 2012-12-05 Thomson Licensing Method and computer program product for interactive changing distances between multidimensional data with structural properties projected on a display device
CN103824062A (zh) * 2014-03-06 2014-05-28 西安电子科技大学 基于非负矩阵分解的分部位人体运动识别方法
CN106338385A (zh) * 2016-08-25 2017-01-18 东南大学 一种基于奇异谱分解的旋转机械故障诊断方法
CN110826017A (zh) * 2019-09-25 2020-02-21 中国地质大学（武汉） 一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5506272B2 (ja) * 2009-07-31 2014-05-28 富士フイルム株式会社 画像処理装置及び方法、データ処理装置及び方法、並びにプログラム

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011013610A1 (ja) * 2009-07-31 2011-02-03 富士フイルム株式会社 画像処理装置及び方法、データ処理装置及び方法、並びにプログラム及び記録媒体
EP2530607A1 (en) * 2011-06-01 2012-12-05 Thomson Licensing Method and computer program product for interactive changing distances between multidimensional data with structural properties projected on a display device
CN103824062A (zh) * 2014-03-06 2014-05-28 西安电子科技大学 基于非负矩阵分解的分部位人体运动识别方法
CN106338385A (zh) * 2016-08-25 2017-01-18 东南大学 一种基于奇异谱分解的旋转机械故障诊断方法
CN110826017A (zh) * 2019-09-25 2020-02-21 中国地质大学（武汉） 一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
EEMD and Multiscale PCA-Based Signal Denoising Method and Its Application to Seismic P-Phase Arrival Picking;Kang Peng 等;Sensors 2021;第21卷(第16期);1-15 *
Noise suppression for gpr data based on svd of window-length-optimized hankel matrix.《Sensors》.2019,3807. *
一种适于存在极性反转的微震初至到时拾取方法;毕丽飞 等;《石油物探》;第59卷(第03期);344-355 *
基于 EMD 的胶合板损伤声发射信号特征提取及神经网络模式识别;徐锋 等;《振动与冲击;第31卷(第15期);30-35 *
基于 Hankel 矩阵 SVD 算法的去噪研究;崔少华 等;《实验室研究与探索》;第37卷(第02期);32-34 *
基于EEMD_Hankel_SVD 的矿山微震信号降噪方法;李伟 等;《煤炭学报;第43卷(第07期);1910-1917 *
李伟 等.基于EEMD_Hankel_SVD 的矿山微震信号降噪方法.《煤炭学报.2018,第43卷(第07期),1910-1917. *

Also Published As

Publication number Publication date
CN111881858A (zh) 2020-11-03

Similar Documents

Publication Publication Date Title
CN111881858B (zh) 一种微震信号多尺度去噪方法、装置及可读存储介质
Boashash et al. Time–frequency features for pattern recognition using high-resolution TFDs: A tutorial review
Hao et al. A joint framework for multivariate signal denoising using multivariate empirical mode decomposition
CN109164483B (zh) 多分量地震数据矢量去噪方法及多分量地震数据矢量去噪装置
Üstündağ et al. Denoising of weak ECG signals by using wavelet analysis and fuzzy thresholding
Traore et al. Structure analysis and denoising using singular spectrum analysis: application to acoustic emission signals from nuclear safety experiments
Kaloop et al. De-noising of GPS structural monitoring observation error using wavelet analysis
CN113392732B (zh) 一种局部放电超声信号抗干扰方法及系统
CN102184451B (zh) 一种基于多小波融合特征的神经元动作特征提取方法
CN111259269A (zh) 基于时间序列数据的异常点检测方法、设备和存储介质
CN113238193A (zh) 一种多分量联合重构的sar回波宽带干扰抑制方法
CN115797318A (zh) 一种光谱数据预处理方法、装置、计算机设备及存储介质
CN115563480A (zh) 基于峭度比系数筛选辛几何模态分解的齿轮故障辨识方法
Liu et al. The measurement and elimination of mode splitting: from the perspective of the partly ensemble empirical mode decomposition
Patil et al. Wavelet denoising with ICA for the segmentation of bio-acoustic sources in a noisy underwater environment
Quadri A review of noise cancellation techniques for cognitive radio
Yao et al. Railway axle box bearing fault identification using LCD-MPE and ELM-AdaBoost
Dawton et al. Proposal for a compressive measurement-based acoustic vehicle detection and identification system
Chen Rolling bearing fault diagnosis with compressed signals based on hybrid compressive sensing
US11936405B2 (en) Method for compressing digital signal data and signal compressor module
CN117251737B (zh) 闪电波形处理模型训练方法、分类方法、装置及电子设备
Neocleous et al. Identification of possible Δ14C anomalies since 14 ka BP: A computational intelligence approach
CN117033908A (zh) 一种基于多尺度小波的实时时频分析方法
Vaghmare BSS in Underdetermined Applications Using Modified Sparse Component Analysis
CN115877453A (zh) 一种三维地震信号提取方法、装置、终端及存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant