基于数学形态滤波的大地电磁强干扰分离方法

汤井田，李晋，肖晓，徐志敏，李灏，张弛

（中南大学 地球科学与信息物理学院，湖南 长沙，410083）

基于数学形态滤波的大地电磁强干扰分离方法

汤井田，李晋，肖晓，徐志敏，李灏，张弛

（中南大学 地球科学与信息物理学院，湖南 长沙，410083）

针对矿集区大地电磁信号采集过程中常引入强噪声干扰等问题，采用数学形态滤波对大地电磁强干扰分离方法进行研究。在仿真信号中加入常见的强干扰来检验形态滤波的降噪能力，根据噪声类型选取不同结构元素尺寸及大小，并将形态滤波应用于实测大地电磁数据的降噪处理。采用非线性共轭梯度法进行反演，考查形态滤波对提高大地电磁测量数据质量的改善情况。研究结果表明：数学形态滤波能有效消除大地电磁强干扰中的大尺度干扰和基线漂移现象，重构信号基本保留原始大地电磁信号特征，改善大地电磁测深数据质量。由于该方法原理简单、并行运算速度快，具有较好的应用价值，适合于矿集区海量大地电磁强干扰分离。

形态滤波；大地电磁；强干扰分离；结构元素

矿产资源是国民经济和社会发展的重要物质基础，随着地表资源勘查程度的提高，矿床发现的难度越来越大，必须增加勘探深度，向深部索取资源，以不断满足我国经济高速发展对资源的需求。大地电磁测深法（Magnetotelluric，MT）自20世纪50年代初提出以来，随着科学技术的不断进步，尤其是计算机技术及现代数字信号处理技术的发展，使得MT得到了飞速发展，已成为矿区找矿、海洋地质调查、工程勘

探、地震预报、油气田普查勘探、岩石圈结构探测、地下水和寻找地热资源等领域的重要手段。大地电磁测深理论提出至今，噪声问题一直困扰着广大MT研究者，如何抑制噪声，提高MT测量数据质量，是国内外共同关注的课题[1−3]。目前大地电磁去噪方法存在一定的局限性，如远参考处理法从理论与实践结果看，能成功消除局部相关噪声，但是，经远参考处理后，单点数据的误差棒变大，特别是在处理受电磁干扰严重和校正量较大的数据段时，该现象尤为明显[4]；Robust处理法无法消除输入端的噪声，且Robust法无法剔除噪声较多和能量较强时的近源电磁相关噪声对数据的干扰[5]；小波变换则过分依赖于小波函数的选取，有时随着尺度增大，相应正交基函数的频谱局部性变差，使其对大地电磁信号更精细分解受到限制[6]；HHT法虽然不要选择基函数，与小波变换相比具有更强的时频刻画能力，但因EMD分解是自适应的，无法揭示每时段的频率特性和能量差异所具有的细微变化，且算法占用大量运算时间，不适合对实测大地电磁信号进行处理[7]；人机联作法虽能较好地改善数据质量，但操作时参与了太多的人为因素，而且耗费很多时间和精力，且操作者必须具备丰富的噪声识别经验，否则效果适得其反[8]。由于矿集区人口稠密、现代通讯设备先进、交通发达及矿产资源丰富、民间开采广泛等因素导致MT数据采集与处理相当困难。大地电磁场具有很宽的频率范围，鉴于其激发机制不同，大地电磁信号表现出的形态和频谱特征也各有差异，且极易受各种噪声干扰，是典型的非线性和非平稳信号。由于野外观测数据不可避免地受到各种噪声的污染，且电磁噪声干扰日益严重，导致大地电磁采集数据的一致性和重复性变差，曲线参数的鲁棒性也变差，使估算分散不能客观地反映地下电性分布情况，甚至得到错误的结论。如何应用现代数字信号处理技术剖析大地电磁强干扰信号特征，得到无偏的阻抗估计是取得良好勘探结果的关键。尤其是在矿集区，矿山的开采、冶炼及与其配套的重工业密集形成的电磁噪声和人文噪声非常复杂，导致测量得到的视电阻率失真严重。因此，研究MT噪声的特点和MT测量数据中的信噪识别和处理技术，对改善MT测量数据质量、探测地壳精细结构及寻找深部控矿构造具有非常重要的实际应用价值。在此，本文作者针对大地电磁噪声的上述特点，引入数学形态滤波进行探讨，针对不同的噪声类型选择不同的结构元素尺寸及大小，构造适合大地电磁强干扰分离的形态滤波器，并对矿集区实测大地电磁强干扰进行降噪处理。

1 数学形态滤波

数学形态学（Mathematical morphology，MM）是基于积分几何和随机集合论等数学理论建立起来的1种非线性信号处理方法。该算法是用集合来描述目标，集合各部分之间的关联表明目标信号的结构特点，即考察目标信号时设计1个称为结构元素的“探针”，通过探针在信号中不停地移动来考察信号各部分之间的关联，从而提取出有用的信息进行结构分析和描述。该方法利用结构元素对信号的几何特征进行局部匹配或修正，同时保留目标信号主要的形状特征，以达到提取有用信息、保持细节成分和抑制噪声的目的。形态滤波器是从数学形态学发展起来的一种新型非线性滤波技术，近年来，随着形态学理论的飞速发展，形态学滤波被逐步推广到一维信号处理领域[9−11]。

1.1 基本形态变换

数学形态学的基本变换为膨胀和腐蚀，膨胀和腐蚀变换等价于离散函数在滑动滤波窗（相当于结构元素）内的最大值和最小值滤波。其中，膨胀变换表示1种扩张变换，使目标肢体扩张，孔洞收缩，用来填平边界不平滑的凹陷部分，增大了谷值，扩展了峰顶。腐蚀变换表示1种收缩变换，使目标肢体收缩，孔洞扩张，用来剔除边界不平滑的凸起部分，减少了峰值，加宽了谷域。其定义如下。

设输入信号f（n）为定义在Df={x0，x1，…，xn}上的离散函数，结构元素g（m）为定义在Dg={y0，y1，…，ym}上的离散函数，且n＞＞m。则f（n）关于g（m）的膨胀和腐蚀运算分别定义为：

式中：⊕和⊙分别表示膨胀和腐蚀运算。

形态开和形态闭运算是膨胀和腐蚀运算的衍生运算，f（n）关于g（m）的形态开和形态闭运算分别定义为：

式中：“◦”和“·”分别表示形态开和形态闭运算。其中：形态开运算是对结构元素先腐蚀后膨胀，用来消除目标信号中的细节、孤立点及叠加在信号上窄的“毛刺”，目的是使目标信号的轮廓光滑，从而剔除尖峰，抑制正脉冲（峰值）噪声；形态闭运算是对结构元素先膨胀后腐蚀，用来弥补目标信号中的孔洞及很窄的“裂缝”，目的是填平谷底，抑制负脉冲（低谷）噪声。开、闭运算均具有低通特性，在实际应用中，常采用形态开、闭运算或形态开、闭运算的级联形式相结合来滤除目标信号中的正、负脉冲噪声。

十六大以来，2002年到2012是我国GDP增长最快的十年，国内GDP增速超过两位数，也一跃成为世界第二大经济体。然而从2012年开始，中国经济已然告别了高速增长期，开始了被称之为＂新常态＂的新阶段，GDP增长率也从2011年接近百分之十下降到2016年的6.7%，这已经是自上世纪九十年代以来所出现的最低增速。

1.2 构建形态滤波器

Maragos和Schafer采用同一类型和尺寸的结构元素，通过不同顺序级联形态开、闭运算，定义了形态开−闭FOC（f（n））和闭−开滤波器FOC（f（n））[12−13]：

由式（5）和（6）可知：形态开−闭和形态闭−开滤波器都能同时去除目标信号中的正、负脉冲干扰，具有形态学中开、闭运算的所有性质。但是，由于形态开运算的收缩性导致形态开−闭滤波器输出偏小，而形态闭运算的扩展性导致形态闭−开滤波器输出偏大。因此，目标信号在滤波过程中存在统计偏倚现象[14]，直接影响到滤波器的噪声抑制性能。为了有效去除各种噪声干扰和抑制统计偏倚现象，采用级联开、闭运算构造形态开−闭和形态闭−开组合形态滤波器作为输出，用于大地电磁强干扰噪声的分离：

式中：y（n）表示组合形态滤波器的输出结果。重构的大地电磁信号ε（n）定义为：

1.3 选取结构元素

形态滤波的质量取决于所选择的形态变换和结构元素。结构元素在形态运算中的作用类似于一般信号处理中的滤波窗口或参考模板，其尺寸和形状都将对形态学基本变换产生很大影响。采用不同的结构元素可以提取出目标信号中不同的形状特征，选取的结构元素要尽可能接近待分析信号的形状特点。常见的结构元素有直线型、三角型、圆盘型、正弦型、抛物线型以及其他多边形组合。待处理信号的形状通常决定了结构元素的形状设计，一般一种结构元素对一类噪声有较好的滤除效果。相对而言，结构元素越复杂，其滤除噪声的能力就越强，但所耗费的时间也就越长[15]。

2 矿集区典型强干扰类型分析

图1所示为安徽省庐枞矿集区某测点的一段实测大地电磁原始数据，测量数据由加拿大凤凰公司生产的V5-2000仪器采集。对该段数据出现强干扰的时间序列进行分析可知：电道和磁道均不同程度地受到了周期性的突跳和波动等信号干扰，这些信号与稳定的天然电磁场信号相比，具有振幅大、能量强和周期性明显等特征。因此，将时域信号中出现的明显非天然电磁场的信号定为噪声信号。对大量MT测量数据的时间序列特征进行分析可知：电道数据常被强度大的方波噪声干扰，造成电道曲线整体基线漂移严重，对低频数据的影响巨大；磁道数据常被大尺度充放电三角波噪声干扰，造成磁道曲线突跳明显且噪声局部能量强，将正常磁场信号几乎完全湮灭，且电道与磁道干扰出现的时刻具有一定的相关性[16]。这2类典型的强干扰噪声大量出现在矿集区大地电磁测深数据中，得到的视电阻率曲线往往表现为明显的近源效应。

图1 矿集区典型强干扰类型Fig.1 Typical strong interference types in ore district

3 数值仿真试验

针对矿集区大地电磁信号的噪声类型，模拟典型的含大尺度方波和充放电三角波的强干扰信号进行仿真试验。根据形态滤波器的降噪性能，选择不同的结构元素及尺寸进行分析。降噪性能评价标准采用滤波误差E和信噪比RSN2组参数进行衡量，分别定义如下：

误差E定义为：

式中：s（n）为原始信号；y（n）为形态滤波输出信号。

图2所示为含大尺度方波的MT信号的形态滤波降噪效果。图3所示为含大尺度充放电三角波的MT信号的形态滤波降噪效果。

图2 含大尺度方波的MT信号形态滤波效果图Fig.2 Morphology filtering effect chart for large-scale square wave MT noise

图3 含大尺度充放电三角波的MT信号形态滤波效果图Fig.3 Morphology filtering effect chart for large-scale charging and discharging triangle wave MT noise

由图2和图3可知：含大尺度方波和充放电三角波的大地电磁强干扰的能量幅值往往达到正常大地电磁信号的几十倍，正常有用信号几乎被完全湮没。经形态滤波器后，提取出的形态轮廓清晰、平滑，几乎完整地勾勒出含方波和充放电三角波的大尺度噪声轮廓曲线，将小于或等于结构元素的信号滤除，只保留了比结构元素大的信号单元，重构后的大地电磁信号有效地剔除了叠加在有用信号上的大尺度干扰和由噪声引起的突跳波形，抑制了基线漂移现象，凸出了大地电磁有用信号的相关局部特征，大尺度干扰与正常信号得到了有效分离。

表1和表2所示为含大尺度方波的MT信号采用不同结构元素及同一结构元素不同尺寸的滤波性能对比结果。

表1 含大尺度方波的MT信号采用不同结构元素滤波性能对比Table 1 Filtering performance comparison by different SEs for large-scale square wave MT noise

表2 含大尺度方波的MT信号采用不同尺寸的矩形结构元素滤波性能对比Table 2 Filtering performance comparison by different sizes of rectangular type SE for large-scale square wave MT noise

表3和表4所示为含大尺度充放电三角波的MT信号采用不同结构元素及同一结构元素不同尺寸的滤波性能对比结果。

结合滤波误差E和信噪比RSN，分析表1～ 4的试验结果可见：待处理信号的形状与结构元素的形状选取之间有一定的相似性，当电道中出现大尺度方波干扰时，选择长度为15点和幅值为0.2的矩形结构元素的滤波效果较好，当磁道中出现大尺度充放电三角波干扰时，选择长度为9点和幅值为0.2的三角形结构元素的滤波效果较好。而且，试验发现：在选取结构元素的长度时，并不是结构元素的长度越长滤波效果就越好，当结构元素长度越长时，计算速度反而越慢。以上数值仿真试验表明：数学形态滤波对大地电磁强干扰具有较好的去噪能力。

表3 含大尺度充放电三角波的MT信号采用不同结构元素滤波性能对比Table 3 Filtering performance comparison by different SEs for large-scale charging and discharging triangle wave MTnoise

表4 含大尺度充放电三角波的MT信号采用不同尺寸的三角形结构元素滤波性能对比Table 4 Filtering performance comparison by different sizes of triangle type SE for large-scale charging and discharging triangle wave MT noise

4 实测资料应用

4.1 大地电磁信号数据提取

图4（a）所示为该矿集区采集的实测大地电磁信号在Synchro Time Series View 中观测到的Ex低频段数据。V5-2000采集数据时在每个测点记录Ex，Ey，Hx，Hy和Hz5个分量，并将原始数据记录在TSH（高频）和TSL（低频）文件中。由于仪器不提供读取时间序列的软件，资料的处理相当于一个黑匣子。因此，首先必须获取原始时间序列，才能进一步分析和处理大地电磁强噪声干扰。

图4（b）所示为根据V5-2000的数据采集格式[17]，在VC环境下进行编程，从TSH和TSL文件中读取出的原始数据时间序列在Matlab中的观测结果。对比图4（a）和（b）可知：读取的原始数据时间序列与仪器采集的数据一致，这样就保证了后续资料处理的可靠性。

4.2 实测数据分析

鉴于大地电磁信号数据量庞大，仅选取其中1个测点4个分量（Ex，Ey，Hx和Hy）中的1个电场分量（Ex）和1个磁场分量（Hy）进行分析。

图5所示为庐枞矿集区中某测点Ex信号经形态滤波后的去噪效果图。

图6所示为该测点Hy信号的去噪效果图。

图4 实测大地电磁测深数据Fig.4 Measured MT signal

图5 实测大地电磁信号Ex的去噪效果图Fig.5 Morphology filtering effect chart for measured MT signal in Ex

图6 实测大地电磁信号Hy的去噪效果图Fig.6 Morphology filtering effect chart for measured MT signal in Hy

分析图5和图6可知：含大尺度方波和充放电三角波的大地电磁强噪声干扰经形态滤波后，提取出强干扰的大尺度形态轮廓，且曲线自然、光滑，较好地保持了原始信号的局部特征。重构后的大地电磁信号则基本去除了由噪声引起的突跳波形，有效地抑制了大地电磁信号中的大尺度干扰和基线漂移，保留了较为丰富的细节成分，基本重现了大地电磁正常信号的原始特征。

图7所示为庐枞矿集区某测线采用非线性共轭梯度法TE模式形态滤波前后的反演效果对比图，图8所示为采用非线性共轭梯度法TM模式形态滤波前后的反演效果对比图。

分析图7和图8可知：形态滤波前后地下介质基本形态一致，形态滤波前由于个别测点受到近源干扰导致局部极值产生假高阻异常，从而影响周围测点的反演结果，形态滤波后尖锐轮廓变得光滑，大大消除了近源干扰，从而使假高阻异常体减少，反演结果更加真实合理。

图7 非线性共轭梯度法TE模式反演效果图Fig.7 Nonlinear conjugate gradient method TE model inversion effect chart

图8 非线性共轭梯度法TM模式反演效果图Fig.8 Nonlinear conjugate gradient method TM model inversion effect chart

5 结论

（1） 数学形态滤波能几乎完整地勾勒出大地电磁强干扰中大尺度噪声的形态轮廓，通过分析待处理信号的形态特征，合理地选取结构元素的类型和尺寸，能有效抑制大尺度干扰和基线漂移现象。

（2） 形态滤波法能较大改善矿集区大地电磁测深数据质量，对电磁法探测结果的处理和反演解释具有重要意义。

（3） 数学形态滤波法原理简单，具有并行运算速度快的优势，适合于矿集区海量大地电磁强干扰分离，为提高大地电磁测深的深部探测能力提供了新的强有力的技术支持，应用前景广阔。

[1] Cai J H, Tang J T. An analysis method for magnetotelluric data based on the Hilbert-Huang transform[J]. Exploration Geophysics, 2009, 40（2）：197−205.

[2] Trad D O, Travassos J M. Wavelet filtering of magnetotelluric data[J]. Geophysics, 2000, 65（2）： 482−491.

[3] 汤井田, 化希瑞, 曹哲民, 等. Hilbert-Huang变换与大地电磁噪声压制[J]. 地球物理学报, 2008, 51（2）： 603−610.

TANG Jing-tian, HUA Xi-rui, CAO Zhe-ming, et al. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J]. Chinese J Geophys, 2008, 51（2）： 603−610.

[4] Gamble T M, Gouban W M, Clarke J. Magnetotelluric with a remote magnetic reference[J]. Geophysics, 1979, 44： 53−68.

[5] Larsen J C, Maekie R L. Robust smooth magnetotelluric transfer functions[J]. Geophys J Int, 1996, 124： 801−819.

[6] 宋守根, 汤井田, 何继善. 小波分析与电磁测深中静态效应的识别、分离及压制[J]. 地球物理学报, 1995, 38（1）：120−128.

SONG Shou-gen, TANG Jing-tian, HE Ji-shan. Wavelets analysis and the recognition, separation and removal of the static shift in electromagnetic soundings[J]. Chinese J Geophys, 1995, 38（1）： 120−128.

[7] 蔡剑华, 龚玉蓉, 王先春. 基于Hilbert-Huang变换的大地电磁测深数据处理[J]. 石油地球物理勘探, 2009, 44（5）： 617−621, 629.

CAI Jian-hua, GONG Yu-rong, WANG Xian-chun. Magnetotelluric data processing based on Hilbert-Huang transform in oil and gas exploration[J]. OGP, 2009, 44（5）： 617−621, 629.

[8] 范翠松. 矿集区强干扰大地电磁噪声特点及去噪方法研究[D].长春： 吉林大学地球探测与科学技术学院, 2009： 43−55.

FAN Cui-song. The strong noise characteristics of MT in ore concentration area and research of denoise method[D]. Changchun： Jilin University. College of Geo-Exploration Science and Technology, 2009： 43−55.

[9] 张文斌, 杨辰龙, 周晓军. 形态滤波方法在振动信号降噪中的应用[J]. 浙江大学学报： 工学版, 2009, 43（11）： 2096−2099.

ZHANG Wen-bin, YANG Chen-long, ZHOU Xiao-jun. Application of morphology filtering method in vibration signal de-noising[J]. Journal of Zhejiang University： Engineering Science, 2009, 43（11）： 2096−2099.

[10] 陈辉, 郭科, 胡英. 数学形态学在地震信号处理中的应用研究[J]. 地球物理学进展, 2009, 24（6）： 1995−2002.

CHEN Hui, GUO Ke, HU Ying. A study on application of mathematical morphology to seismic signal processing[J]. Progress in Geophys, 2009, 24（6）： 1995−2002.

[11] 岳蔚, 刘沛. 基于数学形态学消噪的电能质量扰动检测方法[J]. 电力系统自动化，2002, 26（7）： 13−17.

YUE Wei, LIU Pei. Detection of power quality disturbances based on mathematical morphology filter[J]. Automation of Electric Power Systems, 2002, 26（7）： 13−17.

[12] Maragos P, Schafer R W. Morphological filters-Part I： Their set theoretic analysis and relation to linear shift invariant filters[J]. IEEE Trans on ASSP, 1987, 35（8）： 1153−1169.

[13] Maragos P, Schafer R W. Morphological filters-Part II： Their relation to median, order-statistic and stack filters[J]. IEEE Trans On ASSP, 1987, 35（8）： 1170-−1184.

[14] 赵静, 何正友, 钱清泉. 利用广义形态滤波与差分熵的电能质量扰动检测[J]. 中国电机工程学报, 2009, 29（7）： 121−126.

ZHAO Jing, HE Zhen-you, QIAN Qing-quan. Detection of power quality disturbances utilizing generalized morphological filter and difference entropy[J]. Proceedings of the CSEE, 2009, 29（7）： 121−126.

[15] 张建成, 吴新杰. 形态滤波在实时信号处理中应用的研究[J].传感技术学报, 2007, 20（4）： 828−831.

ZHANG Jian-cheng, WU Xin-jie. Research on applications of morphological filtering in real-time signal processing[J]. Chinese Journal of Sensors and Actuators, 2007, 20（4）： 828−831.

[16] 王大勇. 长江中下游矿集区综合地质地球物理研究-以九瑞、铜陵矿集区为例[D]. 长春： 吉林大学地球探测与科学技术学院, 2010： 64−70.

WANG Da-yong. The integrated geophysical and geological study in the ore belt of the middle and lower reach of the Yangtze river-The cases study of Tongling and Jiurui ore district[D]. Changchun： Jilin University. College of Geo-Exploration Science and Technology, 2010： 64−70.

[17] 柳建新, 刘春明, 马捷. V5-2000大地电磁测深仪文件头数据格式研究[J]. 物探与化探计算技术, 2007, 29（4）： 359−362.

LIU Jian-xin, LIU Chun-ming, MA Jie. The data format of MT sounding instrument V5-2000[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2007, 29（4）： 359−362.

（编辑 何运斌）

Magnetotelluric sounding data strong interference
separation method based on mathematical morphology filtering

TANG Jing-tian, LI Jin, XIAO Xiao, XU Zhi-min, LI Hao, ZHANG Chi
（School of Geosciences and Info-Physics, Central South University, Changsha 410083, China）

Aimed at the problem that magnetotelluric signals are frequently influenced by strong noise interference during acquisition, a new magnteotelluric sounding data strong interference separation method based on mathematical morphology filtering was proposed. Firstly, the noise reduction capability was checked through simulated signal with the common strong interference and then the structure elements size selection program was analysed according to different noise types. Finally, the morphological filtering method was applied to the measured magnetotelluric signals noise reduction process. Experiments using nonlinear conjugate gradient method for inversion were conducted to check the improvement of magnetotelluric signals quality. The results indicate that the proposed method can effectively eliminate large scale disturbance and baseline drift of magnetotelluric signals. The method can better keep the original features of magnetotelluric signals, and the quality of magnteotelluric sounding data is improved. Because the principle of the method is simple and the parallel computing speed is fast, it has good application value and is suitable for the strong interference separation in ore district.

morphology filtering; magnetotelluric sounding; strong interference separation; structuring element

P631

A

1672−7207（2012）06−2215−07

2011−06−05；

2011−08−02

国家科技专项“深部探测技术与实验研究专项”资助项目（SinoProbe-03）；国家自然科学基金资助项目（41104071）；湖南省教育厅资助项目（11B074）；中南大学中央高校基本科研业务费专项资金资助项目（2011QNZT012）

李晋（1981−），男，湖南衡东人；博士研究生，从事大地电磁强噪声压制及信号处理研究；电话：13574170783；E-mail：geologylj@163.com