地球物理学报  2017, Vol. 60 Issue (9): 3642-3654   PDF    
基于压缩感知重构算法的大地电磁强干扰分离
汤井田1,2, 李广1,2, 肖晓1,2 , 李晋1,2,3, 周聪1,2, 朱会杰4     
1. 中南大学地球科学与信息物理学院, 长沙 410083;
2. 有色金属成矿预测与地质环境监测教育部重点实验室(中南大学), 长沙 410083;
3. 湖南师范大学物理与信息科学学院, 长沙 410081;
4. 总装工程兵科研一所, 江苏无锡 214035
摘要:为压制大地电磁信号中的强人文干扰,提出一种基于压缩感知重构算法的大地电磁信号去噪方法.通过构建与常见典型强干扰相匹配而对有用信号不敏感的冗余字典原子,利用改进的正交匹配追踪算法,分离出大地电磁信号中的强干扰成分.为了验证所述方法的强干扰分离效果,首先通过在实测大地电磁信号中加入理想的强干扰信号进行了仿真分离实验,然后从大量实测数据中选取三种含有不同类型强干扰的时间域片段,用所述方法对实测数据中的强干扰进行分离,最后将所述方法应用于青海试验点以及庐枞矿集区某测点实测数据的综合处理.仿真实验结果表明,该方法在分离出强干扰的同时,能够较好地保留有用信号.实测数据处理结果表明,该方法能够有效压制强干扰,改善强干扰区大地电磁数据的质量.
关键词: 大地电磁信号处理      去噪      正交匹配追踪      压缩感知      冗余字典      形态滤波     
Strong noise separation for magnetotelluric data based on a signal reconstruction algorithm of compressive sensing
TANG Jing-Tian1,2, LI Guang1,2, XIAO Xiao1,2, LI Jin1,2,3, ZHOU Cong1,2, ZHU Hui-Jie4     
1. Institute of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring(Central South University), Ministry of Education, Changsha 410083, China;
3. Institute of Physics and Information Science, Hunan Normal University, Changsha 410081, China;
4. The First Engineering Scientific Research Institute of General Armaments Department, Jiangsu Wuxi 214035, China
Abstract: To suppress strong noise in raw magnetotelluric (MT) data, we propose a new time-series denoising method based on a signal reconstruction algorithm of compressive sensing. A redundant dictionary which matches with strong noise but insensitive to useful MT signal is designed, then the strong noise is separated from raw MT data by using this dictionary and an improved orthogonal matching pursuit algorithm. In order to verify the effect of the proposed method, firstly, a strong noise separation simulation experiment is carried out by adding ideal strong noise into the measured MT signal. Secondly, three time domain data segments with different types of typical strong noise selected from a large number of measured data are used for strong noise separation test. Finally, the proposed method is applied to measured data collected in the Qaidam Basin and Lu-Zong ore-concentration area. Simulation results show that the proposed method can effectively separate the strong noise from raw data while keep the useful part. The results of the applications to measured data show that our proposed method is an effective method to suppress strong noise and therefore improve the quality of MT data collected in areas with strong interference.
Key words: Magnetotelluric signal processing      Denoising      Orthogonal matching pursuit      Compressive sensing      Redundant dictionary      Morphology filtering     
1 引言

大地电磁(Magnetotelluric, MT)法利用天然电磁场作为场源,是探测深度最大的电磁勘探方法.由于利用天然电磁场作为场源,大地电磁法有两个固有的缺点,一是地表观测信号微弱,二是场源激发随机(陈乐寿等, 1990; 何继善,2010).这两个缺点使得大地电磁法极易受到其他噪声的干扰,因此在大地电磁法的发展过程中,强干扰压制与分离一直是地球物理学者所关注的问题(肖晓等, 2011; 汤井田等, 2015).

随着人文干扰的增强以及探测深度的增加,现有的大地电磁数据处理手段仍然存在一定的局限性.例如互功率谱法对相关噪声抑制能力有限(Goubau et al., 1978);远参考法存在处理后误差棒明显增大的现象, 且选择参考点的位置也是一个难点(Gamble et al., 1979); Robust法对输入端噪声无能为力,且近源干扰压制效果不佳(Egbert, 1997; 柳建新等, 2003);小波变换能够有效压制局部相关噪声,但过分依赖小波基函数的选取,有时随着尺度增大,相应正交基函数的频谱局部性变差,使其对大地电磁信号更精细的分解受到限制(Trad and Travassos, 2000; 汤井田等, 2012a);Hilbert-Huang变换法能有效抑制工频干扰,但无法揭示每时段的频率特性和能量差异所具有的细微变化,且该算法效率低(汤井田等, 2008, 2009);数学形态滤波在MT信号处理中取得良好的效果,但是形态滤波主要作用于大尺度的低频干扰,且去噪的同时会去除部分有用信号(汤井田等,2012a李晋等,2014; 汤井田等,2015);近年来有学者采用匹配追踪算法较好地压制了大地电磁信号中谐波干扰及尖峰干扰,但对方波、阶跃、充放电等类型的强干扰则没有取得好的效果(郑思伟, 2015).总的来说,现有方法对于近源干扰的处理效果不佳.

从复杂信号中提取信息,通常需要对信号进行分解,在冗余字典稀疏分解思想出现以前,信号的分解都是基于正交基展开,但基于正交基的分解灵活性不够,在很多情况下无法达到令人满意的效果(Mallat et al., 1993; Qian and Chen, 1994; 王天荆等, 2011; 武国宁等, 2012; 蔡涵鹏等, 2013).为了更好地与信号结构匹配,Mallat和Zhang (1993)提出了一种叫做匹配追踪(Matching Pursuit, MP)的算法.匹配追踪算法能够将任何信号表示成从冗余字典中选择的原子组成的线性膨胀,由于字典的选择不受限制,因此匹配追踪是一种非常灵活的信号表示方法.在MP的基础上,Pati等(1993)提出了正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法,OMP比MP多了正交化的步骤,能够保证算法在有限的迭代次数中收敛.

2004年Candès等人发现可以通过小部分傅里叶变换系数精确重构原始信号(Donoho, 2006; Candes and Tao, 2006; Candès et al., 2006),随后正式提出了压缩感知的概念(Donoho,2006; Candès, 2006),压缩感知突破了奈奎斯特采样定理对采样速率的约束,将采样与压缩一步完成,数据量大大减少,可有效降低传感单元与存储空间等物理资源的开销(Candè et al., 2008;石光明等,2009; 刘亚新等,2010),因此压缩感知提出后在各领域引起了强烈反响(李树涛和魏丹, 2009; 唐刚和杨慧珠, 2010).由于重构算法是压缩感知的关键技术之一,国内外学者针对压缩感知信号的重构,在OMP的基础上提出了许多改进算法(杨真真等, 2013),如压缩采样匹配追踪(Compressive Sampling Matching Pursuit, CoSaMP)算法(Needell and Tropp, 2009)、子空间追踪(Subspace Pursuit, SP)算法(Dai and Milenkovic, 2009)、正则化正交匹配追踪(Regularized Orthogonal Matching Pursuit, ROMP)算法(Needell and Vershynin, 2010)、回溯自适应正交匹配追踪(Backtracking-based Adaptive Orthogonal Matching Pursuit, BAOMP)算法(Huang and Makur, 2011)、分段正交匹配追踪(Stagewise Orthogonal Matching Pursuit, StOMP)算法(Donoho et al., 2012)、多匹配正交追踪(Multi-Matching Orthogonal Pursuit, MMOP)算法(李凤莲等, 2015)等等,这些算法与最初的匹配追踪算法相比在重构精度和效率上有明显的提高(Tropp and Gilbert, 2007; Davenport and Wakin, 2010; 刘亚新等, 2010).

压缩感知重构算法已经广泛应用于图像与信号去噪,如语音信号去噪(李凤莲等, 2015)、地震资料去噪(Hennenfent et al., 2008; Tang and Ma, 2011; 武国宁等,2012).在上述背景下,本文基于压缩感知信号重构原理,尝试构造与各类强干扰相匹配的冗余字典原子,利用改进的正交匹配追踪(Improved Orthogonal Matching Pursuit, IOMP)算法,对含有强干扰的大地电磁时间序列进行信噪分离,以期改善大地电磁数据质量.

2 方法原理 2.1 信噪分离依据

从观测到的时域信号中准确识别出天然的大地电磁信号具有较大的难度,但强干扰信号常具有幅值大、波形规则、极化方向稳定等特征,因此可利用这些特征对强干扰信号进行识别(Weckmann et al., 2005; 李晋等,2014; 汤井田等,2015).根据这一机理,我们把强干扰作为目标信号,用构造的原子与强干扰进行匹配,将匹配部分从观测到的大地电磁信号中分离,留下残差作为有用的大地电磁信号.冗余字典原子的构造决定了信号重构的效果,构造的原子与目标信号越相似,信号分离效果就越好(Donoho et al., 2012; 朱会杰等, 2015).对于谐波干扰,可以用正弦原子和余弦原子进行匹配,对于方波、尖峰等干扰,构造的原子如下:

(1)

在时间段[τ1, τ2],原子g (τ1, τ2)的幅值为A,其余时间段幅值为0,那么使用若干个原子就可以构成一个方波信号.阶跃信号可以看成是不同幅值的方波信号构成,脉冲信号可以看成是周期很短的方波信号,因此构造的方波原子可以很好地与阶跃信号和脉冲信号匹配.当时间段[τ1, τ2]较窄时,通过有限个数的原子可以近似构成类尖峰以及其他复杂的强干扰信号.

2.2 IOMP算法流程

IOMP算法是根据SP和OMP改进而来(朱会杰等,2015).OMP算法是一种将被处理信号一次一次逐步分解的过程.首先,在原子库中选择与被处理信号最匹配的原子,然后将匹配的成分从原信号中减去得到残差; 在之后的迭代中,每次选择与残差最匹配的原子后,先将所有选中的原子正交化,之后重新在原信号上进行投影,提取匹配部分并留下残差; 依次反复迭代,直到残余信号的能量小于设定的阈值或达到最大迭代次数(Tropp and Gilbert, 2007Cai and Wang, 2011).

OMP每次仅选一个原子,并且选定的原子不会在以后的迭代中更换,这样将导致如果在某次迭代中选择了一个错误的原子,将需要多次的迭代才能填补该误差.SP是在OMP基础上结合回溯思想的一种改进算法,每次选择K个原子作为候选,并从候选原子中选择K个原子加入已选原子集合,SP可以将错误的原子剔除并替换为正确的原子,因此比OMP有更好的重构效果,但SP需要事先知道目标信号的稀疏度,这在实际应用中很难满足.结合OMP与SP两者优点提出的改进型正交匹配追踪算法,每次选择多个原子加入候选原子集,并将整个候选支撑集投影到原信号上,选择内积最大的n个原子作为支撑集,对全部已选原子进行更新,因此对上次已选原子没有依赖.且IOMP依据信号与残差的能量比Eg能够自适应地确定信号的稀疏度.IOMP算法的步骤如下:

D={gγ}γΓ为冗余字典的原子库,γ为原子gγ的序号,且‖gγ2=1,Γ为序号γ组成的集合,y为观测信号,rn为残差,Λn为候选原子集,ψn为已选原子集,则IOMP步骤如下:

输入:观测信号y,最小原子宽度γmin,最大原子宽度γmax,残差能量比Eg,原子个数K,迭代次数n

输出:重构信号,残差rn.

(1) 初始化:观测信号y,迭代次数n=0,初始残差r0y,已选原子集合ψ0为空集,候选原子集合Λ0为空集.

(2) 迭代次数n=n+1,依次从冗余字典中找出与当前残差最匹配的K个原子放入候选原子集合中.

(2)

(3) 将候选原子集Λn的所有原子投影到观测信号y上,选出n个具有最大内积的原子作为已选原子集合ψn.

(3)

(4) 更新重构信号与残差:

(4)

(5)

(5) 依据噪声阈值标准和残差能量判断是否已经达到原信号与残差的临界点.如果已经达到,则保持稀疏度为n,继续迭代若干次,若临界点不再发生变化,则输出重构信号;若临界点发生变化,或者第n次迭代没有达到临界点,则转到步骤(2) 继续执行.

2.3 数据处理流程

MT测点的数据量很大,通常一个测道的数据为几十万到上百万个采样点,如果采用IOMP算法一次性全部处理,由于计算量大(Donoho et al., 2012),耗时过长(后文有具体分析),且实测数据通常并非所有时段均含有强干扰,因此实际操作时,仅选取有明显强干扰的数据段进行处理.本文所述大地电磁数据处理的流程如图 1所示,首先根据大地电磁强干扰信号的特征将全时段数据分为若干个强干扰段和非强干扰段,非强干扰段不作处理予以保留,强干扰段通过IOMP算法分离出强干扰后与非强干扰段一起根据原采样顺序进行拼接,之后用于视电阻率和相位计算或其他用途.

图 1 大地电磁数据处理流程 Fig. 1 Flow chart of MT data processing
3 强干扰分离仿真

为了测试算法对大地电磁信号中强干扰的分离效果,以青海某地实测无明显干扰的大地电磁原始信号为基础,根据大地电磁信号常见的强干扰类型(汤井田等, 2012b),人为添加仿真产生的强干扰信号得到含噪信号,然后通过IOMP算法从含噪信号中分离出添加的强干扰,含噪信号减去强干扰得到大地电磁重构信号.

3.1 仿真结果评价标准

实验时,通过计算原始信号与重构信号的重构误差E和曲线相似度NCC两个指标来定量评价强干扰分离效果,并与形态滤波法的强干扰分离效果进行对比.重构误差E定义为:

(6)

曲线相似度NCC定义为:

(7)

式(6) 和(7) 中y (x)表示原始信号,r (x)表示重构信号,x为采样点的序号,N为信号长度.重构误差E为0时表示重构信号与原始信号完全一样,因此E越接近0表示强干扰分离效果越好.NCC值为1时表示两个信号完全一样,为0时表示两个信号正交,为-1时表示两个信号绝对值相同而符号相反,因此NCC值越接近1表示强干扰分离效果越好(汤井田等, 2014).

3.2 强干扰分离仿真结果

图 2为在实测数据中加入谐波干扰后用本文方法进行强干扰分离的结果,其中图 2a为实测原始信号,图 2b为含噪信号,图 2c为形态滤波法分离出的强干扰信号(结构元素为4点直线型),图 2d为本文方法分离出的强干扰信号,图 2e为形态滤波法重构的大地电磁信号,图 2f为本文方法重构的大地电磁信号.图 2g2l为与图 2a2f各信号一一对应的频域效果.表 1为该仿真结果的定量评价.

表 1 谐波干扰分离效果定量评价 Table 1 Quantitative evaluation of harmonic wave separation
图 2 谐波干扰分离仿真结果 Fig. 2 Simulation results of harmonic wave separation

图 2中时域效果和频域效果均表明,两种方法均能够有效分离出谐波干扰,与原始信号相比,形态滤波法处理时损失了部分有用信号,而本文方法更好地保留了有用信号.由表 1可知,本文方法的重构误差更小,曲线相似度更高.

图 3为在实测数据中加入伪随机方波和尖峰干扰后用本文方法进行强干扰分离的结果,图中各曲线的含义与图 2类似(形态滤波结构元素为2点抛物线),表 2为该仿真结果的定量评价.如图 3所示,形态滤波法能够有效分离出含噪信号中的方波干扰,但对于尖峰干扰去除不彻底,且提取的方波干扰轮廓不光滑,损失了部分有用信号;本文方法能分离出各种宽度的方波干扰以及尖峰干扰,且噪声轮廓光滑.由表 2可知,本文方法重构的大地电磁信号与原始信号相似度为0.9660,重构误差为0.2651,与形态滤波法相比具有明显的优势.

图 3 方波和尖峰干扰分离仿真结果 Fig. 3 Simulation results of square waves and spike interference separation
表 2 方波和尖峰干扰分离效果定量评价 Table 2 Quantitative evaluation of square waves and spike interference separation

图 3还说明,尽管采用的原子是单一的方波原子,但是对于方波干扰和尖峰干扰均能在保留有用信号的前提下取得较好的分离效果,因此方波原子具有较强的适应性.

对于方波原子构成的冗余字典,其冗余程度由方波原子的最小宽度γmin和最大宽度γmax决定(原子的幅值由算法根据实际干扰信号自适应匹配获得),实际MT信号中含有大量的脉冲干扰,故γmin通常设置为最低值1,因此大地电磁强干扰分离时,字典的冗余度仅由γmax决定,γmax的值越大,意味着需要从更多的原子中追踪最匹配的原子.为分析冗余度对算法性能的影响,通过改变γmax的值、保持其他参数不变,对图 3b所示含噪信号(γmax为155) 进行强干扰分离.表 3所示为不同冗余度下仿真结果统计,其中T表示算法所需的运算时间,取10次重复实验的平均值(本文所有数据处理采用同一台计算机完成,其CPU为i5-4590,主频3.3 GHz,RAM为8 GB,操作系统为64位).

表 3 不同冗余度下仿真结果 Table 3 Simulation results under different redundancies

表 3统计的结果可知,随着γmax的值降低,运算耗时急剧下降且重构效果有所改善;当γmax等于干扰信号最大宽度时,重构效果最佳且耗时较少;随着γmax的值继续降低,曲线相似度迅速下降,重构误差迅速增加,而重构时间减少有限.因此在大地电磁信号处理时,为提高效率,应该适当降低原子的最大宽度,即降低字典的冗余度,为保证重构效果, γmax的值不宜低至实际干扰信号的最大宽度以下.

同样对于图 3b所示含噪信号,保持最大原子宽度等参数不变,仅改变设定的原子个数K(即估计的稀疏度),得到不同预设原子个数时强干扰分离仿真效果如图 4所示.由图 4可知,随着K值的增大,运算耗时略有上升;曲线相似度随着K值的增大而迅速上升,当K等于某一个数值时达到最高,随后缓慢下降;重构误差随着K值的增大而显著下降,当K等于某一个数值时重构误差降至最低,随后略有上升.实际上,曲线相似度最高时对应的K即为目标信号在设定字典下的实际稀疏度,此时重构误差也最低.综上可知,当设定的原子个数小于实际稀疏度时,重构效果显著变差,而大于实际稀疏度时重构效果不会明显变差,因此实际处理时,设定的原子个数应该遵循“宁滥勿缺”的原则,略大于实际稀疏度.

图 4 不同预设原子个数时仿真效果 Fig. 4 Simulation results under different preset atomic numbers
4 实测数据处理 4.1 不同干扰类型分类处理

由于理想的干扰信号与实际干扰信号存在一定差异,因此从实测数据中选取了三种典型的受强干扰影响的数据段,用所述方法对其中的强干扰进行分离.实际处理时,对于谐波等幅值变化较平缓的干扰,采用正弦原子、余弦原子构成的冗余字典,对于方波、尖峰以及其他幅值突变的强干扰,则采用方波原子构成的冗余字典.首先估算强干扰信号的最大宽度,并通过试处理估计信号的稀疏度,然后根据仿真所得规律,将原子的最大宽度设置为略大于干扰信号的最大宽度,原子个数设置为略大于估计的稀疏度.

4.1.1 谐波干扰分离

图 5所示为谐波干扰分离时域效果(a—c)和频域效果(d—f),该部分数据来源于庐枞矿集区某测点Ey道,采样频率为2560 Hz.由图 5可知,原始信号以及噪声轮廓均呈现出明显的周期性,而IOMP重构信号随机并且幅度较小,呈现出明显的天然大地电磁信号特征.频域效果表明,经过处理后,谐波干扰被显著压制,同时没有引入新的干扰.

图 5 谐波干扰分离效果 Fig. 5 Results of harmonic interference separation
4.1.2 方波干扰分离

图 6所示为方波干扰分离时域效果(a—c)和频域效果(d—f),该部分数据来源于庐枞矿集区某测点的Ex道,采样频率为15 Hz.由时域效果可知,本文方法分离出的方波干扰轮廓清晰且光滑,保留了信号的细节部分,重构信号随机并且幅度较小,呈现出明显的天然大地电磁信号特征,由频域效果可知,经过处理后,强干扰信号显著衰减.

图 6 方波干扰分离效果 Fig. 6 Results of square wave interference separation
4.1.3 尖峰干扰分离

图 7所示为尖峰干扰分离时域效果(a—c)和频域效果(d—f),该部分数据来源于青海某测点的Ex道,采样频率为2400 Hz.由时域效果可知,本文方法有效地分离出了尖峰干扰,无强干扰时段噪声轮廓光滑且幅值为0,重构信号随机并且幅度较小,呈现出明显的天然大地电磁信号特征,由频域效果可知,原始信号频谱和噪声频谱均呈现出一定的规律性,而重构信号频谱呈现出天然大地电磁信号的随机特征.

图 7 尖峰干扰分离效果 Fig. 7 Results of spike interference separation
4.2 青海试验点数据处理

为研究大地电磁信号的强干扰分离方法,课题组于2012年在青海柴达木盆地进行了野外试验(李晋等,2014).试验思路是在大地电磁测点附近布置广域电磁发送机,以广域电磁发送机发送的伪随机强干扰信号作为干扰源,以未受到广域电磁发射源干扰时的数据作为对比标准.本文利用青海相关试验获得的数据进行对比.某试验点编号为QH4015,采集时间共约19 h (全时段).其中,在开始约1.5 h (干扰段)内,距该测点2 km处的广域电磁发送机正在向地下注入电流为80 A的伪随机信号,使MT数据采集的同时也接收广域电磁发射源的干扰;而随后的近17.5 h (无干扰段)内,广域电磁发射机停止工作,MT采集的数据未受其影响.由于测点所在地基本无人文干扰且地形地貌情况简单,因此在广域电磁发送机停止工作时采集的数据是无明显干扰的数据.

4.2.1 时间域效果分析

图 8为有广域电磁发射源干扰时实测原始信号时间域片段,采样频率为15 Hz.图 9图 8所示时间域片段经过IOMP法重构后的信号.分析图 8图 9可知,处理前时间域信号呈现出周期性与规则性,有用的大地电磁信号几乎被完全湮没,而经过IOMP法处理后时间域信号呈现出天然大地电磁信号的随机性与无规律性,且幅度较小.

图 8 实测数据时间域片段 Fig. 8 Time series segments of measured data
图 9 实测数据处理后时间域片段 Fig. 9 Time series segments after processing of measured data

表 4图 8所示Ey道数据在不同长度下的实际处理耗时分析,保持原子的最大宽度等参数不变,设定相同的残差能量比,仅改变数据长度和设定的原子个数,且设定的原子个数与信号长度的比值相同.由表 4可知,随着信号长度的增加,处理耗时急剧增加,因此实际处理时,通常通过分段处理以提高效率.

表 4 不同长度数据处理耗时分析 Table 4 Time-consumption analysis of different lengths of data
4.2.2 视电阻率与相位分析

图 10所示为该试验点处理前后视电阻率和相位曲线,图中W表示无干扰段的视电阻率或相位曲线,Q表示处理前全时段的视电阻率或相位曲线,H表示处理后全时段的视电阻率或相位曲线.表 5为以无干扰段作为评价标准计算得到的曲线相似度与重构误差结果,其中Q-W表示处理前相对于无干扰,H-W表示处理后相对于无干扰,NCC_Rxy表示XY方向上视电阻率的曲线相似度,E_Rxy表示XY方向上视电阻率的重构误差,其余变量的含义以此类推.

表 5 测点QH4015的视电阻率和相位曲线定量评价 Table 5 Quantitative evaluation of apparent resistivity and phase curves of site QH4015
图 10 测点QH4015视电阻率和相位曲线 Fig. 10 Apparent resistivity and phase curves of site QH4015

分析图 10可知,处理前全时段的视电阻率和相位曲线整体连续性较差,不光滑.XY方向的视电阻率和相位曲线在0.3 Hz到0.05 Hz之间跳变剧烈.YX方向的视电阻率在0.5 Hz到0.01 Hz之间呈45°上升,低于0.01 Hz时剧烈下掉,相位曲线在0.5 Hz到0.01 Hz之间大多接近-180°.综上可知,处理前,该测点全时段的数据表现为明显的近源干扰.

经过本文方法处理后,全时段的视电阻率和相位曲线在整个频段内都具有良好的连续性与光滑度,0.5 Hz到0.01 Hz之间YX方向视电阻率的45°上升趋势得到明显好转,相位回归到正常,无明显的近源干扰特征.由表 5可知,处理后数据的视电阻率和相位曲线与处理前相比,曲线相似度显著提高,重构误差显著降低.

4.2.3 极化方向分析

电磁场的极化方向是评价大地电磁信号受干扰程度的重要指标之一.当大地电磁信号中没有强干扰时,其极化方向是随机的,而受到强干扰时,由于干扰信号明显强于天然电磁场,导致极化方向较为集中,表现出有规律性(Weckmann et al., 2005; 汤井田等,2015).图 11所示为干扰段处理前与处理后低频部分(TS5) 3 Hz电磁场的极化方向,其中图 11a为处理前电场极化方向, 图 11b为处理前磁场极化方向, 图 11c为处理后电场极化方向, 图 11d为处理后磁场极化方向.从图 11可以看出,处理前干扰段电磁场极化方向高度集中,呈明显的规律性,说明该部分数据受到了非常强烈的干扰.经过IOMP法分离出强干扰后,不论是电场极化方向还是磁场极化方向,都表现出明显的无规律性.从而说明经过处理后,强干扰被极大衰减.

图 11 干扰段处理前后电磁场极化方向 Fig. 11 Polarization directions of disturbed section before and after processing
4.3 庐枞矿集区数据处理

庐枞矿集区C5132测点总采集时间大约为21 h (全时段),通过分析采集到的时间域数据可知,前11 h (干扰段)受干扰严重,而在后10 h (无干扰段)无明显强干扰.图 12所示为C5132测点视电阻率和相位曲线,图中各曲线含义同图 10.

分析图 12可知,经过IOMP法处理后,全时段的视电阻率和相位曲线的连续性和光滑度得到明显改善,电阻率和相位曲线相对于处理前更为接近无明显强干扰时数据段的视电阻率和相位曲线.

图 12 测点C5132视电阻率和相位曲线 Fig. 12 Apparent resistivity and phase curves of site C5132
5 结论

本文以观测信号中的强干扰为目标信号,将稀疏分解引入大地电磁信号处理,构建了与常见典型强干扰相匹配而对有用信号不敏感的冗余字典原子,提出一种基于压缩感知重构算法的大地电磁强干扰分离方法,以改善强干扰区大地电磁数据质量.通过系统的仿真实验及实测数据处理,得出以下五点结论:第一,所述方法能够显著提高强干扰区大地电磁数据的质量,改善视电阻率和相位曲线,压制近源干扰;第二,所述方法的主要优点在于仅对强干扰信号敏感而对有用信号不敏感,能够在较好地保留有用信号的前提下,有效分离出不同频率、不同类型的强干扰信号;第三,为减少计算量,提高处理效率,应该适当降低原子的最大宽度以减少字典的冗余度,但不宜低于干扰信号的最大宽度;第四,实测大地电磁数据中,干扰类型众多,预先设定的冗余字典在处理某些复杂干扰时会产生一定的误差;第五,压缩感知重构算法具有计算量大、重构耗时的缺点,海量数据处理时可以通过并行计算或者使用高性能计算机等方式提高效率.

参考文献
Cai H P, He Z H, Gao G, et al. 2013. Seismic data matching pursuit using hybrid optimization algorithm and its application. Journal of Central South University (Science and Technology Edition), 44(2): 687-694.
Cai T T, Wang L. 2011. Orthogonal matching pursuit for sparse signal recovery with noise. IEEE Transactions on Information Theory, 57(7): 4680-4688. DOI:10.1109/TIT.2011.2146090
Candès E J. 2006. Compressive sampling.//Proceedings of the International Congress of Mathematicians. Madrid:EMS Publishing House, 3:1433-1452.
Candès E J, Romberg J, Tao T. 2006. Robust uncertainty principles:Exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2): 489-509. DOI:10.1109/TIT.2005.862083
Candè E J, Wakin M B. 2008. An introduction to compressive sampling. IEEE Signal Processing Magazine, 25(2): 21-30. DOI:10.1109/MSP.2007.914731
Candes E J, Tao T. 2006. Near-optimal signal recovery from random projections:Universal encoding strategies?. IEEE Transactions on Information Theory, 52(12): 5406-5425. DOI:10.1109/TIT.2006.885507
Chen L S, Wang G E. 1990. Magnetotelluric Sounding Method. Beijing: Geological Publishing House.
Dai W, Milenkovic O. 2009. Subspace pursuit for compressive sensing signal reconstruction. IEEE Transactions on Information Theory, 55(5): 2230-2249. DOI:10.1109/TIT.2009.2016006
Davenport M A, Wakin M B. 2010. Analysis of orthogonal matching pursuit using the restricted isometry property. IEEE Transactions on Information Theory, 56(9): 4395-4401. DOI:10.1109/TIT.2010.2054653
Donoho D L. 2006. Compressed sensing. IEEE Transactions on Information Theory, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
Donoho D L, Tsaig Y, Drori I, et al. 2012. Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit. IEEE Transactions on Information Theory, 58(2): 1094-1121. DOI:10.1109/TIT.2011.2173241
Egbert G D. 1997. Robust multiple-station magnetotelluric data processing. Geophysical Journal International, 130(2): 475-496. DOI:10.1111/gji.1997.130.issue-2
Gamble T D, Goubau 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
He J S. 2010. Wide Field Electromagnetic Method and Pseudo-Random Signal Electrical Method. Beijing: Higher Education Press.
Hennenfent G, Herrmann F J. 2008. Simply denoise:Wavefield reconstruction via jittered undersampling. Geophysics, 73(3): V19-V28. DOI:10.1190/1.2841038
Huang H L, Makur A. 2011. Backtracking-based matching pursuit method for sparse signal reconstruction. IEEE Signal Processing Letters, 18(7): 391-394. DOI:10.1109/LSP.2011.2147313
Li F L, Chang J, Zhang X Y, et al. 2015. Reconstruction algorithm of blind sparse and its de-noising application in speech based on compressed sensing. Journal of Central South University (Science and Technology Edition), 46(1): 164-170.
Li J, Tang J T, Xiao X, et al. 2014. Magnetotelluric data processing based on combined generalized morphological filter. Journal of Central South University (Science and Technology Edition), 45(1): 173-185.
Li S T, Wei D. 2009. A survey on Compressive Sensing. Acta Automatica Sinica, 35(11): 1369-1377. DOI:10.3724/SP.J.1004.2009.01369
Liu J X, Yan J B, He J S. 2003. Robust estimation method of sea magnetotelluric impedance based on correlative coefficient. Chinese Journal of Geophysics, 46(2): 241-245. DOI:10.3321/j.issn:0001-5733.2003.02.018
Liu Y X, Zhao R Z, Hu S H, et al. 2010. Regularized adaptive matching pursuit algorithm for signal reconstruction based on compressive sensing. Journal of Electronics & Information Technology, 32(11): 2713-2717.
Mallat S G, Zhang Z F. 1993. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41(12): 3397-3415. DOI:10.1109/78.258082
Needell D, Tropp J A. 2009. CoSaMP:Iterative signal recovery from incomplete and inaccurate samples. Applied and Computational Harmonic Analysis, 26(3): 301-321. DOI:10.1016/j.acha.2008.07.002
Needell D, Vershynin R. 2010. Signal recovery from incomplete and inaccurate measurements via regularized orthogonal matching pursuit. IEEE Journal of Selected Topics in Signal Processing, 4(2): 310-316. DOI:10.1109/JSTSP.2010.2042412
Pati Y C, Rezaiifar R, Krishnaprasad P S. 1993. Orthogonal matching pursuits:Recursive function approximation with applications to wavelet decomposition. Proceedings of the 27th Asilomar Conference on Signals, Systems and Computers. Pacific Grove, CA, USA:IEEE.
Qian S E, Chen D P. 1994. Signal representation using adaptive normalized Gaussian functions. Signal Processing, 36(1): 1-11. DOI:10.1016/0165-1684(94)90174-0
Shi G M, Liu D H, Gao D H, et al. 2009. Advances in theory and application of compressed sensing. Acta Electronica Sinica, 37(5): 1070-1081.
Tang G, Yang H Z. 2010. Seismic data compression and reconstruction based on Poisson Disk sampling. Chinese Journal of Geophysics, 53(9): 2181-2188. DOI:10.3969/j.issn.0001-5733.2010.09.018
Tang G, Ma J W. 2011. Application of total-variation-based curvelet shrinkage for three-dimensional seismic data denoising. IEEE Geoscience and Remote Sensing Letters, 8(1): 103-107. DOI:10.1109/LGRS.2010.2052345
Tang J T, Hua X R, Cao Z M, et al. 2008. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data. Chinese Journal of Geophysics, 51(2): 603-610. DOI:10.3321/j.issn:0001-5733.2008.02.034
Tang J T, Cai J H, Ren Z Y, et al. 2009. Hilbert-Huang transform and time-frequency analysis of magnetotelluric signal. Journal of Central South University (Science and Technology Edition), 40(5): 1399-1405.
Tang J T, Li J, Xiao X, et al. 2012a. Mathematical morphology filtering and noise suppression of magnetotelluric sounding data. Chinese Journal of Geophysics, 55(5): 1784-1793. DOI:10.6038/j.issn.0001-5733.2012.05.036
Tang J T, Xu Z M, Xiao X, et al. 2012b. Effect rules of strong noise on magnetotelluric (MT) sounding in the Luzong ore cluster area. Chinese Journal of Geophysics, 55(12): 4147-4159. DOI:10.6038/j.issn.0001-5733.2012.12.027
Tang J T, Li H, Li J, et al. 2014. Top-Hat transformation and magnetotelluric sounding data strong interference separation of Lujiang-Zongyang ore concentration area. Journal of Jilin University (Earth Science Edition), 44(1): 336-343.
Tang J T, Liu Z J, Liu F Y, et al. 2015. The denoising of the audio magnetotelluric data set with strong interferences. Chinese Journal of Geophysics, 58(12): 4636-4647. DOI:10.6038/cjg20151225
Trad D O, Travassos J M. 2000. Wavelet filtering of magnetotelluric data. Geophysics, 65(2): 482-491. DOI:10.1190/1.1444742
Tropp J A, Gilbert A C. 2007. Signal recovery from random measurements via orthogonal matching pursuit. IEEE Transactions on Information Theory, 53(12): 4655-4666. DOI:10.1109/TIT.2007.909108
Wang T J, Zheng B Y, Yang Z. 2011. A speech signal sparse representation algorithm based on adaptive overcomplete dictionary. Journal of Electronics & Information Technology, 33(10): 2372-2377.
Weckmann U, Magunia A, Ritter O. 2005. Effective noise separation for magnetotelluric single site data processing using a frequency domain selection scheme. Geophysical Journal International, 161(3): 635-652. DOI:10.1111/gji.2005.161.issue-3
Wu G N, Cao S Y, Sun N. 2012. Matching pursuit method based on complex seismic traces and its application of hydrocarbon exploration. Chinese Journal of Geophysics, 55(6): 2027-2034. DOI:10.6038/j.issn.0001-5733.2012.06.024
Xiao X, Tang J T, Zhou C, et al. 2011. Magnetotelluric sounding in the Lujiang-Zongyang ore-district and preliminary study of electrical structure. Acta Geological Sinica, 85(5): 873-886.
Yang Z Z, Yang Z, Sun L H. 2013. A survey on orthogonal matching pursuit type algorithms for signal compression and reconstruction. Journal of Signal Processing, 29(4): 486-496.
Zheng S W. 2015. A study of matching pursuit on the suppression of magnetotelluric noise[Master's thesis] (in Chinese). Beijing:China University of Geosciences.
Zhu H J, Wang X Q, Rui T, et al. 2015. Implication of improved matching pursuit in de-noising for square wave. Journal of PLA University of Science and Technology (Natural Science Edition), 16(4): 305-309.
蔡涵鹏, 贺振华, 高刚, 等. 2013. 基于混合优化算法的地震数据匹配追踪分解. 中南大学学报(自然科学版), 44(2): 687–694.
陈乐寿, 王光锷. 1990. 大地电磁测深法. 北京: 地质出版社.
何继善. 2010. 广域电磁法和伪随机信号电法. 北京: 高等教育出版社.
李凤莲, 畅江, 张雪英, 等. 2015. 基于压缩感知的语音盲稀疏重构算法及其去噪应用. 中南大学学报(自然科学版), 46(1): 164–170. DOI:10.11817/j.issn.1672-7207.2015.01.023
李晋, 汤井田, 肖晓, 等. 2014. 基于组合广义形态滤波的大地电磁资料处理. 中南大学学报(自然科学版), 45(1): 173–185.
李树涛, 魏丹. 2009. 压缩传感综述. 自动化学报, 35(11): 1369–1377.
柳建新, 严家斌, 何继善, 等. 2003. 基于相关系数的海底大地电磁阻抗Robust估算方法. 地球物理学报, 46(2): 241–245. DOI:10.3321/j.issn:0001-5733.2003.02.018
刘亚新, 赵瑞珍, 胡绍海, 等. 2010. 用于压缩感知信号重建的正则化自适应匹配追踪算法. 电子与信息学报, 32(11): 2713–2717.
石光明, 刘丹华, 高大化, 等. 2009. 压缩感知理论及其研究进展. 电子学报, 37(5): 1070–1081.
唐刚, 杨慧珠. 2010. 基于泊松碟采样的地震数据压缩重建. 地球物理学报, 53(9): 2181–2188. DOI:10.3969/j.issn.0001-5733.2010.09.018
汤井田, 化希瑞, 曹哲民, 等. 2008. Hilbert-Huang变换与大地电磁噪声压制. 地球物理学报, 51(2): 603–610. DOI:10.3321/j.issn:0001-5733.2008.02.034
汤井田, 蔡剑华, 任政勇, 等. 2009. Hilbert-Huang变换与大地电磁信号的时频分析. 中南大学学报(自然科学版), 40(5): 1399–1405.
汤井田, 李晋, 肖晓, 等. 2012a. 数学形态滤波与大地电磁噪声压制. 地球物理学报, 55(5): 1784–1793. DOI:10.6038/j.issn.0001-5733.2012.05.36
汤井田, 徐志敏, 肖晓, 等. 2012b. 庐枞矿集区大地电磁测深强噪声的影响规律. 地球物理学报, 55(12): 4147–4159. DOI:10.6038/j.issn.0001-5733.2012.12.027
汤井田, 李灏, 李晋, 等. 2014. Top-Hat变换与庐枞矿集区大地电磁强干扰分离. 吉林大学学报(地球科学版), 44(1): 336–343.
汤井田, 刘子杰, 刘峰屹, 等. 2015. 音频大地电磁法强干扰压制试验研究. 地球物理学报, 58(12): 4636–4647. DOI:10.6038/cjg20151225
王天荆, 郑宝玉, 杨震. 2011. 基于自适应冗余字典的语音信号稀疏表示算法. 电子与信息学报, 33(10): 2372–2377.
武国宁, 曹思远, 孙娜. 2012. 基于复数道地震记录的匹配追踪算法及其在储层预测中的应用. 地球物理学报, 55(6): 2027–2034. DOI:10.6038/j.issn.0001-5733.2012.06.024
肖晓, 汤井田, 周聪, 等. 2011. 庐枞矿集区大地电磁探测及电性结构初探. 地质学报, 85(5): 873–886.
杨真真, 杨震, 孙林慧. 2013. 信号压缩重构的正交匹配追踪类算法综述. 信号处理, 29(4): 486–496.
郑思伟. 2015. 匹配追踪抑制大地电磁场噪声影响的研究[硕士论文]. 北京: 中国地质大学.
朱会杰, 王新晴, 芮挺, 等. 2015. 改进的匹配追踪在方波信号滤波中的应用. 解放军理工大学学报(自然科学版), 16(4): 305–309.