2018, Vol. 39
矿山微震技术常用来检测矿山岩体受力变形, 以及进行微震源定位.在矿山微震信号的接收过程中, 由于环境和仪器本身的影响, 常伴有强烈的随机噪声, 对于后期的数据处理和分析带来很多不便[1].自从1982年小波变换首次被引进地球物理研究领域之后, 小波变换在矿山微震信号去噪方面取得很多成果, 和以往研究方法相比, 去噪效果有了明显提高.
传统小波去噪利用的是全局阈值理论, 该阈值只是简单地根据信号的长度设定, 因此, 其设定条件的单一性必然无法满足信号的复杂特点.Donoho和Johnstone对全局阈值进行了改进, 根据信号在小波分解的一层系数取一中值, 与全局阈值做乘积运算, 这种方法考虑到了小波变换对信号去噪的实质是根据小波分解系数去噪这一特点, 但是对不同层小波分解系数的区别以及矿山微震信号在不同层的特点没有考虑进去.后人又对Donoho等的阈值设定方法作了进一步的完善, 将系数中值的取定与不同层小波分解系数联系起来, 但仍没考虑矿山微震信号在不同层的特点.
本文首先简要阐述了小波变换固定阈值及现存的改进阈值微震去噪方法, 并在此基础上进行改进, 提出了一种新的自适应阈值方法, 使其能够根据不同的分解层给出不同的阈值, 并全面地分析矿山微震信号在不同层的含噪特点.其次, 利用理论和实际数据进行测试, 同时, 与其他两种现存的阈值去噪方法对比, 分别计算信噪比、峰值信噪比和均方差.最后, 对所得结果、差剖面, 以及去噪评价标准进行分析, 验证所提方法的有效性与优越性.
1 小波变换及去噪原理 1.1 小波变换小波变换的实质是原信号与小波基函数的相似性.小波系数就是小波基函数与原信号相似的系数[2].设ψ(t)为平方可积函数, 且满足t∈L2(R), 对其傅里叶变换ψ(ω)满足如下条件:
|
(1) |
ψ(t)为一个基本小波函数或母函数[3].
将小波母函数ψ(t)进行伸缩和平移, 即
|
(2) |
式中:a为伸缩因子或尺度因子; b为平移因子.而经过伸缩平移后的ψa, b(t)称为连续小波基函数[4].
这里如果信号f(t)满足f(t)∈L2(R), 而将f(t)在式(2)小波基下展开, 这种展开即可称为连续小波变换, 简称CWT, 如式(3)所示:
|
(3) |
在连续小波变换的定义式中, 将其尺度参数a和平移参数b进行如下的离散化采样[5], 如式(4)和式(5)所示:
|
(4) |
|
(5) |
离散小波变换为
|
(6) |
设小波变换去噪模型为
|
(7) |
式中:x为含噪声信号; c为有效信号; e为噪声信号; σ为噪声强度.小波变换阈值去噪过程按照以下3个步骤进行[6]:
1) 对含噪信号进行小波变换.选择一个小波基并确定一个小波分解的层次N, 然后对信号x进行N层小波分解.
2) 对小波系数进行阈值处理.为保持信号的整体形状不变, 保留有效信号, 对分解后的每一层系数, 采用硬阈值、软阈值或其他阈值方法进行量化处理.
3) 进行小波逆变换, 对信号进行重构.
1.3 阈值理论本文选取硬阈值去噪函数.硬阈值处理是把信号小波变换系数的绝对值与阈值比较, 小于或等于阈值的小波系数变为零, 大于阈值的系数保持不变[7], 具有更好的保幅特点[8], 表达式为
|
(8) |
式中:T为阈值; s为小波分解系数.
2 阈值选取由于传统全局阈值不能灵活变化的缺陷, Donoho和Johnstone集中在均方误差准则上提出在全局阈值的基础上加上信号本身的特征, 对全局阈值进行改进[9], 将
|
(9) |
改进为
|
(10) |
式中, ε=(median(s))/0.674 5, median(s)为原信号s的中值, 反映信号本身的特征.
由于原信号本身的中值仅反映时间域特点, 没有考虑信号本身的频率特性, 后人对Donoho和Johnstone方法进行进一步的改进, 结合小波变换, 将原始信号进行小波分解求取中值, 见式(11):
|
(11) |
式(11)阈值改进方法是先将原信号用小波一层分解, 其系数储存方式为[c, l], c为系数空间集合, l为系数长度.结合系数和长度来取T2中的ε值.c1集合取自分解后的系数空间中, 其决定因素为分解后的低频系数个数m.若系数总个数为M个, 那么c1则取第m+1到第M个.
矿山微震信号的特点为低频部分有效信号所占比例大, 高频部分有效信号所占比例小, 因此对矿山微震信号进行小波分层过程中, 随着层数的增加, 高频信号越来越少, 低频信号越来越多.如果用上述阈值进行去噪的话, 可以发现该阈值并没有考虑到分层不同有效信号所占比例不同的这一特点.
因此, 本文在T3阈值算法基础上再次进行改进, 增加一个分层自适应因子r, 改进算法为
|
(12) |
式中:r=ln(j+j/J)，j为所要去噪的层级，J为去噪总层数；K为调谐常数.与传统方法相比, 本文提出方法具有以下优点:
1) 充分利用矿山微震信号的低频特性和随机噪声的高频特性;
2) 结合小波变换, 将原信号进行小波分解, 使矿山微震信号和随机噪声得到有效的分离, 层级越小所含随机噪声越多, 层级越大所含有效信号越多;
3) 自适应因子与小波分层数相关, 随着分层数的变大会使阈值减小, 可以最大限度地去除随机噪声, 且保护矿山微震信号的低频信息.
3 实际数据去噪利用某地实际矿山数据验证本文所提算法的准确性与优越性.图 1为所选实际矿山微震信号, 观察可知, 微震信号中包含了大量的随机噪声, 严重影响后续的监测分析与微震源定位.对所选微震信号进行傅里叶变换, 分析其频谱特征.由图 2可知, 其主频主要集中在10 Hz左右, 信号在高频部分分布较宽.由于矿山微震信号本身具有低频特性, 所以, 高频部分的表现基本为随机噪声, 需要进行压制或去除.
|
图 1 实际信号图 Fig.1 Actual signal data |
|
图 2 实际信号频谱 Fig.2 Actual signal spectrum |
利用小波变换将矿山微震信号进行分解, 所选小波变换尺度参数为4.图 3给出了矿山微震信号在不同尺度下的分解示意图, 每一个尺度下又分为低频部分和高频部分.通过分解图, 可以观察到信号所包含的不同频率成分.对比可知, 有效信号主要集中在低频系数上, 分布在高频系数上的有效信号较少.而随机噪声主要分布在高频系数上, 分布在低频系数上的随机噪声较少.
|
图 3 小波分解系数 Fig.3 Wavelet decomposition coefficients (a)—第一层分解低频系数；(b)—第二层分解低频系数；(c)—第三层分解低频系数；(d)—第四层分解低频系数；(e)—第一层分解高频系数；(f)—第二层分解高频系数；(g)—第三层分解高频系数；(h)—第四层分解高频系数. |
分别利用Donoho等提出的阈值计算方法式(9)和本文提出的阈值计算方法式(11), 通过硬阈值函数, 对所选矿山微震信号进行随机噪声去除处理.图 4为Donoho等提出的阈值计算方法所得结果, 随机噪声得到了较好的压制, 但由于所选阈值为全局阈值, 并没有考虑信号的分尺度特点以及矿山信号的低频特性, 部分有效信号损失.图 5为其去噪结果傅里叶变换的频谱, 观察可知, 高频部分基本无信号存在, 证明随机噪声得到了较好的压制; 然而低频区域部分信号衰减, 说明在随机噪声去除过程中, 损失了部分低频有效信号.图 6为其对应的差剖面, 显示了随机噪声的随机特性.图 7为本文提出的阈值计算方法去噪结果, 随机噪声同样得到较好的压制.由于提出的阈值计算方法考虑了信号在不同尺度的特点以及矿山微震信号的低频特性, 所以在随机噪声去除的同时, 最大限度地保留了有效信号, 与图 4相比, 许多细节得到了较好的保留.同样给出了去噪结果的傅里叶频谱, 如图 8所示, 高频部分基本无信号存在, 证明随机噪声得到了充分的压制, 而低频部分也得到了最大限度的保留, 说明基本没有损失有效信号.图 9为其对应的差剖面，与图 6相比, 随机性更强，更加符合随机噪声的特点.
|
图 4 小波去噪结果 Fig.4 Wavelet denoising result |
|
图 5 小波去噪频谱 Fig.5 Wavelet-denoising spectrum |
|
图 6 小波去噪差剖面 Fig.6 Wavelet-denoising differential profile |
|
图 7 小波自适应阈值去噪结果 Fig.7 Wavelet adaptive-denoising result |
|
图 8 小波自适应阈值去噪频谱 Fig.8 Wavelet adaptive-denoising spectrum |
|
图 9 小波自适应阈值去噪差剖面 Fig.9 Wavelet adaptive-denoising differential profile |
通过上述两种方法的去噪结果、频谱以及差剖面的对比, 证明本文提出的分层自适应阈值方法对随机噪声压制具有较强的能力, 同时还最大限度地保留了有效信息, 说明本文方法的有效性与准确性.
4 结论1) 本文通过对现存小波变换阈值去噪方法的研究, 以及矿山微震信号的分析, 充分地考虑了信号在小波域的分尺度特点以及矿山微震信号的低频特性, 提出了一种新的分层自适应阈值方法, 计算所得阈值能够随着小波变换的尺度变化, 更好地保护低频的矿山微震有效信号, 压制高频的随机噪声.
2) 通过实际矿山微震数据的验证, 以及与常规的全局阈值方法对比, 体现了本文所提出的分层自适应阈值方法的有效性与优越性, 对于矿山微震监测的数据分析与矿山微震源的准确定位, 具有重大的实际意义.
| [1] |
刘玉春. 小波变换在矿震信号滤波和识别的应用研究[D]. 阜新: 辽宁工程技术大学, 2008.
( Liu Yu-chun. Application of wavelet transform in mine seismic signal filtering and recognition[D]. Fuxin: Liaoning University of Engineering and Technology, 2008. http://cdmd.cnki.com.cn/Article/CDMD-10147-2010013097.htm ) |
| [2] |
Awala M A, Mostafaa S S, Ahmada M, et al.
An adaptive level dependent wavelet thresholding for ECG denoising[J]. Biocybernetics and Biomedical Engineering, 2014, 34(4): 238–249.
DOI:10.1016/j.bbe.2014.03.002 |
| [3] |
Han X C, Zhang P, Zhang J, et al. Application of wavelet transform in signal denoising[C]//2003 International Conference on Machine Learning and Cybernetics. Xi'an, 2003: 436-441.
|
| [4] |
Lang M, Guo H, Odegard J E, et al.
Noise reduction using an undecimated discrete wavelet transform[J]. IEEE Signal Processing Letters, 1996, 3(1): 10–12.
DOI:10.1109/97.475823 |
| [5] |
管亮, 冯新泸.
基于小波变换的信号消噪效果影响因素研究及其Matlab实践[J]. 自动化与仪器仪表, 2004(6): 43–46.
( Guan Liang, Feng Xin-lu. Study on influencing factors of signal denoising effect based on wavelet transform and its Matlab practice[J]. Automation and Instrumentation, 2004(6): 43–46. DOI:10.3969/j.issn.1001-9227.2004.06.016 ) |
| [6] |
Tomic M. Wavelet transforms with application in signal denoising[R]. Vienna: DAAAM International Vienna, 2008: 1-3.
|
| [7] |
Al-Ajlouni A M, Abo-Zahhad M, Ahmed J, et al.
ECG signal performance denoising assessment based on threshold tuning of dual-tree wavelet transform[J]. Journal of Medical Engineering, 2008, 16: 425–433.
|
| [8] |
Varady P. Wavelet-based adaptive denoising of phonocardiographic records[C]//Proceedings of the 23rd Annual International Conference of the IEEE. Istanbul, 2001: 1846-1849.
https://ieeexplore.ieee.org/document/1020582/ |
| [9] |
Alireza G, Mohammad A R.
Seismic coherent and random noise attenuation using the undecimated discrete wavelet transform method with WDGA technique[J]. Journal of Geophysics and Engineering, 2012(9): 619–631.
|