2. 中南大学地球科学与信息物理学院, 长沙 410012
2. School of Geosciences and Info-Physics, Central South University, Changsha 410012, China
0 引言
在矿集区，复杂的环境噪声和工业密集造成的人文噪声给大地电磁勘测带来很大的困扰[1]。如矿区作业机械的突然开关，通信中的电磁辐射，电力传输造成的电磁场，轨道、车辆等产生的噪声等等[2-3]。其特点是噪声能量大，频率集中在有限频带内，且这些噪声是多维的，并含有多次谐波，给大地电磁信号的响应函数估计带来很大偏差；随着经济的发展，这类噪声的影响也越来越严重。因此，在强噪声中提取出有用的大地电磁信号，对提高矿集区大地电磁方法的探测能力具有重要意义。
大地电磁工作者提出了很多诸如Robust处理、远参考技术、小波变换等行之有效的大地电磁信号去噪方法[4-9]。Robust方法能降低“飞点”的权，但无法消除输入端的噪声以及电磁相关噪声对数据的干扰[4-5]。远参考法将远参考点与实测点的资料进行相关处理，压制人文噪声对大地电磁资料的影响，但同步和测点的选择给作业带来了麻烦[6-7]。小波变换被广泛应用到大地电磁信号去噪，但分解层数和母小波的选择使得小波的应用受人为因素的影响较大[8-9]。1998年，美籍华人科学家Huang N. E.提出的经验模态分解(EMD)具有全自动的自适应分解能力；以膨胀、腐蚀、开启和闭合4种运算为基本特征的数学形态学方法可在强噪声中提取有用信号。经验模态分解和数学形态学的组合结合了两种算法的优点，也被广泛应用于故障诊断、光谱预处理和心电信号处理等[10-11]，也有学者对它们进行了改进，提出了基于改进的形态学滤波和集合经验模态分解的方法[11]。本文针对含强噪声的大地电磁(MT)信号，提出一种结合了EMD分解和数学形态学滤波的组合去噪方法，对矿集区大地电磁信号的时域信号进行滤波处理，研究其应用的方法与步骤，并用小波处理的结果作对比，评价方法的可行性。
1 方法原理 1.1 经验模态分解经验模态分解，即Huang变换，它将信号自适应地分解为满足以下两个定义条件的一组模态函数(intrinsic mode function，IMF)[12-13]：1) 对于一列数据，极值点和过零点数目必须相等，或至多相差一点；2) 在任意点，由局部极大点和极小点构成的两条包络线的平均值为零。
信号x(t)由EMD分解得到IMF的算法叙述繁琐，详尽步骤可参阅文献[12]，此处不再累述。经EMD分解后，原序列x(t)可表示为IMF分量和一个残余项的和：
(1) 式中：n为分解层数；ci(t)为各阶IMF；rn(t)为残余项。EMD是自适应地从高频到低频的分解过程，经过EMD分解后的各阶模态函数能够完全重构，几乎没有能量损失[12-13]。
1.2 数学形态学滤波数学形态学包括腐蚀、膨胀、开启和闭合4种基本运算。在实际应用中，开启和关闭操作经常结合形成形态学滤波算法[11]，文中的广义形态滤波器是基于广义开-闭和闭-开运算来的。假设存在一个一维信号f(n)，它的域是F={0, 1, …, N-1}；结构元素是g(n)，域是G={0, 1, …, M-1}。这里，N>M，则腐蚀运算的定义为
(2) 膨胀运算定义为
(3) 式中：Θ和⊕分别代表数学形态学中的腐蚀和膨胀算子。依照腐蚀和膨胀算子，开启和关闭操作可分别列式如下：
(4) 式中，∘和∙分别代表数学形态学中的开启和关闭操作。则广义的数学形态学闭-开和开-闭操作可表示为
(5) 最后，广义的数学形态学滤波可写成[14]
(6) 按照统计学上的3σ规则，如果噪声是零均值且其方差是σ2，那么它的值几乎都集中在区间[-3σ, 3σ]。采用3σ规则提取信号x(t)的阈值，可用来分离特征信号和噪声。对于一些受噪的MT信号，高频噪声和工频干扰往往是主要污染源，且基本呈高斯分布，所以阈值处理可以用来消除噪声而保留有用信号。如果噪声的方差保持不变，提取有用信号时，就可以采用一个统一的阈值。但对于一些实测的大地电磁信号，噪声方差可能会发生缓慢变化或突变。因此，有必要根据噪声方差的变化来确定时变阈值。把信号x(t)(t=1, 2, ..., T)分成u段，T为时间长度，每段有Q个采样点，则uQ=T。对于每段，采用绝对均值法来估计噪声标准差，忽略有用信号对它的影响，表示为[15]
(7) 式中：σi是第i段信号的噪声标准差；xi(j)是被分离的第i段信号；Median是计算均值的函数。此处，采用三次样条插值法插值噪声标准差σi。那么，原始时变噪声的标准差就可以表示为
(8) 式中，Spline为三次样条插值函数。最后，用阈值方法处理含噪信号x(t)(t=1, 2, ..., T)。如果信号的瞬时绝对值小于三倍方差，则判定为噪声，并被强置为0；否则判为有用信号予以保留。提取的有用信号为
(9) 图 1给出了组合方法去噪的流程图。首先，MT信号x(t)被分解为一系列的模态函数IMFs，包括高阶模态函数h(t)和低阶模态函数l(t)；再将高频部分与低频部分分开，分别进行去噪处理。对于低阶模态函数l(t)，根据公式(6)，先用数学形态学滤波方法对它进行去噪处理，得到滤波结果z(t)：
(10) 式中，k(t)(t=1, 2, ..., V, V≪T)是结构元素。这个结果保留了MT信号低阶模态函数的分类信息。用l(t)减去z(t)得到峰值信号a(t)：
(11)
|
| Ⅰ.Hilbert-Huang变换；Ⅱ.数学形态学滤波；Ⅲ.基于形态学信号的提出；Ⅳ.自适应阈值滤波；Ⅴ.低阶IMF消噪。 图 1 基于组合方法的MT信号的去噪流程图 Figure 1 De-noising diagram of MT signal based on combinational filter |
|
|
通常情况下，被分离出来的a(t)信号依然保留了一些有用信息，特别是一些低频信号保留了一些深层的地质结构信息，丢失这部分信息会大大减小MT的深层勘测能力。所以对a(t)进行二次信噪分离是非常必要的，以尽可能多地保留真实MT信号的信息。根据公式(9)，对a(t)进行基于3σ规则的自适应阈值滤波处理，获得去噪后的信号a′(t)，这就是滤波后的峰值信号。a′(t)和来自数学形态学滤波后的z(t)相加得到去噪后的低阶模态函数l′(t)。
对于高阶模态函数h(t)，采用平滑滤波来消除噪声。对于MT信号来说，基线漂移干扰主要存在于高阶模态函数。对高阶模态函数执行多点平滑滤波，基线漂移有望得到抑制，计算格式如下[16-17]:
(12) 式中，V≪T。
(13) 现在，我们就可以构造整个去噪后的MT。来自低阶模态函数部分的去噪结果l′(t)加上来自高阶模态函数部分的去噪结果h′(t)，其和就被视为受噪MT信号x(t)的去噪结果x′(t)。
(14) 实验信号y(t)为
(15) 式中：w(t)为仿真的有用信号；n(t)为噪声。为验证方法的有效性，笔者在仿真信号中加入了几种不同的噪声进行处理。
图 2a、b分别为加噪前后的仿真信号(信噪比为-10 dB)，噪声包含三角波干扰、工频干扰、高频噪声和基线漂移。图 3a是含噪信号的EMD分解图，图 3b是对应的被分解后各阶IMF分量的频谱图。从图 2、图 3中可以看出，信号自适应地从高频到低频一层层地被分解，每一层都有不同的幅度和频率值。从图 3还可以看出：高频噪声和工频干扰主要存在于低阶模态函数(主要在IMF1-IMF4)；基线漂移和工频干扰的高次谐波存在于高阶模态函数；三角波由于其丰富的频率成分，分布在多阶模态函数上。这些就为消噪处理带来了方便。
|
| 图 2 原始仿真信号(a)和加噪后的信号(b) Figure 2 Qriginal simulated signal (a) and its noised signal (b) |
|
|
|
| 图 3 含噪仿真MT信号的EMD分解(a)和对应各阶IMF分量的频谱(b) Figure 3 EMD decomposition of the noised simulated MT signal (a) and Fourier spectrum of each of IMF (b) |
|
|
按照前面提出的方法，对于低阶模态函数，先采用数学形态学滤波方法提取有用信号，然后再用自适应滤波方法进行二次分离；对于高阶模态部分，采用多点平滑进行去噪处理；最后两部分的处理结果相加得到去噪后的信号。图 4展示了去噪前后和原始信号3种信号的时频分布。从图 4a中可看出，加噪后信号的能量幅值增加了好几个数量级，且在50 Hz和低频带的频率上出现了突出的奇异值；这些均来自仿真加噪的三角波干扰、工频干扰和基线漂移的影响。而在去噪后的时频分布(图 4b)中，对比原始信号的时频图(图 4c)，这些噪声都得到了有效抑制。
|
| a.去噪前；b.去噪后；c.原始的。 图 4 去噪前后仿真信号Hilbert时频能量谱 Figure 4 Comparative of 3-D Hilbert spectrum of simulated signal before and after de-noising |
|
|
采用SNR(信噪比)、NSR(噪声抑制率)[18-19]和SDR(信号失真比)[20]对去噪效果进行评价，它们的定义如下：
(16) 
(17) 
(18) 式中：Pw和Pn分别代表信号和噪声的功率谱；y(ti)和(ti)是去噪前后的信号。SNR反映了信号与噪声的能量大小比；而NSR越大、SDR越小，则说明去噪效果越好[21-22]。
笔者采用了小波软阈值法与本文方法进行对比，以评估本文方法的优越性。评价参数的对比结果如表 1所示。从表 1中可以看出：小波软阈值去噪方法和本文提出的组合滤波方法都有较好的消噪效果，在不同的SNR条件下，SDR都不大，但组合滤波方法有更低的SDR；就较高信噪比SNR=10 dB时，NSR也较小，为0.170 3，SDR为0.020 1；随着信噪比的减小，当SNR=-10 dB时，NSR由0.170 3增加到0.792 7，尽管SDR也从0.020 1增大到0.040 9，但值仍旧很小。这些参数说明EMD+数学形态学的滤波方法有很好的去噪和信号细节保持能力。
| 方法 | SNR/dB | NRS | SDR |
| 小波软阈值方法 | 10 | 0.187 9 | 0.046 7 |
| 5 | 0.534 1 | 0.040 2 | |
| -5 | 0.613 2 | 0.067 3 | |
| -10 | 0.834 7 | 0.141 6 | |
| 本研究方法 | 10 | 0.170 3 | 0.020 1 |
| 5 | 0.421 5 | 0.029 6 | |
| -5 | 0.598 8 | 0.031 8 | |
| -10 | 0.792 7 | 0.040 9 |
在矿区，井下大功率牵引设备产生的强电磁干扰、开采时的震动等都是常见的三角波干扰源，甚至构成由多个三角波连续出现而形成的锯齿波状干扰。它们一般只影响磁场信号，在时间序列中，通常具有如下表达式：
(19) 式中：a(t)为信号函数；u(t)为单位函数。实测点位于正在开采的集矿区，矿区内井下有很多正在运行的直流电力牵引矿车，单机工作电流为84 A；此外，还有多台碎石机、风钻等大型设备，它们的开关、负荷改变，均可产生强幅度的三角波干扰。
图 5所示为实测MT数据。显然，在图 5中，磁道Hy，Hx在记录时间段内出现明显的三角状噪声，且相关性较好，噪声局部能量较强，将正常磁场信号完全湮灭，而电道信号(Ex, Ey)却无此类波形出现。单个噪声幅值为正常信号幅值的数倍到数十倍。Hy信号的统计参数：最大值为6 912 nT，最小值为-8 588 nT，均值为39.65 nT，方差为5.31×106，能量为4.35×1010 nT2；Hx信号的统计参数：最大值为7 951 nT，最小值为-7 808 nT，均值为44.50 nT，方差为3.36×106，能量为2.75×1010 nT2。
|
| 图 5 受三角波干扰的MT信号 Figure 5 MT signal disturbed by triangle wave |
|
|
将受干扰的信号进行EMD分解如图 6所示。用本文提出的组合方法对各阶IMF进行消噪处理，
|
| 图 6 受三角波干扰Hy分量(a)和Hx分量(b)的EMD分解 Figure 6 EMD decomposed of MT signal Hy (a) and Hx (b)disturbed by triangle wave |
|
|
最后重构得到消噪后的信号，图 7为消噪后的4个大地电磁信号的时域波形图。消噪后Hy信号的统计参数：最大值为6 205 nT，最小值为-6 980 nT，均值为-9.06 nT，方差为4.54×105，能量为3.72×109 nT2；Hx信号的统计参数：最大值为6 182 nT，最小值为-7 080 nT，均值为0.34 nT，方差为4.01×105，能量为6.03×109 nT2。消噪后锯齿形状没有了，消噪后信号的方差减小到近原来的10%，信号变得平稳。
|
| 图 7 三角波干扰抑制后的MT信号 Figure 7 De-noised MT signal disturbed by triangle wave |
|
|
用去噪前后的数据分别计算电阻率曲线和相位曲线。采用方差和均值来评估去噪前后两参数曲线的平稳性，用曲线相似性参数(NCC)来评价去噪前后两参数曲线的相似程度。NCC定义如下：
(20) 式中，e(n)和d(n)为两离散序列。NCC的值为-1~1：-1代表变换前后两数据波形反向；0代表两波形正交；1代表完全相同[6, 21]。
去噪前计算的xy方向视电阻率ρxy曲线的方差是2.21×104，均值是313.64，去噪后计算的曲线方差是1.96×104，均值是312.78，两曲线的相似性参数为0.905；去噪前计算的yx方向视电阻率ρyx曲线的方差是1.13×104，均值是201.58，去噪后计算的曲线方差是1.11×104，均值是202.23，两曲线的相似性参数为0.913。去噪前计算的TE模式下相位曲线Φxy的方差是41.3，均值是42.2，去噪后计算的曲线方差是32.3，均值是42.2，两曲线的相似性参数为0.917；去噪前计算的的TM模式下相位曲线Φyx的方差是147.83，均值是40.68，去噪后计算的曲线方差是123.34，均值是40.70，两曲线的相似性参数为0.925。
图 8显示了去噪前后视电阻率曲线和相位曲线的对比。由图形和数据的对比可知，由消噪前后计算的曲线相似性参数都较大，即整体趋势是一致的，但细节发生了明显的改善，去噪后的曲线方差减小到了去噪前的一半，曲线的突变点得到了有效抑制(如ρxy曲线的3 278 Hz频点电阻率由原来的983 Ω·m下降到743 Ω·m)。
|
| 图 8 去噪前后大地电磁视电阻率曲线和相位曲线对比 Figure 8 Comparative plot of magnetotelluric response curves obtained for a measured site data before and after de-noising |
|
|
1) 在矿集区，大地电磁场污染严重，导致大地电磁观测信号信噪比低下，为了从受噪的时间序列中得到可靠的MT响应参数，把有用的MT序列与噪声分开是必要的。
2) 先对MT数据序列进行EMD分解，然后再采用数学形态学滤波，自适应阈值滤波和多点平滑对分解得到的各阶IMF进行去噪处理。这种组合滤波方法充分结合了EMD的多尺度分解及其可重构特性和其他几种去噪方法的优点。
3) 仿真和实测数据处理结果表明，用本文提出的组合滤波方法消除噪声是有效的，再消噪的同时尽可能多地保留了有用信息，提高了大地电磁信号的信噪比。
| [1] | Szarka. Geophysical Aspects of Man-Made Electro-magnic Noise in the Earth Review[J]. Geophysics, 1988, 9: 287-318. |
| [2] |
汤井田, 李灏, 李晋, 等. Top-hat变换与庐枞矿集区大地电磁强干扰分离[J].
吉林大学学报(地球科学版), 2014, 44(1): 336-343.
Tang Jingtian, Li Hao, Li Jin, et al. Top-Hat Transformation and Magnetotclluric Sounding Data Strong Interference Separation of Luzong Ore Concentration Area[J]. Journal of Jilin University(Earth Science Edition), 2014, 44(1): 336-343. |
| [3] |
汤井田, 李晋, 肖晓, 等. 数学形态滤波与大地电磁噪声压制[J].
地球物理学报, 2012, 55(5): 1784-1793.
Tang Jingtian, Li Jin, Xiao Xiao, et al. Mathematical Morphology Filtering and Noise Suppression of Magnctotclluric Sounding Data[J]. Chinese Journal of Geophysics, 2012, 55(5): 1784-1793. DOI:10.6038/j.issn.0001-5733.2012.05.036 |
| [4] |
朱威, 范翠松, 姚大为, 等. 矿集区大地电磁噪声场源分析及噪声特点[J].
物探与化探, 2011, 35(5): 658-662.
Zhu Wei, Fan Cuisong, Yao Dawei, et al. Noise Source Analysis and Noise Characteristics Study of MT in an Ore Concentation Area[J]. Geophysical and Geochemical Exploration, 2011, 35(5): 658-662. |
| [5] |
蔡剑华, 汤井田, 王先春. 基于经验模态分解的大地电磁资料人文噪声处理[J].
中南大学学报(自然科学版), 2011, 42(6): 1786-1790.
Cai Jianhua, Tang Jingtian, Wang Xianchun. Human Noise Elimination for Magnetotelluric Data Based on Empirical Mode Decomposition[J]. Journal of Central South University (Science and Technology), 2011, 42(6): 1786-1790. |
| [6] | Cai J H, Tang J T, Hua X R. An Analysis Method for Magnetotelluric Data Based on the Hilbert-Huang Transform[J]. Exploration Geophysics, 2009, 40: 197-205. DOI:10.1071/EG08124 |
| [7] | Oliver R, Junge A, Dawes J K. New Equipment and Processing for Magnetotelluric Remote Reference Observations[J]. Geophys, 1998, 132: 535-548. DOI:10.1046/j.1365-246X.1998.00440.x |
| [8] | Suto N, Harada M, Izutsu J. Time Variation of the Electromagnetic Transfer Function of the Earth Estimated by Using Wavelet Transform[J]. Proc Jpn Acad: Ser B, 2006, 82: 175-180. DOI:10.2183/pjab.82.175 |
| [9] | Trad D O, Travassos J M. Magnetotelluric Data Ana-lysis: Robust Filter in Discrete Wavelet Analysis[C]//13th Workshop on Electromagnetic Induction in the Earth. Hokkaido: [s. n.], 1996: paper 7.13. |
| [10] |
宗永涛, 沈艳霞, 纪志成. 基于改进的形态学滤波和EEMD方法的滚动轴承故障诊断[J].
江南大学学报(自然科学版), 2015, 14(5): 533-537.
Zong Yongtao, Shen Yanxia, Ji Zhicheng. Fault Diagnosis for Rolling Bearing Based on the Improved Morphological Filter and EEMD Method[J]. Journal of Jiangnan University (Natural Science Edition), 2015, 14(5): 533-537. |
| [11] |
蔡剑华, 胡惟文, 王先春. 基于组合滤波的鱼油二十碳五烯酸含量近红外光谱检测[J].
农业工程学报, 2016, 32(1): 312-317.
Cai Jianhua, Hu Weiwen, Wang Xianchun. Near-Infrared Spectrum Detection of Fish Oil Eicosapentaenoic Acid Content Based on Combinational Filtering[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(1): 312-317. DOI:10.11975/j.issn.1002-6819.2016.01.043 |
| [12] | Huang N E, Shen Z, Long S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis[J]. Proceedings of the Royal Society of London Series A, 1998, 454: 903-998. DOI:10.1098/rspa.1998.0193 |
| [13] |
蔡剑华, 熊锐. 基于频率切片小波变换的时频分析与MT信号去噪[J].
石油物探, 2016, 55(6): 904-912.
Cai Jianhua, Xiong Rui. Magnctotclluric Data De-Nosing Based on Time-Frequency Analysis of the Frequency Slice Wavelet Transform[J]. Geophysical Prospecting for Petrolcum, 2016, 55(6): 904-912. |
| [14] |
蔡剑华, 肖晓. 基于小波自适应阈值去噪的MT信号处理方法[J].
地球物理学进展, 2015, 30(6): 2433-2439.
Cai Jianhua, Xiao Xiao. Method of Processing Magnetotelluric Signal Based on the Adaptive Threshold Wavelet[J]. Progress in Geophysics, 2015, 30(6): 2433-2439. DOI:10.6038/pg20150601 |
| [15] | Donoho D L. De-Noising by Soft-Shresholding[J]. IEEE Transactions on Information Theory, 1995, 41: 613-627. DOI:10.1109/18.382009 |
| [16] |
黄进文, 王威廉. 平滑滤波对ECG的R波幅值影响研究[J].
云南大学学报, 2008, 30(6): 559-563.
Huang Jinwen, Wang Weilian. Research on the Effect of R-Peak in ECG Wave by Smoothing Filtering[J]. Journal of Yunnan University, 2008, 30(6): 559-563. |
| [17] |
刘鹏, 庞澄清, 吕颖. 改进平滑滤波算法在车辆组合导航系统中的仿真研究[J].
机械与电子, 2016, 34(1): 24-27.
Liu Peng, Pang Chengqing, Lü Ying. Simulation Study of Smoothing Filtering Algorithm in Vehicle Navigation System[J]. Machinery & Electronics, 2016, 34(1): 24-27. |
| [18] |
蔡剑华, 王先春. 基于经验模态分解的近红外光谱预处理方法[J].
光学学报, 2010, 3(1): 267-271.
Cai Jianhua, Wang Xianchun. Near Infrared Spectrum Pretreatment Based on Empirical Mode Decomposition[J]. Acta Optica Sinica, 2010, 3(1): 267-271. |
| [19] |
蔡剑华, 龚玉蓉, 王先春. 基于Hilbert-Huang变换的大地电磁测深数据处理[J].
石油地球物理勘探, 2009, 10(5): 617-620.
Cai Jianhua, Gong Yurong, Wang Xianchu. Magnetotelluric Data Processing Based on Hilbert-Huang Transform[J]. Oil Geophysical Prospecting, 2009, 10(5): 617-620. |
| [20] |
李夕兵, 张义平, 刘志祥, 等. 爆破震动信号的小波分析与HHT变换[J].
爆炸与冲击, 2005, 25(6): 528-535.
Li Xibing, Zhang Yiping, Liu Zhixiang, et al. Wavelet Analysis and Hilbert-Huang Transform of Blasting Vibration Signal[J]. Explosion and Shock Waves, 2005, 25(6): 528-535. DOI:10.11883/1001-1455(2005)06-0528-08 |
| [21] |
蔡剑华, 王先春. 基于Teager-Huang的滚动轴承故障诊断方法[J].
机械设计, 2012, 29(7): 54-58.
Cai Jianhua, Wang Xianchun. Rolling Bearing Fault Diagnosis Method Based on Teager Huang transform[J]. Journal of Machine Design, 2012, 29(7): 54-58. |
| [22] |
李继业, 武晓军, 刘长生. 前兆观测干扰信号频谱特征分析[J].
中国地震, 2015, 31(3): 574-583.
Li Jiye, Wu Xiaojun, Liu Changsheng. Analysis of Frequency Spectrum Characteristics of Disturbed Signal in Earthquake Precursory Observation[J]. Earthquake Research in China, 2015, 31(3): 574-583. |