2017, Vol. 60
2. 湖南师范大学物理与信息科学学院, 长沙 410081;
3. 西北综合勘察设计研究院, 西安 710003
2. Institute of Physics and Information Science, Hunan Normal University, Changsha 410081, China;
3. Northwest Research Institute of Engineering Investigations and Design, Xi'an 710003, China
大地电磁测深法(Magnetotelluric，简称MT法)诞生于20世纪50年代，是一种以天然交变电磁场为场源，通过测量地表相互正交的电场和磁场，获得地下电性结构信息的地球物理方法(Tikhonov，1950).与有源的电磁勘探方法相比，天然大地电磁场频带范围宽且本身信号非常微弱，野外观测到的大地电磁信号极易受到各种噪声的污染(魏文博，2002；Trad and Travassos, 2000).随着我国“深部探测技术与实验研究”(SinoProbe)专项的深入，不可避免地需要在长江中下游成矿带及典型矿集区开展大地电磁探测工作(吕庆田等, 2011, 2015；董树文等，2012).由于这些地区经济发达、人烟稠密、矿山密布、矿山开采的大功率直流电机车、高压电网、各种金属管网、广播电台、通讯电缆及信号发射塔等造成的电磁干扰严重污染了实际大地电磁信号，引起大地电磁阻抗估算偏差严重及测量获得的视电阻率值过度失真等状况，导致不能客观反映地下电性分布.这些都给地球物理勘探工作带来了巨大困扰，极大地影响了地下物性结构和电性分析的可解释性及采集数据本身的可靠性；而无论是研究深部地质构造还是寻找深部隐伏盲矿体，都不可能完全避开强干扰地区(肖晓等，2014).因此，如何运用现代信号处理技术剖析矿集区大地电磁信号与强干扰的特性，对矿集区中各种强干扰进行有效压制和精细分离是一项极具挑战性的工作，也一直是大地电磁测深领域最重要的研究前沿和热点之一.
大地电磁测深理论提出至今，噪声问题一直困扰着广大大地电磁研究者，如何消除大地电磁信号中的噪声干扰，提高大地电磁测深数据质量，是国内外长期瞩目并不断取得进展的研究课题.近三十年来，国内外大地电磁工作者都一直在不断寻求新方法排除干扰因素，以获得无偏的阻抗估计.所幸大地电磁处理技术与其它学科息息相关，尤其是随着信息技术的飞速发展和广泛应用，极大地推动了大地电磁数据处理的蓬勃发展：在傅里叶变换基础上提出的最小二乘法、自功率谱法和互功率谱法(Vozoff，1972；Kao and Rankin, 1977；Goubau et al., 1978)；从大地电磁阻抗张量估算出发提出的Robust法(Sutamo and Vozoff, 1991)；为解决噪声与信号相关性提出的远参考道法(Gamble et al., 1979)；基于功率谱和自相关本身存在的问题提出的高阶统计量法(王书明和王家映，2004)；为压制静态效应和工频干扰引入的小波变换和Hilbert-Huang变换法(宋守根等，1995；汤井田等，2008；Cai et al., 2009；Chen et al., 2012；Cai, 2013, 2014).同时，还有一些学者从不同角度提出的噪声压制方法，如广义S变换(景建恩等，2012)、方差比维纳滤波(Kapple，2012)、同步时间序列依赖(王辉等，2014)等，均在一定程度上对改善大地电磁数据质量起到了积极作用.
分析国内外相关文献可知，现有的各种去噪技术均具有一定的优势，同时也存在一定的局限性.由于矿集区面临复杂的噪声干扰环境，在矿集区开展大地电磁测深往往会导致采集到的大地电磁时间序列中包含众多明显不是天然大地电磁信号的异常波形，且噪声干扰的频谱通常分布在宽频带、甚至是全频带范围内，现有的频率域去噪方法对该类强噪声干扰几乎无能为力.针对上述一系列不利因素及实际情况，近年来围绕大地电磁信号和强干扰的时间域序列特征，运用数学形态滤波进行大地电磁噪声压制的研究已初露锋芒(Li et al., 2011；汤井田等, 2012a, 2012b, 2015；李晋等，2014b).由于数学形态学只涉及加、减等简单运算，且对大尺度的强噪声干扰具有明显的压制效果.因此，该方法适合研究和处理矿集区海量大地电磁数据.但是，该方法与现有其他去噪方法一样，仅是对采集到的大地电磁数据进行整体处理，并没有对大地电磁信号和强干扰进行信噪辨识，处理后的结果往往也损失了大地电磁低频段的缓变化信息，导致“过处理”(李晋等，2014a).鉴于大地电磁测深采集时间长，受噪声干扰程度也不完全一致，如果能寻求到鲁棒性的特征参数来定量辨识未受到噪声干扰或受干扰程度较小的数据段，并将这些较为“真实”的大地电磁信号保留下来不被“过处理”，仅对那些明显不是天然大地电磁信号的异常波形进行噪声压制，势必对提升大地电磁噪声压制精度、获取高分辨率的大地电磁有用信号，以及评价大地电磁测深曲线的客观真实性起到显著效果.
分形几何是分析复杂非线性信号的有效手段，为复杂信号提供了一种几何结构分析方法.它反映了复杂形体占有空间的有效性，是复杂形体不规则性的度量.该方法由于能定量描述信号形态的复杂度和规则度，已在描述非线性、非平稳信号的特征检测方面取得了一定的成果，并在故障识别、图像分割、目标检测、岩石电传导特性等领域得到了广泛应用(Xia et al., 2006；火元莲等，2013；Wei et al., 2015).数学形态谱是在数学形态学颗粒分析基础上提出来的一种有效分析复杂形体形状特征的方法.它基于数学形态学多尺度运算，并能刻画不同尺度下信号的微小变化，已成功应用于故障诊断、纹理特征描述、信号奇异变化等(Torre et al., 2013；张亢等，2013).本文结合数学形态学、分形维数和形态谱，针对矿集区典型强干扰类型和大地电磁微弱信号，尝试性地从信号处理的角度引入鲁棒性的特征参数--形态分形维数和形态膨胀谱熵，提出基于信噪辨识的矿集区大地电磁噪声压制方法.首先，对大地电磁信号和强干扰进行定量辨识；然后，有选择性地对明显不是天然大地电磁时间序列的异常波形进行噪声压制；最后，重构大地电磁有用信号，并对算法进行性能评价.仿真结果表明，与形态滤波整体处理相比，该方法能更精细地保留大地电磁信号低频段的缓变化信息，其结果更加真实地“逼近”于天然大地电磁信号的本质特征；在压制近源干扰的同时，该方法能有效避免形态滤波在大地电磁噪声压制中的“过处理”，卡尼亚电阻率曲线的整体形态更为光滑、连续，低频段的大地电磁测深数据质量得到了明显改善.
2 形态分形维数 2.1 分形维数分形可用来描述自然界中传统欧几里得几何不能描述的复杂无规律的几何对象，而分形维数是度量分形的重要参数，可有效反映系统的非线性特性和混沌吸引子的奇异程度，并能定量描述不同尺度下被分析对象的复杂度与规则度.分形维数有多种定义方式，大多数定义都是基于“尺度ε下的度量”，主要包括Hausdorff维、自相似维、盒维、容量维、信息维、关联维、Lyapunov维等(李兵等，2010).由于盒维数的计算直观易行，且复杂度较小，因而得到最为广泛的应用.以盒维数为例，信号的分形维数估计如下：
假设信号X为Rn中一个非空有界子集，N(ε)表示采用边长为ε的网格覆盖信号X所需的最少个数，则信号的分形维数D定义为
|
(1) |
实际计算时对lnN(ε)和ln (1/ε)进行最小二乘线性拟合，并求其直线的斜率即可得到对信号的分形盒维数估计.由此可知，盒维数的基本思想就是用尽可能小的网格划分Rn，通过计算覆盖信号X所需的最少盒子数N(ε)来对信号X进行度量.
2.2 形态分形维数估计由于盒维数采用了规则划分的网格，因而存在所需盒子数估计不准确、影响分形维数估计精度的问题.从式(1)可知，分形维数的估计实质为在不同尺度下对分形集，即不规则几何形状集合进行度量的过程，而数学形态学中的多尺度运算正好提供了一种在不同尺度下对被分析对象的几何形态进行度量的思路.因此，Maragos结合数学形态学提出了形态分形维数估计方法(Maragos and Sun, 1993).
以一维离散信号为例，结合分形维数和多尺度形态学的定义，形态分形维数DM的估计如下：
假设f(n)和g(m)分别为定义在F={0, 1, 2, …, N-1}和G={0, 1, 2, …, M-1 }上的离散函数，且N≥M.其中，f(n)表示待处理信号，g(m)表示结构元素，则ε尺度下的结构元素εg定义为
|
(2) |
式(2)中，符号⊕表示膨胀运算，εmax表示信号分析的最大尺度.从式(2)可知，ε尺度下的结构元素εg，即为单位结构元素g(m)自身膨胀ε次.
f(n)关于g(m)的ε尺度形态腐蚀和形态膨胀运算分别定义为
|
(3) |
|
(4) |
式(3)中，符号Θ表示腐蚀运算.从式(3)和式(4)可知，ε尺度下的形态腐蚀和形态膨胀运算，即为在每个尺度下将一维离散函数g(m)作为单位结构元素对待处理信号f(n)进行一次腐蚀和膨胀运算.
ε尺度下f(n)关于g(m)的覆盖面积Ag(ε)定义为
|
(5) |
当ε→0时，Maragos证明Ag(ε)满足如下条件：
|
(6) |
式(6)中，c为常数，DM为信号f(n)的形态分形维数估计，实际计算时采用最小二乘法对
由于形态分形维数不需要划分网格，仅采用一维形态学(膨胀和腐蚀运算)覆盖进行估计.因此，该维数不受信号幅值范围和旋转等情况的影响，具有更高的计算效率和计算精度.
3 形态谱熵 3.1 形态谱数学形态学的多尺度运算能对待处理信号中包含多种尺度下的不同形态特征成分进行描述，Maragos在多尺度形态学的基础上提出了形态谱的概念，目的是进一步直观反映多尺度形态分析的结果、刻画待处理信号不同尺度下的变化(Maragos，1989).
同样以上述一维离散信号为例，f(n)关于g(m)的开运算形态谱PSo(ε)和闭运算形态谱PSc(ε)分别定义如下：
|
(7) |
|
(8) |
式(7)和式(8)中，符号
将更宽泛的形态学算子(腐蚀运算和膨胀运算)来替代形态开运算和形态闭运算，可得到腐蚀运算形态谱PSe(ε)和膨胀运算形态谱PSd(ε)：
|
(9) |
|
(10) |
形态谱反映了待处理信号中具体存在哪种尺度的形态特征成分，若含有ε尺度的形态特征成分，则在形态谱ε尺度处会存在谱线；谱线值越大则说明待处理信号中包含与该尺度结构元素相匹配的形态特征成分越多.因此，形态谱能较好地描述信号在不同尺度下的精细变化，且能更好地体现信号的形态特征.在信息论中，熵表示信源平均不确定性的度量；熵值越大，不确定的信息则越多.借鉴熵的定义，给出形态谱熵的定义，即将不同尺度的形态谱整理为一个概率分布序列，通过求其熵值来表征形态谱的稀疏程度，从而揭示不同形态形状概率分布下的有序程度信息(郝如江等，2008).
为此，f(n)关于g(m)的形态谱熵H(f/g)定义如下：
|
(11) |
式(11)中，q(ε)=PS(ε)/A(f·Kg)表示信号f(n)中ε尺度形态特征成分出现的概率，即形态谱中s尺度处的谱线值与整个形态谱谱线值和的比值.如果信号的形态相似，则形态复杂度和形态谱熵也近似相等；如果形态特征不相似，信号的形态复杂度和形态谱熵则完全不同.因此，通过形态谱熵的大小可定量描述信号形态的复杂度，并对信号进行分类.
4 仿真实验分析 4.1 数据源为了剖析矿集区典型强干扰与天然大地电磁微弱信号之间的定量辨识关系，首先必须获取矿集区典型强干扰类型和相对“真实”的大地电磁“纯净”信号.
图 1所示为庐枞矿集区某测点的电道和磁道大地电磁测深数据时间域波形.
|
图 1 庐枞矿集区某测点的大地电磁测深数据时间域波形 Fig. 1 Time domain waveform of magnetotelluric sounding data in Luzong ore concentration area |
从图 1可知，由于该测点面临复杂的人文电磁干扰环境，导致采集到的电道Ex和Ey中出现大量类方波干扰，相应时刻的磁道Hy和Hx中则伴随出现大量类充放电三角波干扰和类脉冲干扰.这些异常干扰波形的幅值和能量远高于正常大地电磁有用信号，天然大地电磁信号几乎被完全湮没，结果往往会引起视电阻率曲线呈45°渐近线上升的近源效应(汤井田等，2012c).通过观测长江中下游矿集区(庐枞、铜陵等)海量大地电磁实测数据发现，类方波干扰、类充放电三角波干扰和类脉冲干扰经常出现在电道和磁道的时间域波形中.因此，文中将在这三类最常见的干扰类型作为矿集区典型强干扰类型进行分析和研究.
为了获取相对“真实”的大地电磁“纯净”信号，在人烟稀少、基本无电磁干扰的青海省海西蒙古族藏族自治州开展大地电磁测深工作.
图 2所示为青海地区某测点的电道和磁道大地电磁测深数据时间域波形.从图 2可知，电道和磁道的时间域波形中均未观测到矿集区中常见的非天然大地电磁场的异常波形.分析该测点周围的采集环境可知，测点周围地势平坦、未发现其他明显的噪声干扰源.
|
图 2 青海地区某测点的大地电磁测深数据时间域波形 Fig. 2 Time domain waveform of magnetotelluric sounding data in Qinghai area |
图 3所示为该测点几乎未受到噪声干扰时间段的卡尼亚电阻率-相位曲线.从图 3可知，该测点的视电阻率曲线光滑、连续，并没有出现45°上升的近源效应，且视电阻率值平稳，相位呈50°左右，这些均符合天然大地电磁场的正常逻辑.
|
图 3 青海地区某测点的卡尼亚电阻率-相位曲线 Fig. 3 Cagniard resistivity-phase curve of Qinghai area |
鉴于该地区几乎未受到明显的人文电磁干扰源影响，且采集到的大地电磁数据具备天然大地电磁场信号的本质特征.因此，文中将在该地区采集的大地电磁测深数据作为相对“真实”的大地电磁“纯净”信号来逼近天然大地电磁微弱信号进行分析.上述数据源的信号和典型强干扰均来自低频采样率的MT数据.接下来，本文将从形态分形维数和形态谱熵两方面研究矿集区典型强干扰类型和几乎未受到电磁干扰的大地电磁信号之间的定量辨识关系.
4.2 形态分形维数定量辨识将数据源中任选一组几乎未受到电磁干扰和含类方波、类充放电三角波及类脉冲干扰的大地电磁信号作为测试样本，其时域波形如图 4所示.
|
图 4 测试样本 Fig. 4 Test sample |
图 5所示为测试样本中含类方波干扰的大地电磁信号采用尺度ε为20和40的单位结构元素g={1, 1, 1 }进行腐蚀和膨胀运算的效果图.分析图 5可知，腐蚀运算是一种收缩变换，待处理信号的尖峰和边界不平滑的凸起部分被剔除、谷域被加宽，目标整体呈收缩状态.膨胀运算是一种扩张变换，待处理信号的谷值和边界不平滑的凹陷部分被填平、尖峰被扩展，目标整体呈扩张状态.腐蚀和膨胀之间的区域即为形态学覆盖的面积Ag(ε).
|
图 5 腐蚀膨胀运算效果图 Fig. 5 Effect of erosion operation and dilation operation |
图 6所示为测试样本中四种信号对ln (Ag(ε)/ε2)和ln (1/ε)的双对数图.其中，直线斜率即为信号的形态分形维数估计DM，通过最小二乘法对ln (Ag(ε)/ε2)和ln (1/ε)进行线性拟合得到，实验中ε均取20.
|
图 6 ln (Ag(ε)/ε2)和ln (1/ε)的双对数图 Fig. 6 Double logarithmic of ln (Ag(ε)/ε2) and ln (1/ε) |
分析图 6可知，几乎未受到电磁干扰的时间序列其形态分形维数最大(DM1=1.6906).根据分形维数的定义可以解释为，该段信号应最为复杂、最无规律可循.含类脉冲干扰的时间序列由于与几乎未受到电磁干扰的时间序列在波形形态上非常相似，除了突然受到尖脉冲干扰外，时域波形几乎看不出规律，因而其形态分形维数DM4(DM4=1.5374)与DM1最为接近.含类方波干扰和类充放电三角波干扰的时间序列由于具有比较明显的类周期特性，其形态分形维数DM2(DM2=1.4987)和DM3(DM3=1.3861)明显低于几乎未受到电磁干扰时间序列的形态分形维数，且DM3最小，从分形维数的含义可以认为，该段信号所包含的规律性和确定性最为强烈.
为进一步探讨矿集区典型强干扰类型和大地电磁微弱信号之间的定量辨识关系，任选含类方波干扰、类充放电三角波干扰和类脉冲干扰的大地电磁时间序列以及相对“真实”的大地电磁“纯净”信号中各30个样本，共计120个样本进行分析与研究，每个样本的采样时间均为1 min.
图 7所示为上述所选四类大地电磁时间序列样本的形态分形维数分布图.从图 7可知，几乎未受到电磁干扰的大地电磁时间序列其形态分形维数最大(DM1=1.6687~1.7635)，且与其他三类典型强干扰的形态分形维数区分明显；含类方波干扰(DM2=1.3427~1.5093)和类充放电三角波干扰(DM3=1.3801~1.4905)的大地电磁信号其形态分形维数较小，分形维数曲线出现混叠现象，导致噪声干扰类型很难区分.根据形态分形维数的定义可知，几乎未受到电磁干扰的信号应不存在明显的类周期性，且形态最为复杂、毫无规律，接近于随机噪声状态；含类方波干扰和类充放电三角波干扰的信号本身具有一定的规律性和确定性，导致形态分形维数较小；含类脉冲干扰的信号由于与几乎未受到电磁干扰的信号在形态上具有相似性，且又不具备另两类信号明显的类周期性，因而形态分形维数位于两者之间(DM4=1.5096~1.5913).该结论与测试样本(图 6)得到的结论一致.
|
图 7 形态分形维数分布图 Fig. 7 Distribution of morphological fractal dimension |
表 1所示为图 2中电道Ex、Ey分量的形态分形维数.该段数据是采样率为15 Hz的低频MT信号，每1 min为1帧，共分为10帧.
|
表 1 Ex、Ey信号的形态分形维数 Table 1 Morphological fractal dimension of Ex and Ey |
结合图 7分析表 1可知，该段低频MT信号的形态分形维数明显大于随机选取的其他三类典型强干扰形态分形维数，因此该参数可以较好地辨识低频MT信号和典型强干扰类型.
4.3 形态谱熵定量辨识图 8所示为采用形态膨胀、形态腐蚀、形态开和形态闭运算对测试样本中四种信号获取形态谱的分布情况.分析图 8可知，形态膨胀和形态腐蚀运算获得的形态谱在定性区分几乎未受到电磁干扰和受到典型强干扰的效果上明显优于形态开和形态闭运算获得的形态谱.形态开运算和形态闭运算获得的形态谱曲线跳变紊乱，在不同尺度处出现明显的谱线混叠；分析形态谱的定义可知，此时无法判断信号中不同特征成分的分布特点，从形态谱曲线的整体特征很难区分信号类型.形态膨胀运算获得的形态谱曲线其分布特征尤为明显，并未出现谱线混叠的现象，且具有较好的区分度.
|
图 8 形态谱分布情况 Fig. 8 Distribution of morphological spectrum |
由于形态膨胀运算获得的形态谱能达到分辨信号不同特征的目的，在此基础上对图 7所选样本进一步获取形态膨胀谱熵的分布情况如图 9所示.
|
图 9 形态膨胀谱熵分布图 Fig. 9 Distribution of morphological dilation spectrum entropy |
从图 9可知，由于四类信号采用形态膨胀运算得到的形态谱曲线复杂度不同，所获得的形态膨胀谱熵也不一样.几乎未受到电磁干扰的大地电磁信号其形态膨胀谱熵(E.U.1=1.8187~2.1527)最大，从熵的定义可以解释为，该类信号的复杂度最高、含不确定的信息最多.虽然含类脉冲干扰(E.U.4=0.8738~1.2344)的形态膨胀谱熵高于其他两类干扰，且含类方波干扰(E.U.2=0.3707~0.8838)和类充放电三角波干扰(E.U.3=0.5147~0.7937)的形态膨胀谱熵仍存在一定程度的重叠，但这三类受典型强干扰的大地电磁信号其形态膨胀谱熵明显低于几乎未受到电磁干扰的大地电磁信号的形态膨胀谱熵，且曲线具有明显的区分度.究其原因，形态膨胀谱熵能根据信号形状的变化准确地反映信号的形态特征信息，从整体上能定量描述形态膨胀形态谱曲线的复杂程度，并能对不同类型信号的特征分布信息进行区分和甄别.
表 2所示为图 2中电道Ex、Ey分量的形态膨胀谱熵，每1 min为1帧，共分为10帧.
|
|
表 2 Ex、Ey信号的形态膨胀谱熵 Table 2 Morphological dilation spectrum entropy of Ex and Ey |
结合图 9分析表 2可知，该段低频MT信号的形态膨胀谱熵明显大于随机选取的其他三类典型强干扰的形态膨胀谱熵.因此，该参数对低频MT信号和典型强干扰类型同样具有较好的噪声鲁棒性，可以用来定量辨识低频MT信号和大尺度典型强干扰类型.
5 实测资料处理 5.1 算法流程及步骤针对矿集区实测大地电磁数据，首先将测点的时间序列进行分帧，然后对每帧信号提取鲁棒性的特征参数进行定量辨识，最后对辨识出有强干扰的时间段进行噪声压制.以加拿大凤凰公司V5-2000采集到的大地电磁数据(TSL和TSH格式)为例，算法的具体步骤如下：
步骤1：将采样率为24 Hz的低频时间序列(Ex、Ey、Hx、Hy)进行分帧，每帧为1 min (1440个采样点)，帧与帧之间不重叠；
步骤2：计算电道(Ex和Ey)中每1帧信号的形态分形维数和形态膨胀谱熵；
步骤3：设置形态分形维数和形态膨胀谱熵的阈值参数Th1和Th2进行定量辨识，若小于其中任意一个阈值，则辨识为典型强干扰时间段；
步骤4：将辨识为强干扰时间段的电道(Ex和Ey)信号做形态滤波处理，磁道(Hx和Hy)信号遵循相关性的原则，即当电道Ex中某1帧辨识为典型强干扰时，磁道Hy相同时刻对应的帧同样做形态滤波处理；
步骤5：将形态滤波处理后的帧信号与该帧绝对值的算数平均值进行比较，进一步剔除残余的尖脉冲干扰；
步骤6：采样率为320 Hz和2560 Hz的高频时间序列则按24 Hz低频时间序列所对应的时刻辨识出噪声时段，并重复步骤4和步骤5进行相应的形态滤波处理和残余尖脉冲压制；
步骤7：重构大地电磁有用信号，获取卡尼亚电阻率-相位曲线，评价算法性能.
文中所选Th1经验值范围为1.45~1.55，Th2经验值范围为1.00~1.20.由于TSL格式(采样率为24 Hz)为连续采集的低频信号，一般仅需对TSL文件进行处理.
5.2 时间域滤波效果图 10和图 11分别所示为两段Ex和Hy大地电磁时间序列经形态滤波整体处理和本文所提方法处理的效果对比图.其中，Hy道大地电磁信号与相应的Ex道大地电磁信号在同一时刻采集，分别选取8640个采样点进行分析，共分为6帧.
|
图 10 一段Ex道和Hy道MT信号时间域处理效果对比图 Fig. 10 Effect comparison of Ex and Hy in the time domain |
|
图 11 一段Ex道和Hy道MT信号时间域处理效果对比图 Fig. 11 Effect comparison of Ex and Hy in the time domain |
从图 10可知，由于相关性原则，原始MT时间序列(Ex)中出现类方波干扰的时间段对应Hy道同样出现类充放电三角波干扰.经形态滤波整体处理后，Ex道中的类方波干扰和Hy道中的类充放电三角波干扰虽然得到了有效压制，但整个时间序列中几乎完全损失了曲线的缓变化信息.
表 3所示为图 10中原始MT时间序列(Ex)的形态分形维数和形态膨胀谱熵.分析表 3可知，第3帧信号的形态分形维数和形态膨胀谱熵明显低于阈值参数Th1和Th2的经验范围.因此，通过对大地电磁信号进行定量辨识，仅需将该帧含强干扰类型的大地电磁信号进行噪声压制.对比分析形态滤波整体处理，本文所提方法在压制上述强干扰的同时，能较好地保留几乎未受到电磁干扰或含干扰程度较少时间段的缓变化信息，重构的大地电磁信号中包含了更多丰富的低频细节成分.
|
|
表 3 每帧Ex信号的形态分形维数和形态膨胀谱熵 Table 3 Morphological fractal dimension and morphological dilation spectrum entropy of each frame signal of Ex |
从图 11可知，形态滤波整体处理方法由于缺少信噪辨识环节，仅通过对待处理信号直接进行形态滤波处理，处理后的Ex和Hy信号均在基线附近，本不含典型强干扰类型的微弱大地电磁信号的缓变化信息也被随之滤除.这些被形态滤波“过处理”的低频有用信号将严重影响卡尼亚电阻率-相位曲线低频段的数据质量.
表 4所示为图 11中原始MT时间序列(Ex)的形态分形维数和形态膨胀谱熵.分析表 4可知，Ex道的MT信号其形态分形维数和形态膨胀谱熵均大于阈值参数Th1和Th2的经验范围，由本文所提方法可知该时间段采集的大地电磁信号几乎未受到强干扰的影响.观测原始MT时间序列(Ex和Hy)也可发现，该段信号的时间域波形中并未出现类似于典型强干扰的异常波形.因此，采用本文所提方法对该段Ex道的MT信号进行信噪定量辨识，可以很好地保留Ex道及相同时刻Hy道MT信号中的低频缓变化信息.
|
|
表 4 每帧Ex信号的形态分形维数和形态膨胀谱熵 Table 4 Morphological fractal dimension and morphological dilation spectrum entropy of each frame signal of Ex |
图 12所示为庐枞矿集区B42742测点的卡尼亚电阻率-相位曲线.从图 12可知，原始数据的视电阻率曲线大约在100~0.05 Hz频段以近45°渐进线快速上升，且曲线较为光滑、几乎没有误差棒；在5~0.05 Hz频段，对应的相位呈0°或±180°，近15个频点的数据完全失真；在0.08 Hz时，视电阻率值竟然超过1000000 Ωm左右；在0.05~0.0005 Hz频段，视电阻率值突然下降至1000 Ωm左右，误差棒增大、曲线跳变剧烈，低频段的噪声鲁棒性明显变差.观测该测点电道和磁道的时间域序列可以发现，该测点包含大量类方波干扰、类充放电三角波干扰和类脉冲干扰的异常波形，能量幅值远高于正常大地电磁有用信号几个数量级.这些现象与CSAMT曲线的近区特性几乎完全一致，属于典型的近源效应.该类数据是矿集区最为“头疼”的数据，现有技术几乎失效，功率谱筛选也无能为力.
|
图 12 原始数据的卡尼亚电阻率-相位曲线 Fig. 12 Cagniard resistivity-phase curve of original data |
图 13和图 14分别所示为该测点经形态滤波整体处理和本文所提方法处理后的卡尼亚电阻率-相位曲线图.其中，Th1选为1.55，Th2选为1.20.从图 13可知，该测点经形态滤波整体处理后，视电阻率曲线近45°上升的现象有所缓解，5~0.05 Hz频段的相位曲线从最初的0°或±180°均有明显的抬升，视电阻率的最大值也下降了1个数量级.然而，视电阻率曲线的整体形态呈“倒U型”，在0.05~0.0005 Hz超低频段，曲线一直呈下掉趋势，且全时段的视电阻率-相位曲线几乎没有误差棒，处理后的数据仍无法进行功率谱筛选.究其原因，采用形态滤波整体处理虽能有效压制叠加在微弱大地电磁信号上的大尺度干扰和基线漂移，但由于没有信噪辨识等约束条件，导致原始数据中包含的一些低频缓变化成分也被直接剔除，其结果不能精确反映地下电性结构.
|
图 13 经形态滤波整体处理的卡尼亚电阻率-相位曲线 Fig. 13 Cagniard resistivity-phase curve of overall processing by morphological filtering |
|
图 14 本文所提方法处理的卡尼亚电阻率-相位曲线 Fig. 14 Cagniard resistivity-phase curve of the proposed method |
分析图 14可知，经本文所提方法先对大地电磁信号进行信噪辨识再有针对性地去除强干扰后，“倒U型”的视电阻率曲线形态完全消失，除个别频点外视电阻率-相位曲线的整体形态更加平稳，视电阻率值相对稳定；xy方向和yx方向的视电阻率最大值均下降了近2个数量级，近源干扰得到了有效压制.5~0.05 Hz频段的相位曲线显著抬升，对应的视电阻率曲线光滑、连续，近15个频点的数据质量得到了明显改善；0.05~0.0005 Hz超低频段曲线下掉的现象有明显缓解，且误差棒较原始数据有所减小.由此可知，本文所提方法保留了更为丰富、可靠的低频信息，整个低频段的数据质量得到了有效改善；与形态滤波整体处理相比，其结果更为合理、可靠、置信，且更能真实地反映测点本身所固有的大地电磁深部构造信息.
图 15所示为在图 14的基础上对部分误差棒较大的频点进行功率谱筛选得到的卡尼亚电阻率-相位曲线图.分析图 15可知，受近源干扰严重的实测点数据经本文所提方法处理后，再经简单的功率谱筛选即可获得更为光滑、连续的卡尼亚电阻率-相位曲线，其结果将有助于后续电磁法反演解释水平的提升.
|
图 15 经功率谱筛选得到的卡尼亚电阻率-相位曲线 Fig. 15 Cagniard resistivity-phase curve by power spectrum filtering |
针对形态滤波方法在大地电磁强干扰分离中存在“过处理”问题，本文在剖析矿集区典型强干扰类型和相对“真实”的大地电磁“纯净”信号定量辨识关系的基础上，根据其形态特征复杂度的不同将形态分形维数和形态膨胀谱熵应用于大地电磁信噪辨识，研究了基于信噪辨识的矿集区大地电磁噪声压制方法.通过仿真实验和实测资料处理，结果表明：形态分形维数和形态膨胀谱熵对受到类周期性噪声比较明显的矿集区大地电磁数据具有较好的信噪辨识区分度和噪声鲁棒性；经本文所提方法处理后，低频段的缓变化信息可以得到更为精细的保留，重构的大地电磁时间序列逐步“逼近”于天然大地电磁信号的本质特征；对于受到近源干扰的大地电磁数据，该方法获得的卡尼亚电阻率-相位曲线的整体形态更为光滑、连续，其结果更加真实地反映了测点本身所固有的大地电磁深部构造信息.本文所提方法为今后在矿集区建立鉴定大地电磁信噪辨识的评价准则，以及在矿集区开展高精细、高分辨噪声压制提供了一种新的解决思路，具有重要的理论意义与实际应用价值.
| Cai J H, Tang J T, Hua X R, et al. 2009. An analysis method for magnetotelluric data based on the Hilbert-Huang transform. Exploration Geophysics, 40(2): 197-205. DOI:10.1071/EG08124 | |
| Cai J H. 2013. Magnetotelluric response function estimation based on Hilbert-Huang transform. Pure and Applied Geophysics, 170(11): 1899-1911. DOI:10.1007/s00024-012-0620-3 | |
| Cai J H. 2014. A combinatorial filtering method for magnetotelluric time-series based on Hilbert-Huang transform. Exploration Geophysics, 45(2): 63-73. DOI:10.1071/EG13012 | |
| Chen J, Heincke B, Jegen M, et al. 2012. Using empirical mode decomposition to process marine magnetotelluric data. Geophysical Journal International, 190(1): 293-309. DOI:10.1111/gji.2012.190.issue-1 | |
| Dong S W, Li T D, Cheng X H, et al. 2012. Progress of deep exploration in mainland China:A review. Chinese J. Geophys. (in Chinese), 55(12): 3884-3901. DOI:10.6038/j.issn.0001-5733.2012.12.002 | |
| Gamble T D, Gouban W M, Clarke J. 1979. Magnetotellurics with a remote magnetic reference. Geophysics, 44(1): 53-68. DOI:10.1190/1.1440923 | |
| Goubau W M, Gamble T D, Clarke J. 1978. Magnetotelluric data analysis:removal of bias. Geophysics, 43(6): 1157-1166. DOI:10.1190/1.1440885 | |
| Hao R J, Lu W X, Chu F L. 2008. Multiscale morphological analysis on fault signals of rolling element bearing. Chinese Journal of Mechanical Engineering (in Chinese), 44(11): 160-165. DOI:10.3901/JME.2008.11.160 | |
| Huo Y L, Zhang G S, Lü S H, et al. 2013. Fractal characteristics research of lightning and its application to automatic recognition. Acta Phys. Sin. (in Chinese), 62(5): 059201. | |
| Jing J E, Wei W B, Chen H Y, et al. 2012. Magnetotelluric sounding data processing based on generalized S transformation. Chinese J. Geophys. (in Chinese), 55(12): 4015-4022. DOI:10.6038/j.issn.0001-5733.2012.12.013 | |
| Kao D W, Rankin D. 1977. Enhancement of signal-to-noise ratio in magnetotelluric data. Geophysics, 42(1): 103-110. DOI:10.1190/1.1440703 | |
| Kapple K N. 2012. A data variance technique for automated despiking of magnetotelluric data with a remote reference. Geophysical Prospecting, 60(1): 179-191. DOI:10.1111/gpr.2012.60.issue-1 | |
| Li B, Zhang P L, Ren G Q, et al. 2010. Mathematic morphology-based fractal dimension calculation and its application in fault diagnosis of roller bearings. Journal of Vibration and Shock (in Chinese), 29(5): 191-194. | |
| Li J, Tang J T, Xiao X. 2011. De-noising algorithm for magnetotelluric signal based on mathematical morphology filtering. Noise & Vibration Worldwide, 42(11): 65-72. | |
| Li J, Tang J T, Wang L, et al. 2014a. Noise suppression for magnetotelluric sounding data based on signal subspace enhancement and endpoint detection. Acta Phys. Sin. (in Chinese), 63(1): 019101. | |
| Li J, Tang J T, Xiao X, et al. 2014b. Magnetotelluric data processing based on combined generalized morphological filter. Journal of Central South University (Science and Technology) (in Chinese), 45(1): 173-185. | |
| Lü Q T, Shi D N, Tang J T, et al. 2011. Probing on deep structure of middle and lower reaches of the Yangtze Metallogenic Belt and typical ore concentration area:a review of annual progress of SinoProbe-03. Acta Geoscientica Sinica (in Chinese), 32(3): 257-268. | |
| Lü Q T, Dong S W, Tang J T, et al. 2015. Multi-scale and integrated geophysical data revealing mineral systems and exploring for mineral deposits at depth:A synthesis from SinoProbe-03. Chinese J. Geophys. (in Chinese), 55(12): 4319-4343. DOI:10.6038/cjg20151201 | |
| Maragos P. 1989. Pattern spectrum and multiscale shape representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7): 701-716. DOI:10.1109/34.192465 | |
| Maragos P, Sun F K. 1993. Measuring the fractal dimension of signals:Morphological covers and iterative optimization. IEEE Transactions on Signal Processing, 41(1): 108-121. DOI:10.1109/TSP.1993.193131 | |
| Song S G, Tang J T, He J S. 1995. Wavelets analysis and the recognition, separation and removal of the static shift in electromagnetic soundings. Acta Geophysica Sinica (in Chinese), 38(1): 120-128. | |
| Sutamo D, Vozoff K. 1991. Phase-smoothed robust M-estimation of magnetotelluric impedance functions. Geophysics, 56(12): 1999-2007. DOI:10.1190/1.1443012 | |
| Tang J T, Hua X R, Cao Z M, et al. 2008. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data. Chinese J. Geophys. (in Chinese), 51(2): 603-610. DOI:10.3321/j.issn:0001-5733.2008.02.034 | |
| Tang J T, Li J, Xiao X, et al. 2012a. Mathematical morphology filtering and noise suppression of magnetotelluric sounding data. Chinese J. Geophys. (in Chinese), 55(5): 1784-1793. DOI:10.6038/j.issn.0001-5733.2012.05.036 | |
| Tang J T, Li J, Xiao X, et al. 2012b. Magnetotelluric sounding data strong interference separation method based on mathematical morphology filtering. Journal of Central South University (Science and Technology) (in Chinese), 43(6): 2215-2221. | |
| Tang J T, Xu Z M, Xiao X, et al. 2012c. Effect rules of strong noise on Magnetotelluric (MT) sounding in the Luzong ore cluster area. Chinese J. Geophys. (in Chinese), 55(12): 4147-4159. DOI:10.6038/j.issn.0001-5733.2012.12.027 | |
| Tang J T, Liu Z J, Liu F Y, et al. 2015. The denoising of the audio magnetotelluric data set with strong interferences. Chinese J. Geophys. (in Chinese), 58(12): 4636-4647. DOI:10.6038/cjg20151225 | |
| Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of the Earth's crust. Dok1. Akad. Nauk. SSSR, 73(2): 295-297. | |
| Torre E, Picado-Muiño D, Denker M, et al. 2013. Statistical evaluation of synchronous spike patterns extracted by frequent item set mining. Frontiers in Computational Neuroscience, 7: 132. | |
| Trad D O, Travassos J M. 2000. Wavelet filtering of magnetotelluric data. Geophysics, 65(2): 482-491. DOI:10.1190/1.1444742 | |
| Vozoff K. 1972. The magnetotelluric method in the exploration of sedimentary basins. Geophysics, 37(1): 98-141. DOI:10.1190/1.1440255 | |
| Wang H, Wei W B, Jin S, et al. 2014. Removal of magnetotelluric noise based on synchronous time series relationship. Chinese J. Geophys. (in Chinese), 57(2): 531-545. DOI:10.6038/cjg20140218 | |
| Wang S M, Wang J Y. 2004. Application of higher-order statistics in magnetotelluric data processing. Chinese J. Geophys. (in Chinese), 47(5): 928-934. DOI:10.3321/j.issn:0001-5733.2004.05.027 | |
| Wei W, Cai J C, Hu X Y, et al. 2015. An electrical conductivity model for fractal porous media. Geophysical Research Letters, 42(12): 4833-4840. DOI:10.1002/2015GL064460 | |
| Wei W B. 2002. New advance and prospect of magnetotelluric sounding (MT) in China. Progress in Geophysics (in Chinese), 17(2): 245-254. | |
| Xia Y, Feng D G, Zhao R C. 2006. Morphology-based multifractal estimation for texture segmentation. IEEE Transactions on Image Processing, 15(3): 614-623. DOI:10.1109/TIP.2005.863029 | |
| Xiao X, Wang X Y, Tang J T, et al. 2014. Conductivity structure of the Lujiang-Zongyang ore concentrated area, Anhui Province:constraints from magnetotelluric data. Acta Geologica Sinica (in Chinese), 88(4): 478-495. | |
| Zhang K, Cheng J S, Yang Y. 2013. Rotating machinery fault diagnosis based on local mean decomposition and pattern spectrum. Journal of Vibration and Shock (in Chinese), 32(9): 135-140. | |
| 董树文, 李廷栋, 陈宣华, 等. 2012. 我国深部探测技术与实验研究进展综述. 地球物理学报, 55(12): 3884–3901. DOI:10.6038/j.issn.0001-5733.2012.12.002 | |
| 郝如江, 卢文秀, 禇福磊. 2008. 滚动轴承故障信号的多尺度形态学分析. 机械工程学报, 44(11): 160–165. | |
| 火元莲, 张广庶, 吕世华, 等. 2013. 闪电的分形特征研究及其在自动识别中的应用. 物理学报, 62(5): 059201. | |
| 景建恩, 魏文博, 陈海燕, 等. 2012. 基于广义S变换的大地电磁测深数据处理. 地球物理学报, 55(12): 4015–4022. DOI:10.6038/j.issn.0001-5733.2012.12.013 | |
| 李兵, 张培林, 任国全, 等. 2010. 基于数学形态学的分形维数计算及在轴承故障诊断中的应用. 振动与冲击, 29(5): 191–194. | |
| 李晋, 汤井田, 王玲, 等. 2014a. 基于信号子空间增强和端点检测的大地电磁噪声压制. 物理学报, 63(1): 019101. | |
| 李晋, 汤井田, 肖晓, 等. 2014b. 基于组合广义形态滤波的大地电磁资料处理. 中南大学学报(自然科学版), 45(1): 173–185. | |
| 吕庆田, 史大年, 汤井田, 等. 2011. 长江中下游成矿带及典型矿集区深部结构探测-SinoProbe-03年度进展综述. 地球学报, 32(3): 257–268. | |
| 吕庆田, 董树文, 汤井田, 等. 2015. 多尺度综合地球物理探测:揭示成矿系统、助力深部找矿--长江中下游深部探测(SinoProbe-03)进展. 地球物理学报, 58(12): 4319–4343. DOI:10.6038/cjg20151201 | |
| 宋守根, 汤井田, 何继善. 1995. 小波分析与电磁测深中静态效应的识别、分离及压制. 地球物理学报, 38(1): 120–128. | |
| 汤井田, 化希瑞, 曹哲民, 等. 2008. Hilbert-Huang变换与大地电磁噪声压制. 地球物理学报, 51(2): 603–610. DOI:10.3321/j.issn:0001-5733.2008.02.034 | |
| 汤井田, 李晋, 肖晓, 等. 2012a. 数学形态滤波与大地电磁噪声压制. 地球物理学报, 55(5): 1784–1793. DOI:10.6038/j.issn.0001-5733.2012.05.036 | |
| 汤井田, 李晋, 肖晓, 等. 2012b. 基于数学形态滤波的大地电磁强干扰分离方法. 中南大学学报(自然科学版), 43(6): 2215–2221. | |
| 汤井田, 徐志敏, 肖晓, 等. 2012c. 庐枞矿集区大地电磁测深强噪声的影响规律. 地球物理学报, 55(12): 4147–4159. DOI:10.6038/j.issn.0001-5733.2012.12.027 | |
| 汤井田, 刘子杰, 刘峰屹, 等. 2015. 音频大地电磁法强干扰压制试验研究. 地球物理学报, 58(12): 4636–4647. DOI:10.6038/cjg20151225 | |
| 王辉, 魏文博, 金胜, 等. 2014. 基于同步大地电磁时间序列依赖关系的噪声处理. 地球物理学报, 57(2): 531–545. DOI:10.6038/cjg20140218 | |
| 王书明, 王家映. 2004. 高阶统计量在大地电磁测深数据处理中的应用研究. 地球物理学报, 47(5): 928–934. DOI:10.3321/j.issn:0001-5733.2004.05.027 | |
| 魏文博. 2002. 我国大地电磁测深新进展及瞻望. 地球物理学进展, 17(2): 245–254. | |
| 肖晓, 王显莹, 汤井田, 等. 2014. 安徽庐枞矿集区大地电磁探测与电性结构分析. 地质学报, 88(4): 478–495. | |
| 张亢, 程军圣, 杨宇. 2013. 基于局部均值分解与形态谱的旋转机械故障诊断方法. 振动与冲击, 32(9): 135–140. | |