一种基于频域奇异值分解的矿山单通道微震信号降噪方法及应用

1.本发明涉及矿山微震及图像处理技术领域，尤其涉及一种基于频域奇异值分解的矿山 单通道微震信号降噪方法及应用。


背景技术：

2.采用常规的时域、频域分析方法很难对微震信号进行降噪处理。为此，国内外专家提 出了众多针对非平稳、非线性、瞬态变化信号的去噪方法。目前，常规滤波处理方法包括 带通滤波、emd滤波、小波包阈值滤波以及巴特沃斯滤波法、切比雪夫滤波法等方法。还 有以傅氏变换为基础的带通滤波、fk域滤波方法，小波阈值、emd等。徐宏斌等利用小 波变换对大尺度岩体结构下的微震监测信号进行去噪研究。kimiaefar r.等利用人工神经网 络和小波包分解对地震信号的随机噪声进行了压制处理。但小波包在自适应方面具有一定 的缺陷，在处理微震信号这类非平稳、瞬态信号时，会引起信号的畸变。
3.目前还存在一定的局限性：vmd是一种自适应信号分解方法，其能够克服处理信号 时发生的模态混叠现象，但该方法只是去除了含有噪音的模态分量，没有根据信号频谱特 征进行模态的选择，所以选择出来的模态并不一定包含原始信号的信号特征，另外在选取 有用模态分量后，只是对包含的工频噪音实现了降噪，使得这种降噪方法达到的降噪效果 并不是很好。小波包分解与人工神经网络相结合，对地震信号的随机噪声进行压制处理效 果良好，但小波包在自适应方面具有一定的缺陷，在处理微震信号这类非平稳、瞬态信号 时，会引起信号的畸变。针对弱能量、高频率、瞬态非平稳，且传播介质相对单一的矿山 微震信号，现有的微震信号处理方法难以提取出有效信号。


技术实现要素：

4.矿山微震信号降噪是为了将有效信号从背景干扰中提取出来，进而服务于波形分类识 别、到时拾取、定位计算以及属性特征的挖掘等环节。为了提高矿山微震信号的信噪比， 本发明提供了一种基于fft和svd的单通道微震信号降噪处理方法。
5.为解决上述技术问题，本发明的基于频域奇异值分解的矿山单通道微震信号降噪方法 包括如下步骤：
6.步骤s1：对微震信号的傅里叶变换
7.对微震信号x1进行傅里叶变换，将其转换到频域，得到与之对应的频域信号转 换公式为：
[0008][0009]
步骤s2：关键参数的确立
[0010]
确立单通道微震数据奇异值分解的相关参数：时间延迟量τ、重构阶数k以及hankel 矩阵长度n和维度m；
存在单通道微震时间序列x，该序列的自相关函数r可表述为：
[0028][0029]
其中，n为单通道微震记录的采样点数；r取最小值r
min
时所对应的延迟时间作为τ。
[0030]
进一步的，步骤s2中hankel矩阵长度n和维度m与τ的关系式为：(m-1)
×
τ+n＝n； 假设m＝n，则有：在对m、n的取值上，为了避免构建的hankel矩阵维数过高、 特征值过多，对计算时间和过程要求较高等缺点，可认为m、n近似相等。
[0031]
进一步的，步骤s2中引入奇异值能量差分谱用于确立奇异值分解重构svd阶数，奇 异值能量差分谱e描述了相邻奇异值所代表能量的变化情况，其计算公式可表述为：
[0032][0033]
其中，e(i)为第i个奇异值能量差分谱，i＝1,2,
…
,k-1，e为k-1个奇异值能量差 分谱组成的序列；λi为第i个奇异值，λ
max
和λ
min
分别为奇异值矩阵中的最大值和最小值。
[0034]
进一步的，步骤s6中的判断方法为：利用信号的信噪比snr、能量百分比esn、均 方根误差rmse以及信号平滑度sti分别对降噪前后的微震信号进行定量描述，其公式可 表述为：
[0035][0036][0037][0038][0039]
其中，sn和分别滤波前、后的微震信号，n为信号的采样点数。
[0040]
本发明还提供一种计算机程序产品，所述计算机程序产品包括计算机指令，所述计算 机指令在由处理器运行时使得计算机设备执行前述基于频域奇异值分解的矿山单通道微 震信号降噪方法。
[0041]
本发明还提供一种计算机可读存储介质，其上存储有计算机可执行指令，所述指令在 被处理器执行时用于实现前述基于频域奇异值分解的矿山单通道微震信号降噪方法。
[0042]
本发明还提供一种基于频域奇异值分解的矿山单通道微震信号降噪系统，包括：
[0043]
一个或多个处理器；以及
[0044]
一个或多个存储器，其中存储有计算机可执行程序，当由所述处理器执行所述计
算机 可执行程序时，执行前述基于频域奇异值分解的矿山单通道微震信号降噪方法。
[0045]
本发明的基于fft和svd的单通道微震信号降噪处理方法针对矿山微震系统各通道 信号相关性低、hankel矩阵维度低的特点，提出了单通道微震信号的hankel矩阵构建方 法，并提出了基于傅里叶变换与逆变换以及svd奇异值分解的fsvd频域降噪方法及其 处理流程。与传统方法相比，该方法具有不受阈值限制、兼顾了高低频带特性的特点，在 实现对微震信号内干扰成分的有效剔除的同时，提高了信号的信噪比、平滑度等，确保了 信号的保真度和分辨率。该方法对于矿山微震信号的快速分析处理，以及后续的初始到时 精确拾取、矿山微震信号特征挖掘与分析具有重要意义。
附图说明
[0046]
下面结合附图对本发明的具体实施方式做进一步阐明。
[0047]
图1(a)本发明的矿山微震信号svd频域降噪流程图
[0048]
图1(b)为本发明的降噪原理示意图；
[0049]
图2为本发明实施例中微震事件波形图；
[0050]
图3(a)为本发明实施例中序号1～5区间的奇异值的fsvd奇异值能量分布示意图；
[0051]
图3(b)为本发明实施例中序号1～5区间的奇异值的对应重构信号的频谱特性图；
[0052]
图4(a)为本发明实施例中序号6～10区间的奇异值的fsvd奇异值能量分布示意图；
[0053]
图4(b)为本发明实施例中序号6～10区间的奇异值的对应重构信号的频谱特性图；
[0054]
图5(a)为本发明实施例中序号11～15区间的奇异值的fsvd奇异值能量分布示意图；
[0055]
图5(b)为本发明实施例中序号11～15区间的奇异值的对应重构信号的频谱特性图；
[0056]
图6(a)为本发明实施例中序号16～20区间的奇异值的fsvd奇异值能量分布示意图；
[0057]
图6(b)为本发明实施例中序号16～20区间的奇异值的对应重构信号的频谱特性图；
[0058]
图7(a)为本发明实施例中序号21～25区间的奇异值的fsvd奇异值能量分布示意图；
[0059]
图7(b)为本发明实施例中序号21～25区间的奇异值的对应重构信号的频谱特性图；
[0060]
图8(a)为本发明实施例中高信噪比信号频域svd去噪结果分析示意图；
[0061]
图8(b)为本发明实施例中低信噪比信号频域svd去噪结果分析示意图。
具体实施方式
[0062]
本发明提出针对单个通道的奇异值分解(svd)去噪，利用该方法对矿山微震信号进 行分解降噪，主要是利用了单个通道内微震信号的周期性。利用svd对微震信号进行分解， 可以将该信号按能量大小划分为若干个本征值，并一定程度拓宽了该信号本征值的有
效频 宽，根据能量的分布对本征值进行相应的频率补偿，同时剔除以噪声为主的本征值，进而 对信号进行重构，这是该滤波方法的基本思想。重构信号的信噪比与分辨率与原始信号相 比，得到较大幅度的提升。
[0063]
要实现对单通道微震信号的奇异值分解，首先需要对微震信号进行“定长”划分。假 设单通道微震信号表述为x＝[x1,x2,x3,
…
,xn]，总采样点数为n。为便于svd分析，将单 通道信号划分为m维(道)，每一维(道)为采样点数为n的一列数据，后一列数据为前 一列数据延迟τ后的等长序列，由此构建的分解矩阵dm可表述为：
[0064][0065]
其中，m为矩阵重构的维数，τ为时间延迟量，n为每一维中信号的采样点数。
[0066]
此时，微震信号可以看作是由未受噪声干扰子空间d和噪声子空间n组成，其hankel 矩阵可表述为：
[0067][0068]
由此可以看出，u、v只是对原矩阵的旋转，而s则确立了矩阵的线性(压缩)程度， 线性程度越好则d约越逼近d，而利用svd进行降噪处理，实质上就是寻求对未受干扰 子空间d的最佳逼近，逼近的效果越好，则去噪效果越佳。
[0069]
利用上述分析，对分解矩阵dm进行奇异值分解svd后，奇异值由三部分组成，可表 述为：sw代表强干扰成分对应的奇异值，sn对应随机干扰对应的奇异 值，sd则是有效信号的奇异值。因此，利用svd对单通道微震信号去噪处理，其过程就 是保留sd对应的有效信号，而将sw和sn对应的干扰信号奇异值置为0，再进行svd反 变换得到去噪后的微震信号。对干扰信号奇异值置零的过程实际上是对矩阵dm进行“压 缩”的过程，如所示，
①
对应的是sw和sn部分，
②
则对应sd部分。
[0070]
本实施例的基于频域奇异值分解的矿山单通道微震信号降噪方法包括如下步骤：
[0071]
步骤s1：对微震信号的傅里叶变换
[0072]
对微震信号x1进行傅里叶变换，将其转换到频域，得到与之对应的频域信号转 换公式为：
[0073][0074]
步骤s2：关键参数的确立
[0075]
确立单通道微震数据奇异值分解的相关参数：时间延迟量τ、重构阶数k以及hankel 矩阵长度n和维度m；
[0076]
步骤s3：奇异值变换svd
[0077]
将二维微震信号进行等长度划分，构建hankel矩阵d，并对矩阵进行奇异值分解：假 设二维地震剖面为p，道数为m，每道的采样点数为n，对m
×
n矩阵p进行奇异值分解， 可得到如下关系式：
[0078][0079]
其中，u和v分别表示左右奇异阵(正交矩阵)，其中u∈rm×m、v∈rn×n；p的秩 为k(k＝min(m,n))，一般地m远小于n；s是由pp
t
或p
t
p的特征值σ按递减顺序组建的对 角矩阵；r是矩阵p的秩，奇异值个数k与矩阵的秩r相等。
[0080]
对角矩阵s＝diag(σ1,σ2,
…
,σk)为矩阵r的特征值，其中，σk为第k个特征值，其关系 式可表述为：
[0081][0082]
矩阵pp
t
或p
t
p的奇异值λk与特征值σk的关系可定义为其中λ1≥λ2≥
…ꢀ
≥λk≥0；信号在重构过程中，第k个特征量σk对整个信号的贡献与第k个奇异值λk是呈正 比的。因此，以λk或来表征该通道内地震信号的能量大小，则第j通道内信号的能量贡 献率cj可表述为：
[0083][0084]
可以看出，特征值或奇异值的分量越大，在整个地震信号中的贡献率越高。
[0085]
步骤s4：二维信号重构
[0086]
分析奇异值分布规律，并根据奇异值优选原则，确立合理的重构阶数k和奇异值序号； 利用svd反变换获得去噪后的单通道二维微震信号；
[0087]
一般地，由于奇异值是按从大到小顺序排列的，地震领域会选取最初的几个特征量(贡 献集中的部分)，来对原始信号进行描述，也就是保留对角矩阵s中的前几个有效奇异值， 将其他的奇异值置为0，然后利用奇异值分解的逆过程对信号进行重构。
[0088]
步骤s5：傅里叶逆变换
[0089]
利用傅里叶反变换将重构的频谱变换为期望的目标信号x2，逆变换公式如下：
[0090][0091]
步骤s6：判断去噪后信噪比是否满足要求，若不满足，则返回步骤s1。
[0092]
本实施例优选地，步骤s2中利用自相关函数获得svd相空间矩阵构建的延迟时间τ， 假设存在单通道微震时间序列x，该序列的自相关函数r可表述为：
[0093][0094]
其中，n为单通道微震记录的采样点数；r取最小值r
min
时所对应的延迟时间作为τ。
[0095]
本实施例优选地，步骤s2中hankel矩阵长度n和维度m与τ的关系式为： (m-1)
×
τ+n＝n；假设m＝n，则有：在对m、n的取值上，为了避免构建的hankel矩阵维数过高、特征值过多，对计算时间和过程要求较高等缺点，可认为m、n近 似相等。
[0096]
本实施例优选地，步骤s2中引入奇异值能量差分谱用于确立奇异值分解重构svd阶 数，奇异值能量差分谱e描述了相邻奇异值所代表能量的变化情况，其计算公式可表述为：
[0097][0098]
其中，e(i)为第i个奇异值能量差分谱，i＝1,2,
…
,k-1，e为k-1个奇异值能量差 分谱组成的序列；λi为第i个奇异值，λ
max
和λ
min
分别为奇异值矩阵中的最大值和最小值。
[0099]
进一步的，步骤s6中的判断方法为：利用信号的信噪比snr、能量百分比esn、均 方根误差rmse以及信号平滑度sti分别对降噪前后的微震信号进行定量描述，其公式可 表述为：
[0100][0101][0102][0103][0104]
其中，sn和分别滤波前、后的微震信号，n为信号的采样点数。
[0105]
为验证单通道信号svd去噪方法的有效性，本实施例以山东某矿现场实测的一次微 震事件为例(时间为2014-06-10 20:13:42)，对本发明所提出方法的去噪处理过程进行介 绍，微震事件波形图如图2所示。现场微震监测系统的相关参数如下：采样频率1000hz， 连续采集缓存(连续采集长度15min)，后续采用sta/lta进行事件的拾取与截取；传感 器选用速度型，频率特性为50～5khz，灵敏度为30v/g，采集的频率范围为0～1000hz。为 了完好地采集煤岩体破裂的微震信号，将微震传感器埋设于煤层内部(距孔口20～45m)。
[0106]
1、参数选取与数据处理
[0107]
以图2所示现场微震事件为例，下面就各关键参数的选取过程进行叙述。首先以单通 道9#为例，介绍个关键参数是如何获取，去噪后效果如何评判。
[0108]
(1)计算延迟时间
[0109]
针对单通道微震记录的特点，利用matlab的autocorr()函数求取信号的自相关系数， 延迟时间τ与自相关系数r的计算结果如下：
[0110][0111]
(2)构建hankel矩阵
[0112]
延迟时间计算结果τ＝9，按步骤s2，最终求解出m＝100，n＝109。由此，获得了hankel 矩阵构建的关键参数，利用这些参数可以实现对单通道的微震信号的svd分解；下一步 需要明确的是如何选取合理的svd重构阶数，即如何选择有效信号所对应的奇异值。
[0113]
(3)重构阶数的确立
[0114]
为了说明重构阶数选择的合理性，下面将按不同奇异值范围对信号进行重构，观察各 奇异值对原始微震信号的贡献程度。矿山微震信号的主频集中于50～200hz范围，噪声能量 比较强、分布宽。通过对该频谱进行奇异值分解，确立原始信号的能量谱主要集中在前20 个，因此，研究的对象为前20个奇异值序号所对应的信号成分。
[0115]
对构建好的hankel矩阵进行fsvd分解，并按一定规律对奇异值序列进行选择和信号 重构——分别选取序号1～5、6～10、11～15、16～20、21～25区间的奇异值，最后分别 得到fsvd奇异值能量分布，分别如图3(a)、4(a)、5(a)、6(a)和7(a)所示， 曲线为能量占比，线条对应奇异值编号，以及对应重构信号的频谱特性图，分别如图3(b)、 4(b)、5(b)、6(b)和7(b)所示，从图中可以看出，选择不同的序列(阶数)所重 构的信号特征不同，图3(a)、图3(b)以及图4(a)、图4(b)完整体现了波形信号 本身，而图5(a)、图5(b)；图6(a)、图6(b)以及图7(a)、图7(b)中噪声部 分已占据主导。图3(a)、图3(b)噪声得到有效抑制，但细节信息丢失较多，这主要是 由于选择的奇异值不完整所致；通过对图4(a)、图4(b)分析发现，该系列奇异值既有 有效成分又有干扰成分，属于一个过渡带；而图5(a)、图5(b)；图6(a)、图6(b) 以及图7(a)、图7(b)频谱特征体现出对波形有效成分的影响较小。因此，选择1～10 作为有效奇异值序列。从最终的去噪结果也能看出，图6(a)、图6(b)的去噪效果最佳， 主频成分得到有效保护，底噪压制干净，初时起跳点明显。
[0116]
2、频谱分析对比
[0117]
为了说明利用本发明提供的方法对矿山微震信号进行去噪处理的有效性，以5#、12# 通道为例(分别代表高、低信噪比信号)，对比该通道内资料去噪前、后的波形变化及频 谱变化规律，其结果如图8(a)和图8(b)所示。对比图8(a)和图8(b)中去噪前后 效果，可以看出，微震信号的噪声得到有效压制，时频分布范围更为集中，与此对应的是 明显被压制
的底噪。从时频图图8(a)中看，去噪前微震信号的频率广泛分布于0～400hz， 持续范围为200～800ms；而去噪后，600～800ms的频率成分被去除，320hz以上频段的 成分被削弱。图8(b)低信噪比信号的去噪效果更为明显，尤其是图8(b)中300～400hz 范围的频率完全去除，对应的波形噪声得到很好的抑制。
[0118]
由此说明，本实施例提供的方法从频率角度对微震信号进行去噪，在保留原有信号的 主要成分基础上，充分利用了矿山微震信号频率范围小的特点，对异常的频率成分进行了 抑制和去除，有效降低了背景噪声等干扰成分的影响，得到了更高的信噪比。
[0119]
本实施例还提供一种计算机程序产品，所述计算机程序产品包括计算机指令，所述计 算机指令在由处理器运行时使得计算机设备执行前述基于频域奇异值分解的矿山单通道 微震信号降噪方法。
[0120]
本实施例还提供一种计算机可读存储介质，其上存储有计算机可执行指令，所述指令 在被处理器执行时用于实现前述基于频域奇异值分解的矿山单通道微震信号降噪方法。
[0121]
本实施例还提供一种基于频域奇异值分解的矿山单通道微震信号降噪系统，包括：
[0122]
一个或多个处理器；以及
[0123]
一个或多个存储器，其中存储有计算机可执行程序，当由所述处理器执行所述计算机 可执行程序时，执行前述基于频域奇异值分解的矿山单通道微震信号降噪方法。
[0124]
本发明的基于频谱奇异值分解的单通道微震波形去噪方法实质是将微震信号转化到 频域，借助svd奇异值分解从频域角度对信号特征进行增强并压制干扰噪声。该方法避 免了有效信号损失严重等缺陷，有效实现了波长分离和去噪处理。该方法的提出为后续微 震资料的解译分析提供了保障。本发明通过对矿山单通道波形特征的分析，研究了单通道 波形的分解矩阵dm构建方法，提出了基于傅里叶变换和奇异值变换的fsvd降噪方法， 进而建立了该方法降噪流程及关键参数(延迟时间τ、重构阶数k、hankel矩阵)确立方 法。本发明完善了单通道矿山微震信号频谱去噪的流程及关键参数确立。通过对比分析， 本发明所提供的方法相对于多种常规单一去噪方法，能够更好地保留原始信号的信息(高 相似性、高能量百分比)，并能得到更高的信噪比，有效降低信号中的噪声成分，这对微 震信号的到时拾取、特征提取与分析具有实际意义。
[0125]
在以上的描述中阐述了很多具体细节以便于充分理解本发明。但是以上描述仅是本发 明的较佳实施例而已，本发明能够以很多不同于在此描述的其它方式来实施，因此本发明 不受上面公开的具体实施的限制。同时任何熟悉本领域技术人员在不脱离本发明技术方案 范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案做出许多可能的变动 和修饰，或修改为等同变化的等效实施例。凡是未脱离本发明技术方案的内容，依据本发 明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术 方案保护的范围内。