子空间增强在大地电磁信噪分离中的应用研究

燕 欢， 李 晋,2， 李春阳

(1.湖南师范大学 物理与信息科学学院，长沙 410081; 2.中南大学 有色金属成矿预测与地质环境监测教育部重点实验室， 地球科学与信息物理学院，长沙 410083)

子空间增强在大地电磁信噪分离中的应用研究

燕 欢1， 李 晋1,2， 李春阳1

(1.湖南师范大学 物理与信息科学学院，长沙 410081; 2.中南大学 有色金属成矿预测与地质环境监测教育部重点实验室， 地球科学与信息物理学院，长沙 410083)

为了提高强干扰地区大地电磁测深数据质量、保留大地电磁信号低频段的有用信息，压制大尺度强干扰已刻不容缓。针对矿集区典型的强干扰类型，引入子空间增强方法对大地电磁信号与强干扰进行信噪分离研究。通过计算机模拟常见的大尺度方波和充放电三角波干扰，讨论了子空间增强算法中帧长与特征值判别阈值的最优选取范围，给出了算法流程，并对矿集区实测大地电磁数据进行了信噪分离处理。实验结果表明：该方法是切实可行的，能有效剔除时间域序列中的大尺度强干扰，近源效应得到了有效压制；经处理后，重构的大地电磁信号中包含了更为丰富的细节成分，视电阻率曲线更为光滑、连续，低频段的大地电磁数据质量得到了明显改善。

大地电磁； 信噪分离； 子空间增强； 强干扰

0 引言

大地电磁测深(Magnetotelluric，简称MT)法是由Tikhonov和Cagniard[1-2]提出的以天然交变电磁场为场源的电磁勘探方法。由于天然电磁场信号弱、频带宽，因此极易受到各种噪声的干扰，是典型的非线性、非平稳信号[3-5]。随着人类文明的不断进步和经济的高速发展，人文电磁干扰日益严重，尤其在矿集区，电话线路、无线通讯基站、采矿用的大功率直流电机车等引起强能量的电磁干扰，导致原本微弱的大地电磁信号几乎完全被湮没在这些“大尺度”干扰中。诸多现代信号处理方法(如Hilbert-Huang变换[6]、广义S变换[7]、方差比维纳滤波[8]、数学形态滤波[9]、同步时间序列[10]等),被应用到大地电磁噪声压制，均在一定程度上改善了大地电磁测深数据质量，推动了大地电磁测深法在各种噪声中的发展。

子空间增强是一种非线性、非平稳信号分析方法，上世纪90年代该方法已被应用于语音增强领域。由于语音信号的噪声方差可以估计，根据子空间增强算法原理，带噪语音信号的向量空间可分解为信号子空间和噪声子空间，通过剔除噪声子空间，再利用信号子空间即可恢复较为纯净的语音信号。子空间增强法现已被逐步推广到故障诊断、地震资料处理等领域[11-12]。由于大地电磁信号是典型的非线性、非平稳信号，为此借助语音增强、地震资料处理等领域广泛应用的子空间增强法，提出分离矿集区微弱大地电磁信号和典型强干扰的研究思路，分析了关键参数的最优选取方案，给出了具体的算法处理流程。

1 子空间增强基本原理

子空间分解的降噪方法早期由Dendrinos[13]提出，该算法利用带噪信号组成具有Toeplitz结构的矩阵，随后对该矩阵进行奇异值分解；根据最小二乘估计，忽略较小的奇异值后可将降维得到的数据根据Toeplitz矩阵结构恢复出原信号。然而，子空间增强法具有一定的局限性，只能处理白噪声、不能处理有色噪声。在此基础上，Ephraim等[14]做了进一步完善，利用信号协方差矩阵的特征值分解(Eigen value decomposition，EVD)，将带噪信号的向量空间分解为噪声子空间和信号子空间两个相互正交的子空间再进行语音增强处理。随后，Jensen等[15]将基于SVD(Singular value decomposition，SVD)的子空间算法扩展到有色噪声领域。

子空间增强具有控制信号失真和噪声残留的特征，其基本原理是将受干扰的带噪信号映射到两个子空间①信号子空间；②噪声子空间，通过数学方法将噪声子空间全部置零，滤除信号子空间中所包含的噪声干扰，再利用信号子空间恢复出较为纯净的原始信号[16]。

以一维离散信号为例，假定噪声是加性的，且与纯净信号相互独立，则带噪信号可表示为：

y=x+n

(1)

式中：y、x和n分别表示带噪信号、纯净信号和加性噪声的向量表示,对其进行分帧，可得到相应的K维向量，其带噪信号、纯净信号和噪声的协方差矩阵Ry、Rx和Rn关系如下：

Ry=E[yyT]=E[xxT]+E[nnT]=

Rx+Rn

(2)

式中:E[·]表示均值运算。

假设纯净信号估计为：

(3)

式中:H表示K×K维的线性估计器矩阵，则该估计器的真实值与估计值的误差ε为：

εx+εn

(4)

式中：εx和εn分别表示信号的失真向量和残余噪声向量。子空间增强即通过线性估计投影到信号子空间，从而得到原始信号的较好估计。

tr(HRxHT-HRx-RxHT+Rx)

(5)

tr(HRnHT)

(6)

通过求解以下时域约束条件方程，可获得优化的线性估计器：

(7)

式中:σ2为正常数。在可接受的噪声残差水平σ2下，该估计器能最小化信号失真。

上述方程的解为：

Hopt=Rx(Rx+μRn)-1

(8)

式中:μ为拉格朗日乘子。

由Rx的特征值分解为：

Rx=UΛxUT

(9)

式中：U表示Rx的归一化特征向量矩阵；Λx表示由Rx的特征值组成的对角矩阵：

(10)

Hopt=UΛx(Λx+μUTRnU)-1UT

(11)

当噪声n为方差σn的加性白噪声时，Rn可表示为：

Rn=σnI

(12)

当噪声n为有色噪声时，往往用对角矩阵Λn来近似UTRnU：

(13)

为此，Hopt可优化为式(14)。

Hopt=UΛx(Λx+μΛn)-1UT

(14)

2 仿真信号分析

2.1 子空间增强关键参数分析

观测矿集区海量大地电磁数据的时间域波形可知，大尺度方波和充放电三角波干扰通常出现在各种采样数据中，是最为典型且严重的噪声干扰类型[17]。为了验证文中所提方法的可行性，计算机模拟大尺度方波和充放电三角波进行仿真实验。大尺度方波通过选取Matlab自带的block信号叠加4dB的白噪声来模拟(图1)。大尺度充放电三角波通过选取幅值为5%、峰谷距为10、间距为100的充放电三角波叠加随机噪声来模拟(图2)。

图1 模拟大尺度方波干扰Fig.1 Simulated large-scale square interferences

图2 模拟大尺度充放电三角波干扰Fig.2 Simulated large-scale charging and discharging triangular interferences

由于子空间增强算法中分帧帧长和协方差矩阵特征值的判别阈值对滤波精度尤为关键，为此引入波形相似度对该关键参数进行性能评价。相似度的定义为式(15)。

(15)

式中:f(n)、g(n)分别表示两个离散序列。NCC的值在-1～1之间，-1表示变换前后波形反向，“0”表示正交，“1”表示完全相同；NCC越接近“1”，则表示波形之间越相似。

从图3可知，大尺度方波的帧长取值在0～8范围内相似度逐渐上升；当帧长取值在8～18时相似度趋于平稳，且达到最大值0.9；帧长取值在18以后，相似度逐渐下降。大尺度充放电三角波的帧长取值在0～10范围时，相似度逐渐上升；在10～18时，相似度接近于0.8；帧长为18以后，相似度逐渐下降。分析图4可知，大尺度方波的判别阈值在10～100时，相似度在0.75～0.8之间；大尺度充放电三角波具有同样的变化趋势，其相似度在0.9～0.95之间。由于模拟的大尺度强干扰相比实测大地电磁信号要单一得多，而判别阈值本身即为噪声残留与信号失真的折中，且每一帧噪声的特征值变化没有实测信号复杂。综上分析，当帧长选取10～18、判别阈值选取40～60时，滤波效果最优。

图3 不同帧长滤波性能对比Fig.3 Filtering performance comparison of different frames

图4 不同判别阈值滤波性能对比Fig.4 Filtering performance comparison of different criterion thresholds

2.2 算法流程

根据子空间增强基本原理，将该算法应用到模拟大尺度信噪分离，其具体步骤如下：

1)将模拟带噪信号分段，每段为一帧，帧长在最优范围内选取。

2)构造协方差矩阵，对其进行特征值分解。

3)假设大尺度噪声与微弱信号互不相关，设置最优判别阈值，分离出噪声子空间和信号子空间。

4)将噪声子空间的特征值置零，重构原始信号。

图5所示为模拟大尺度干扰经子空间增强后的去噪效果图。文中帧长为12，判别阈值为40。

从图5可知，叠加在大尺度干扰上的信号极其微弱，经文中所提方法处理后，提取的大尺度干扰轮廓较为光滑、连续，重构信号中较好地保留了更为丰富的细节成分。

图5 模拟大尺度干扰去噪效果图Fig.5 Denoising effect of simulated large-scale interferences(a)模拟大尺度干扰(方波);(b)子空间分离出的大尺度干扰(方波);(c)重构信号(方波);(d)模拟大尺度干扰(充放电三角波);(e)子空间分离出的大尺度干扰(充放电三角波);(f)重构信号(充放电三角波)

3 实测数据处理

3.1 时间域去噪效果

为了验证子空间增强法的实用性，对矿集区受强干扰严重的大地电磁数据进行噪声压制。由于矿集区大地电磁数据量庞大，文中给出某测点在同一时刻电道(EX、EY)和磁道(HY、HX)中四道时间域片段的去噪效果，如图6所示。

从图6可知，由于该测点所处地区噪声干扰源众多，导致时间域波形中出现大量能量强、频带宽的大尺度类方波、类充放电三角波等噪声干扰；波形的整体形态呈不连续状态、跳跃性强，且电道EX/EY和磁道HY/HX出现干扰的时刻具有一定的相关性，微弱的大地电磁有用信号几乎被完全湮没。经子空间增强法处理后，大尺度干扰的轮廓特征几乎被完整地提取出来，分离出的大地电磁微弱信号在剔除基线漂移的同时，保留了有用信号更多的细节信息。由于实测大地电磁信号所含噪声种类众多，即使是同种形态的噪声其幅值、间距也是多种多样的。在如此复杂的干扰环境下，子空间增强算法仍能较清晰地提取出叠加在原始大地电磁信号上的大尺轮廓，说明该方法在时间域对大尺度干扰的分离是可行的。

3.2 视电阻率曲线

选用加拿大凤凰公司生产的V5-2000大地电磁测深仪采集的数据进行分析，该仪器采样频率通常为高频(2 560 Hz)、中频(320 Hz)和低频(24 Hz)，数据的存储格式为TSH、TSL，其中TSH文件记录了2 560 Hz和320 Hz采样率的中高频数据，TSL文件记录了24 Hz采样率的低频数据[18]。由于高频和中频是交替采集，每5 min采集1次高频或中频，其中只有1 s的高频数据和连续8 s的中频数据，所含信息量较少；低频数据为全时间段采集，能记录丰富的大地电磁数据，且更能反映测点深部的电性结构信息，为此我们重点对低频信号(TSL文件)进行去噪处理。

图7所示为庐纵矿集区B40993测点的原始时间域波形(EX、EY、HX、HY)和视电阻率曲线。

从图7(a)可知，原始数据电道EX、EY中出现大量大尺度类方波干扰，磁道HX、HY中出现大量大尺度类充放电三角波干扰。分析图7(b)可知，原始数据的视电阻率曲线整体形态连续性差。在大于40 Hz时，视电阻率曲线比较平稳，且变化趋势一致；在0.05 Hz～50 Hz，视电阻率曲线呈45°左右渐近线快速上升，0.05 Hz左右视电阻率值达到最大值接近于100 000 Ω·m，表现为典型的近源效应；在0.005 Hz～0.05 Hz，视电阻率曲线突然下降至100 Ω·m左右，误差棒增大，并出现不同程度的突跳与畸变。分析原始数据的时间域波形和视电阻率曲线可知，该测点受到了各种复杂干扰源的影响，其数据已不能客观反映真实的地电信息。虽然功率谱筛选具有较强的适用性，但对矿集区高强度近源干扰也无能为力。

图6 实测大地电磁数据去噪效果Fig.6 Denoising effect of measured magnetotelluric data(a)原始实测信号EX;(b)子空间分离出的大尺度干扰EX;(c)重构大地电磁信号EX;(d)原始实测信号EY;(e)子空间分离出的大尺度干扰EY;(f)重构大地电磁信号EY;(g)原始实测信号HY;(h)子空间分离出的大尺度干扰HY;(i)重构大地电磁信号HY;(j)原始实测信号HX;(k)子空间分离出的大尺度干扰HX;(l)重构大地电磁信号HX

图8为子空间增强法得到的时间域波形和视电阻率曲线。其中，图8(a)为经子空间增强法处理后，重新还原至V5 System 2000观测的时间域波形。

从图8(a)可知，充斥在时间域中的大尺度类方波和类充放电三角波已被基本剔除、基线漂移得到了有效压制，处理后的时间域波形基本连续、无明显跳变。分析图8(b)可知，经子空间增强法处理后，0.05 Hz～50 Hz呈近45°上升的近源趋势完全消失，曲线变得光滑、连续、同步；在0.05 Hz～5 Hz，视电阻率最大值下降了近3个数量级，近15个频点的数据质量得到了明显改善。分析对比图7可知，时、频域的处理更加真实、合理、置信，其结果已逐步接近于测点本身所固有的地下介质电性结构。

图9所示为在子空间增强法基础上做简单的功率谱筛选得到的视电阻率曲线。

分析图9可知，经简单功率谱筛选后，视电阻率曲线在0.000 5 Hz～0.05 Hz的突跳频点得到了较好恢复，曲线下掉的趋势抬升显著，且误差棒减小；视电阻率曲线的整体形态更为光滑、连续，视电阻率曲线相对稳定。实验结果表明，经子空间增强法处理后，受近源干扰严重的测点其低频段的整体数据质量得到了明显改善。

图7 测点B40993原始大地电磁数据Fig.7 Original MT data at site B40993(a)时间域波形;(b)视电阻率曲线

图8 子空间增强方法Fig.8 The subspace enhancement method(a)时间域波形;(b)视电阻率曲线

图9 经简单功率谱筛选的视电阻率曲线Fig.9 Apparent resistivity curve through simple power spectrum filtering

4 结论

1)将子空间增强算法应用于矿集区大地电磁信噪分离，通过模拟典型的强干扰类型，分析了算法中影响滤波精度的帧长及特征值判别阈值的最优选值范围，给出了算法流程。

2)通过模拟信号仿真和实测数据处理，验证了方法的可行性。受近源干扰严重的测点经子空间增强法处理后，时间域的大尺度干扰和基线漂移得到了有效压制，视电阻率曲线的整体形态更为光滑、连续，低频段的大地电磁数据质量得到了明显改善。子空间增强法对提升矿集区大地电磁测深深部勘探能力，具有重要的实际应用价值。

[1] TIKHONOV A N. On determining electrical characteristics of the deep layers of the Earth’s crust[J]. Dok1. Akad. Nauk. SSSR, 1950, 73(2): 295-297.

[2] CAGNIARD L.Basic theory of the magnetotelluric method of geophysical prospecting[J].Geophysics,1953,18(3):605-635.

[3] 朱威, 范翠松, 姚大为, 等. 矿集区大地电磁噪声场源分析及噪声特点[J]. 物探与化探, 2011, 35(5):658-662. ZHU W, FAN C S , YAO D W,et al. Noise source analysis and noise characteristic study of MT in an ore concentration area[J]. Geophysical and Geochemical Exploration, 2011, 35(5): 658-662. (In Chinese)

[4] 王书明, 王家映. 大地电磁信号统计特征分析[J]. 地震学报, 2004, 26(6): 669-674. WANG S M, WANG J Y. Analysis on statistic characteristics of magnetotelluric signal[J]. Acta Seismologica Sinca, 2004, 26(6): 669-674. (In Chinese)

[5] 王书明, 王家映. 大地电磁信号非最小相位性的讨论[J]. 地球物理学进展, 2004, 19(2): 216-221. WANG S M, WANG J Y. Discussion on the non-minimum phase of magnetotelluric signals[J]. Progress in Geophysics, 2004, 19(2): 216-221. (In Chinese)

[6] 汤井田, 化希瑞, 曹哲民, 等. Hilbert-Huang变换与大地电磁噪声压制[J]. 地球物理学报, 2008, 51(2):603-610. TANG J T, HUA X R, CAO Z M, et al. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J]. Chinese Journal Geophysics, 2008, 51(2):603-610. (In Chinese)

[7] 景建恩, 魏文博, 陈海燕, 等. 基于广义S变换的大地电磁测深数据处理[J]. 地球物理学报, 2012, 55 (12): 4015-4022. JING J E, WEI W B, CHEN H Y, et al. Magnetotelluric sounding data processing based on generalized S transformation[J]. Chinese Journal Geophysics，2012, 55(12): 4015-4022. (In Chinese)

[8] KAPPLE K N. A data variance technique for automated despiking of magnetotelluric data with a remote reference[J]. Geophysical Prospecting, 2012, 60(1): 179-191.

[9] 汤井田, 李晋, 肖晓, 等. 数学形态滤波与大地电磁噪声压制[J]. 地球物理学报, 2012, 55(5): 1784-1793. TANG J T, LI J, XIAO X, et al. Mathematical morphology filtering and noise suppression of magnetotelluric sounding data[J]. Chinese Journal Geophysics, 2012, 55(5): 1784-1793. (In Chinese)

[10]王辉, 魏文博, 金胜, 等. 基于同步大地电磁时间序列依赖关系的噪声处理[J]. 地球物理学报, 2014, 57(2): 531-545. WANG H, WEI W B, JIN S, et al. Removal of magnetotelluric noise based on synchronous time series relationship[J]. Chinese Journal Geophysics, 2014, 57(2): 531-545. (In Chinese)

[11]唐贵基，庞彬，刘尚坤. 基于奇异差分谱和平稳子空间分析的滚动轴承故障诊断[J]. 振动与冲击, 2015, 34 (11): 83-87. TANG G J, PANG B, LIU S K. Fault diagnosis of rolling bearings based on difference spectrum of singular value and stationary subspace analysis[J]. Journal of Vibration and Shock, 2015, 34 (11): 83-87. (In Chinese)

[12]陆文凯, 丁文龙, 张善文, 等. 基于信号子空间分解的三维地震资料高分辨率处理方法[J]. 地球物理学报, 2009, 52 (8): 2174-2181. LU W K, DING W L, ZHANG S W, et al. A high resolution processing technique for 3_D seismic data based on signal subspace decomposition[J]. Chinese Journal Geophysics, 2009, 52(8): 2174-2181. (In Chinese)

[13]DENDRINOS M, BAKAMIDIS S, CARAYANNI G. Speech enhancement from noise: A regenerative approach[J]. Speech Communication, 1991, 10(2): 45-57.

[14]EPHRAIM Y, VAN TREES H L. A signal subspace approach for speech enhancement[J]. IEEE Transactions on Speech and Audio Processing, 1995, 3(4): 251-266.

[15]JENSEN S H, HANSEN P C, HANSEN S D, et al. Reduction of broad-band noise in speech by truncated QSVD[J]. IEEE Transactions on Speech and Audio Processing, 1995, 3(6): 439-448.

[16]李晋, 汤井田, 王玲, 等. 基于信号子空间增强和端点检测的大地电磁噪声压制[J]. 物理学报, 2014,63(1): 019101. LI J, TANG J T, WANG L, et al. Noise suppression for magnetotelluric sounding data based on signal subspace enhancement and endpoint detection[J]. Acta Physica Sinica, 2014, 63(1): 019101. (In Chinese)

[17]汤井田, 刘子杰, 刘峰屹, 等. 音频大地电磁法强干扰压制试验研究[J]. 地球物理学报, 2015, 58(12): 4636-4647. TANG J T, LIU Z J, LIU F Y, et al. The denoising of the audio magnetotelluric data set with strong interferences[J]. Chinese Journal Geophysics, 2015, 58(12): 4636-4647. (In Chinese)

[18]代小威. 基于V5-2000格式MT时间序列处理与功率谱估计及软件开发[D]. 武汉：长江大学, 2015. DAI X W. Processing and power spectrum estimation of MT time serizes in V5-2000 format and software development[D]. Wuhan: Yangtze University, 2015. (In Chinese)

Application research of subspace enhancement for magnetotelluric signal to noise separation

YAN Huan1, LI Jin1,2, LI Chunyang1

(1.Institute of Physics and Information Science, Hunan Normal University, Changsha 410081, China; 2.School of Geosciences and Info-Physics, Key Laboratory of Metallogenic Prediction of Non-Ferrous Metals and Geological Environment Monitor, Ministry of Education, Central South University, Changsha 410083, China)

In order to improve the quality of magnetotelluric sounding data and retain the useful information of low frequency band for magnetotelluric, suppressing the large-scale strong interference is imperative. Aimed at typical strong interference types of ore concentration area, we introduced the subspace enhancement method to study signal to noise separation of magnetotelluric data and strong interference. Through simulated large-scale square wave and charge and discharge triangle wave by computer, we discussed the optimal range of frame and criterion threshold of eigenvalue in subspace enhancement, and the concrete steps of the algorithm are presented. Then, the measured magnetotelluric data of ore concentration area were processed by signal to noise separation. Experiment results show that the method is feasible, which can effectively eliminate the large-scale strong noise interference in the time domain and the near source interference of ore concentration area has been effectively suppressed. At the same time, the apparent resistivity curve is more smoothly and continuously. Moreover, the quality of magnetotelluric sounding data is improved significantly in the low frequency band.

magnetotelluric sounding data; signal to noise separation; subspace enhancement; strong interference

2016-06-14 改回日期：2016-12-15

国家自然科学基金资助项目(41404111)；湖南省自然科学基金资助项目(2015JJ3088)；中国博士后科学基金资助项目(2015M570687) ;湖南省研究行科研创新项目(CX2017B224)

燕欢(1991-)，女，硕士，主要从事矿集区大地电磁强干扰压制研究，E-mail：aviva163@126.com。

李晋(1981-)，男，博士后，副教授，硕士生导师，主要从事信号处理及电磁勘探研究， E-mail：geologylj@163.com。

1001-1749(2017)03-0346-08

P 631.2

A

10.3969/j.issn.1001-1749.2017.03.08