地球物理学报  2019, Vol. 62 Issue (10): 3934-3949   PDF    
多周期全波形激电抗干扰数据处理方法及在大规模探测中的应用分析
刘卫强1,2, 吕庆田1,2, 林品荣1, 陈儒军3     
1. 中国地质科学院地球物理地球化学勘查研究所, 自然资源部地球物理电磁法探测技术重点实验室, 河北廊坊 065000;
2. 中国地质科学院矿产资源研究所, 自然资源部成矿作用与资源评价重点实验室, 北京 100037;
3. 中南大学地球科学与信息物理学院, 有色金属成矿预测与地质环境监测重点实验室, 长沙 400071
摘要:激电法是金属矿产勘查中一种十分重要的电法勘探分支方法,但是各种电磁干扰的存在限制了激电法在大规模探测中的应用.近年来,国内外先后实现了三维分布式全波形激电探测仪器系统的研发和推广,全波形数据记录为激电信号的抗干扰处理提供了新的空间.本文针对多周期全波形采样的激电数据提出了一套基于统计分析的时间序列抗干扰数据处理方法,主要包括:经验模态分解用于分离低频趋势项干扰;相关分析用于消除突发性强噪声干扰;稳健统计用于多周期时间序列叠加;分段稳健平均和低频相对相位谱用于时/频域激电参数提取.将上述数据处理方法应用于由国产分布式电法系统实测的三维全波形激电数据,并与线性拟合、均值叠加等常处理方法进行对比,发现新方法可以有效识别和压制激电数据中的强噪声干扰,提高大供电极距和低频点探测时的激电数据质量,从而进一步推动激电法在深部矿产资源勘查中的应用.
关键词: 全波形激电      电磁噪声干扰      抗干扰处理      多参数提取     
Anti-interference processing of multi-period full-waveform induced polarization data and its application to large-scale exploration
LIU WeiQiang1,2, LÜ QingTian1,2, LIN PinRong1, CHEN RuJun3     
1. MNR Laboratory of Geophysical EM Probing Technologies, Institute of Geophysical and Geochemical Exploration, CAGS, Hebei Langfang 065000, China;
2. MNR Key Laboratory of Metallogeny and Mineral Assessment, Institute of Mineral Resources, CAGS, Beijing 100037, China;
3. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring, School of Geosciences and Info-Physics, Central South University, Changsha 400071, China
Abstract: Induced Polarization (IP) is one of the most important electrical methods in exploration of metallic ore deposits. However the existence of various electromagnetic interferences restricts its application in large-scale exploration. To solve this problem,3D disturbed full-waveform IP instrument systems have been developed in China and elsewhere in the world,of which data processing needs further improvement. In this paper,a complete anti-interference processing method based on statistical analysis for the multi-period full-waveform IP data is proposed. Its steps are as follows:Empirical mode decomposition is used to separate low-frequency trend term interference. Correlation analysis is made to eliminate the sudden strong noise interference. Robust statistics are performed to multi-period time series stack. The subsection robust average and low frequency relative phase spectrum are employed to extract the time/frequency domain IP parameters. This data processing method is applied to the three-dimensional full-waveform IP data acquired by a domestic distributed electrical system and compared with the common processing methods such as linear fitting and mean stacking. Results show that the new method can effectively identify and suppress the strong noise interference in IP data,improve the quality of IP data when the current electrode space is large and the frequency is low,and further promote the application of the IP method in exploration of mineral resources in the deep subsurface.
Keywords: Full-waveform induced polarization    Electromagnetic noise interference    Anti-interference processing    Multi-parameter extraction    
0 引言

激电法是矿产资源勘查中最常用的方法之一,经过几十年的发展,在仪器研发、正演模拟及反演解释等方面都已取得了突破(Li and Oldenburg, 1992, 1994; Kemna et al., 2012; Karaoulis et al., 2013; McGill et al., 2015; Liu et al., 2017).但强电磁干扰的存在依然限制了激电法在三维大规模探测中的应用.比如在宽频带激电观测中,由天然大地电流引起的低频趋势项漂移,会给激电相位数据带来极大干扰.在矿集区探测中,矿山钻机施工、工业发电、高压电网等人文干扰也会给实测激电序列带来各种突发性强干扰.此外,由天电干扰引起的尖峰脉冲噪声,由环境背景引起的高斯、类高斯随机噪声,以及激电探测中固有的电磁耦合等,都会造成激电原始时间序列的畸变,从而给激电参数的提取带来较大误差.长期以来,国内外专家学者主要从以下方面开展激电抗干扰技术研究:在硬件方面,增大发射电流提高数据采集信噪比,通过选频接收、相干检测等增强接收机抗干扰能力(陈儒军等, 2003;杨洋等, 2013;叶俊麟和罗有春, 2015;付国红等, 2016;肖都等, 2016;冉军林和刘俊岩, 2018).在信号处理方面,主要采用线性拟合、均值叠加、频域滤波、小波分析以及人工手动挑选数据等方法对时间序列进行处理(张春贺等, 1998; Gazoty et al., 2013; Deo and Cull, 2015;李耀恩和吴小平, 2015; Olsson et al., 2015, 2016; Bernard, 2016;徐信等, 2017).还有一种抗干扰思路是采用系统辨识的思想,重新定义和计算介质复电阻率响应(赵璧如等, 2006;罗先中等, 2014;罗延钟等, 2016;罗维斌等, 2017).但是由于电磁干扰的复杂性,上述方法还需在提高工作效率、减小有用信号损失、完善方法理论等方面进一步改进.强干扰背景下的数据处理依然是激电大规模三维探测中的一个重要挑战.

近十年来,激电法的一个重要进展是实现了分布式全波形数据采集,为激电信号和噪声干扰的全程监控与评价提供了新的舞台.国外具有全波形采集能力的仪器主要包括:澳大利亚MIMEX公司于1994年研制出的MIMDAS三维电法系统、美国NEWMONT矿业公司于2010年推出的NEWDAS三维分布式激电采集系统、加拿大QUANTEC公司推出的TITAN 24电法系统和ORION 3D全景三维系统、法国IRIS公司推出的Fullwaver大功率全波形激电系统等(Kingman et al., 2007; Eaton et al., 2010; Gharibi et al., 2012;张亚伟等, 2015).国内,中国科学院地质与地球物理研究所、中国地质科学院地球物理地球化学勘查研究所、中南大学、吉林大学、中国地质大学(北京)、中国地质大学(武汉)、成都理工大学、中石油东方地球物理公司等单位也都推出了分布式电磁法探测系统,可开展包括激电法在内的各种电法、电磁法探测工作,并进行全波形数据记录(何继善,2010张文秀等,2012林品荣等,2015底青云等,2016张锐锋等,2016).

电流、电压数据的全波形记录使得激电勘探的数据量大大增加,对于每一个测点,不再是单一地输出某一采样点的视参数计算结果,而是可以完整的输出激励电流和各测道电压的全波形数据,进而提取时域极化率、积分充电率及频域复电阻率幅值、相位、频散率、Cole-Cole模型参数等信息.在地球物理勘探中,大地电磁、地震勘探等已经形成了相对成熟的数据处理流程,并引入了许多新的时间序列干扰压制技术(Liu et al., 2016; Chen et al., 2017; Li et al., 2017, 2018, 2019; Wang et al., 2017; Tang et al., 2018; Yuan et al., 2018),而激电法由于全波形采集技术刚刚得到推广,目前国内外都尚没有形成完善的时间序列抗干扰数据处理方法.为进一步推动激电法在深部矿产勘查中的应用,有必要研究相应的信号识别及干扰消除技术,从而利用全波形采集的优势,提高数据质量.

1 全波形激电抗干扰数据处理方法

分布式全波形激电勘探需要在多个测点同时进行长周期数据采集,不可避免的会遭遇到大地电流、天电干扰、工业强干扰、高斯随机噪声、电磁耦合等电磁干扰影响.常用的信号处理方法可以参考法国IRIS全波形激电系统(Bernard, 2016),主要包括:采用人工手动识别方法剔除电压序列中的大尺度干扰数据段,采用线性或非线性拟合消除电压数据的趋势项漂移,采用均值叠加对观测序列进行多周期平均,然后利用电压数据和电流数据计算时频域激电参数等.为了提高数据处理效率和精度,本文提出基于统计分析的抗干扰方法处理方法,主要包括:采用相关分析自动识别和剔除强干扰数据段,采用经验模态分解方法消除低频趋势项干扰,采用稳健统计同时压制高斯随机干扰和离群值干扰,采用分段稳健平均和相对相位谱提取时/频域激电参数等.

1.1 相关分析用于消除突发性强噪声干扰

在分布式激电探测中,每个测点的数据都经历了较长的采集时间,数据长度可能达到十几个周期或几十个周期.在不同的时间段内,噪声干扰水平不尽相同,甚至在某一时间段内可能出现由工业干扰、机械施工等引起的突发性强噪声,因此从原始时间序列中识别并剔除强干扰数据段至关重要,由于分布式观测产生数据量很大,靠人工手动识别强干扰数据段效率很低,无法满足实际需求.本文提出通过计算各周期电流电压数据的相关度来评价激电数据质量,识别并剔除强干扰数据段.首先将同步采集的原始电流数据及电压数据按周期分段,每周期为一数据段,对于每一周期数据段,相关系数计算公式为(Gayen, 1951):

(1)

其中,N为一个周期包含的采样点数,为电压或电流一个周期内所有采样点数的平均值,当发射电流为对称分布的双极性矩形波、伪随机2n序列波或扩频m序列波时,不含噪声干扰的纯电流电压信号均值为0.

因为电流信号是一个不含噪声干扰的标准信号,而实测电压信号是纯电压信号叠加了噪声干扰的混合信号,所以推导可得,噪声干扰越强,二者相关系数C(U, I)越低.在以往的文献中,利用相关分析进行大地电磁等数据挑选时多是直接给定阈值,相关度低于该阈值时,则认为该数据段干扰严重,从而剔除.本文提出对于多周期激电数据,可以根据5σ准则,剔除强噪声干扰数据段,即根据所有周期的相关系数分布情况,计算所有数据段相关系数C的中值和中值离差,公式为

(2)

(3)

然后剔除相关系数小于的周期,保留高质量数据用于后续处理,从而提高激电参数的计算精度.

分别以某矿集区实测的含强噪声干扰和不含强噪声干扰的两组信号为例,分别计算各周期电流电压相关系数,然后利用5σ准则设置阈值,剔除强干扰数据段,处理结果如图 1所示.

图 1 相关分析用于剔除突发强噪声干扰示意图 左侧分图从上至下依次为电流波形、含突发强干扰的电压波形、电流电压分段相关度(黑色直线为阈值,相关系数低于该阈值的数据段将被剔除)、剔除低相关性数据段后的电压波形.右侧分图从上至下依次为电流波形、不含突发强干扰的电压波形、电流电压分段相关度、剔除低相关性数据段后的电压波形. Fig. 1 Correlation analysis for removing sudden strong noise On the left, from top to bottom: Current waveforms, voltage waveforms with sudden strong interference, segmented correlation of current and voltage (black line is the threshold, and data segment with the correlation coefficient lower than the threshold will be removed), and the voltage waveform after removing the low-correlation data segment. On the right, from top to bottom: Current waveforms, voltage waveforms without sudden strong interference, segmented correlation, and the voltage data after eliminating low-correlation data segment.

图 1左侧分图可见,对于含有突发强噪声干扰的数据,噪声干扰越强的数据段,相关系数越小,反之,相关系数越高.根据相关系数的整体分布,利用5σ准则可以定位干扰数据段位置,有两个周期的数据因干扰严重被剔除.高品质的数据被保留.由图 1右侧分图可见,当原始数据不含明显干扰时,相关系数整体较高且都在5σ之上,所有数据都被保留,处理前后信号无变化.说明该方法不会造成有用信号的损失.

1.2 经验模态分解用于分离低频趋势项干扰

因为分布式激电观测时间很长,由大地电流、工业游散电流、电极极差等引起的随时间缓慢变化的低频电场也是多周期全波形激电时域序列中的一项重要干扰,变现为叠加在整体激电序列上的趋势项漂移.由于激电信息主要表现在微弱的二次电位信号中,趋势项的存在会使二次电位信号产生形变,对后续激电参数的提取造成不利影响,尤其会对激电低频相位数据的计算带来极大偏差.

常规消除趋势项干扰的方法主要是对激电序列做整体拟合来近似趋势项,包括:线性拟合、多项式拟合、指数函数拟合及形式更加复杂的函数拟合.基于函数拟合的去趋势项方法都要假设趋势项漂移随时间的变化满足某一特定函数.实际应用中该假设不一定成立,一旦选用的拟合函数及阶次不合理,会给激电信号带来新的形变.除了函数拟合之外,还有基于频域低通滤波的去趋势项方法,但趋势项干扰的频率范围往往与激电信号频率范围相重叠,低通滤波处理虽然消除了趋势项,也造成了激电低频信息的损失.

为了压制低频趋势项漂移,并尽量减少处理方法本身带来的附加干扰,本文提出通过经验模态分解(EMD)方法进行多周期激电数据的去趋势项处理.EMD方法是一种基于信号本身特征的非线性分解方法,关于其基本原理,国内外已有众多专家进行了阐述(Huang et al., 1998Li et al., 2019),在以往的文献中多是将该方法用于天然源随机平稳信号处理,如大地电磁信号、天然地震信号等,本文对EMD方法进行改进并引入多周期全波形激电时序处理,步骤如下:

(1) 首先识别出一个周期中信号的极大值和极小值位置,对于激电信号,极值位置位于矩形波上升沿及下降沿位置.

(2) 对于每一个极值点,找到其他所有周期中与其位置对应的极值点,并连线形成激电时间序列的包络线,然后通过线性插值,获得一条采样点数与原始序列相等的完整包络线.对上一步识别出的所有极值点都进行多周期识别及连线插值处理,从而获得多条完整的上下包络线.

(3) 对所有的上下包络线进行稳健平均,从而得到包络线均值,从原始激电时域信号中减去此包络线均值,得到新的激电序列,重复步骤(1)~(3),直至包络线均值近似为幅值为0的直线时,此时得到的为去趋势项后的纯激电信号.原始激电序列减去纯激电序列即为分离得到的趋势项干扰.

图 2展示了某矿集区不同测点No.1~No.4经EMD方法分离得到的纯激电信号及趋势项干扰.

图 2 测点No.1~No.4经EMD方法分离得到的纯激电信号及趋势项干扰 (a)、(b)、(c)分别为测点No.1原始信号序列、纯激电信号序列及趋势项干扰;往下每一行依次为测点No.2、测点No.3及测点No.4处理结果. Fig. 2 The separations of pure IP signal and trend drift using EMD method for survey points No.1~No.4 (a), (b) and (c) are the original IP signal, pure IP signal and trend drift respectively for the No.1; the following subfigures from top to bottom are the separated results of No.2, No.3 and No.4.

图 2可见,测点No.1因为信号较强,趋势项干扰幅值很小,对原始数据的干扰不明显,处理前后信号形态也没有明显差别.但对于测点No.2、No.3和No.4, 因为本身信号很弱,趋势项干扰影响严重,利用EMD方法可以将信号和干扰有效分离.为了进一步验证EMD方法的效果,分别利用干扰弱的测点No.1和干扰强的测点No.2的数据计算复电阻率,图 3展示了测点No.1和测点No.2去趋势项前后得到的复电阻率幅值和相位的计算结果.

图 3 不含趋势项干扰的测点No.1和含趋势项干扰的测点No.2处理前后的视电阻率和视相位对比 (a)测点No.1去趋势项前后视电阻率;(b)测点No.1去趋势项前后视相位;(c)测点No.2去趋势项前后视电阻率;(d)测点No.2去趋势项前后视相位. Fig. 3 Comparison of apparent resistivity and phase before and after processing of No.1 (without trend interference) and No.2 (with trend term interference) (a) Apparent resistivity before and after detrending of No.1; (b) Apparent phase before and after detrending of No.1; (c) Apparent resistivity before and after detrending of No.2; (d) Apparent phase before and after detrending of No.2.

图 3ab可以看到,对于测点No.1,原始数据不含趋势项干扰,处理前后复电阻率幅值和相位基本没有变化,说明该处理方法不会带来附加干扰.由图 3cd可以看到,对于测点No.2,数据趋势项干扰严重,处理后,电阻率谱和相位谱都得到改善,而且频率越低,处理之后效果越明显.对于长周期全波形激电序列,当无明显仪器零漂及电极极差时,经EMD分离得到的低频趋势项干扰主要由大地电流引起,还可以用来分析测区浅部不均匀体的电性分布,篇幅所限,本文不再论述.

1.3 稳健统计用于多周期激电时序叠加

完成原始数据的强突发干扰剔除及趋势项干扰分离之后,就需要对得到的多周期激电序列进行时序叠加,从而压制随机噪声,提高信噪比.均值叠加是各类数据处理中最常用的时序叠加方法,其主要基于最小二乘原理,当噪声干扰呈高斯分布时,可以实现最优估计.但是,对于强干扰区的激电探测,噪声干扰不仅包括呈高斯分布的环境背景噪声,也包括由天电干扰、大功率机械开关等引起的尖峰脉冲噪声,从而使需要叠加的采样点中存在离群值,噪声分布不再满足高斯分布,此时均值叠加就会成为有偏估计,因为少量离群值的存在就会给时序叠加结果带来较大偏差.

为了压制时序叠加中离群值干扰的影响,本文提出将稳健统计方法用于激电周期叠加.对于多次观测的某一采样点{yi},其中:

(4)

ξ1ξ2ξ3、…、ξN为误差序列,N为总观测次数.采用稳健统计方法由多次观测值y1y2y3、…、yN估计真实值θ,就是根据最大似然估计准则求解方程(Huber, 1964),即:

(5)

其中θ称为“位置参数”,表征观测数据的真实值.σ称为“尺度参数”,表征观测值与真实值之间的残差分布,ψ(x)函数为影响函数,有多种类型,如:“Cauchy”、“Huber”、“Hampel”等(Holland and Welsch, 1977).公式通过迭代求解,从而实现对于偏离位置参数的观测点给予低权值,降低其对统计过程的影响.

为了对比均值叠加和稳健叠加的效果,以某测点实测多周期数据为例,进行分析.首先对原始数据进行去突发干扰和去趋势项处理.然后进行时序叠加.分别抽取每一个采样点多周期重复观测的数值,如图 4a所示,周期总数为10.可见虽然大部分采样点在10次重复观测中变化不大,近似为直线,但仍有部分采样点在重复观测中出现了离群值(outliers).选取含离群值干扰的某一采样点,分别进行均值叠加和稳健叠加,结果如图 4b所示,可见,虽然重复观测的数值中只有一个离群值,但仍然使均值叠加结果偏离了群体值,而稳健叠加则因为自动识别并降低了离群值的权重,叠加结果是群体值的真实反映.

图 4 均值叠加和稳健叠加对比示意图 (a)不同采样点在多周期重复观测中的数值;(b)对含离群值的某一采样点分别进行均值叠加和稳健叠加所得结果. Fig. 4 Comparison of mean stack and robust stack (a) Repeated observations for various sampling sites; (b) Mean stack and robust stack results of repeated observations for one sampling site containing outliers.

对激电波形的每一个采样点都进行均值叠加和稳健叠加,分别得到一个完整的周期序列如图 5所示.

图 5 多周期全波形激电数据周期叠加示意图 (a)原始多周期全波形激电序列;(b)均值叠加所得激电波形;(c)稳健叠加所得激电波形序列. Fig. 5 Time series stacking of multi-period full-waveform IP data (a) Original multi-period full-waveform IP data; (b) IP data from mean stack; (c) IP data from robust stack.

图 5a可见,虽然原始激电序列经过相关分析和去趋势项处理大部分干扰已经消除,但仍有尖峰脉冲噪声保留在数据波形中.由图 5b可见,采用均值叠加得到的激电周期序列,依然受到原始多周期序列中脉冲干扰的影响,很多位置存在毛刺.而由图 5c可见,采用稳健叠加之后,脉冲干扰的影响被压制,叠加后的激电波形更加准确.

1.4 分段稳健平均和相对相位谱用于时/频域激电参数提取

对于分布式全波形激电法,无论发射电流波形为矩形波、伪随机2n序列波还是扩频m序列波,都可以基于全波形时域序列,同时提取时域极化率和频域激电相位等参数.提取时域激电参数的优点是可以通过增大延时,来避开电磁耦合的影响,缺点是二次场电位信号微弱,易受噪声干扰影响.提取频域激电参数的优点是,计算基于总场,抗干扰能力强于时域激电,但是因为无法躲开上升沿和下降沿,易受电磁耦合干扰影响.本文针对时/频域激电参数提取的特点,提出了相应的抗干扰技术.

时域激电参数提取主要受到噪声干扰的影响.时域极化率计算通常采用不同延时的二次电位与总电位比值(李金铭,2005),公式为

(6)

其中UpUsUt分别表示一次电位、二次电位及总电位,通常是某一延时的电位采样或是某一延时附近小范围时窗内的电位积分平均.另一种计算时域极化率的方法是采用某一时间段内二次电位积分比上总电位(何继善,2006),即:

(7)

经过前述的相关分析、经验模态分解和稳健统计等处理,虽然大部分强噪声干扰已被压制,但是当测区干扰特别强或观测周期不足,导致噪声干扰及离群值占原始数据序列的比例接近或超过50%时,便仍可能有残留噪声干扰保留到周期叠加后的二次场电位曲线中,第一种计算极化率的方法由于只采用了某一延时的单点或少量点采样数据,所以易受残留噪声干扰的影响.第二种计算方法采用积分方式,提高了计算精度,但是当积分段内二次电位存在严重形变时,同样会造成计算结果的偏差.

为了压制二次电位中的噪声干扰,提高时域极化率的计算精度.本文基于激电二次场电位早期变化快,晚期变化平缓的特点提出通过分段稳健平均的方式提取时域极化率,即将整段二次电位曲线从开始供电到断电之前,按对数间隔依次分为多段,每段为一个时窗,然后对每个时窗内的二次电位采样点进行稳健叠加平均作为该时窗中点的采样值,以压制残留噪声的影响.同时用一个小时窗对断电前的数据进行平均作为充分极化后的总电位值,早期不同时窗内的叠加均值作为不同延时的一次电位值,二者相减得到不同延时的二次电位值,然后利用公式(6)得到不同延时内对应的极化率值.图 6a为某段实测激电二次场电位信号窗口叠加示意图,图 6b为利用上述结果计算的不同延时极化率.由图 6可见利用该方法可有效提取不同延时的极化率,并压制二次场电位曲线中的噪声干扰.

图 6 分段稳健平均用于从激电二次场电位数据中提取不同延时极化率示意图 (a)原始激电二次场电位数据及各时窗稳健平均结果(时窗按对数间隔分布且分别由不同颜色代表);(b)不同延时计算得到的时域极化率. Fig. 6 Polarizability extraction from IP second potential data using robust stacking of various delayed time windows (a) Original IP second potential data and robust stacking results (time windows are distributed on a logarithmic interval represented by different colors); (b) Polarizability of various delay time.

频域激电参数提取主要受到电磁耦合干扰的影响.频域复电阻率幅值和相位的计算公式为(罗延钟和张桂青,1988):

(8)

其中,U(f)、I(f)、ρ(f)分别表示电压、电流及复电阻率频谱.实测U(f)是岩矿石在外电流激发下,激发极化效应和电磁耦合效应共同的频谱响应.在时间序列中,电磁耦合表现为激励电流供电和断电时引起的电压信号的短时突变,时间一般不超过10 ms;激发极化效应则是一个缓慢变化的过程,在激励电流供电和断电之后,往往还要经过几十秒的时间,才能达到饱和(何继善,2006).在频率域,随着频率从低到高,电磁耦合效应逐渐增强,激发极化效应逐渐减弱.当频率低到一定程度时,激发极化效应趋于饱和,随频率缓慢平稳变化,而电磁耦合依然随频率近似线性快速变化(Hallof, 1974; Coggon, 1984; Pelton et al., 1978).电阻率参数主要由信号较强的一次场决定,整体时间序列叠加上瞬时突变的电磁耦合干扰后,所受影响也较小,但相位或极化率等激电参数由信号较弱的二次场决定,叠加上电磁耦合干扰后,受影响较大而产生畸变.

为了压制频域激电中电磁耦合的影响,激电信号应尽量采用低频(10-3~100Hz)段观测,并计算相对相位来代替相位.因为在低频段内,由激电效应引起的复电阻率相位随频率缓慢平稳变化,而由电磁耦合引起的相位随频率近似线性快速变化,计算相对相位谱是对相位数据做线性校正,可以消除电磁耦合的影响.某一个频点的相对相位由该频点及相邻频点的相位计算得到,公式为

(9)

其中,f1f2为由低到高相邻的两个频点,为了对相对相位去除电磁耦合的效果进行分析,分别选取耦合感应较强和较弱的两个测点计算相位和相对相位,结果如图 7所示.

图 7 耦合感应较强和较弱的两个测点相位和相对相位计算结果 (a)耦合感应较强的某测点供电电流;(b)耦合感应较弱的某测点供电电流;(c)耦合感应较强的该测点电压波形;(d)耦合感应较弱的该测点电压波形;(e)耦合感应较强的该测点相位和相对相位对比;(f)耦合感应较弱的该测点相位和相对相位对比. Fig. 7 hase and relative phase of two survey points with strong and weak EM coupling, respectively (a) Current of survey point with strong coupling; (b) Current of survey point with weak coupling; (c) Potential difference waveforms of survey point with strong coupling; (d) Potential difference waveforms of survey point with weak coupling; (e) Phase and relative phase of survey point with strong coupling; (f) Phase and relative phase of survey point with weak coupling.

图 7ace可见,该测点激电效应较弱而电磁耦合较强,电压信号在对应电流上升沿和下降沿处出现突变.复电阻率相位谱受耦合影响较大,随频率近似直线上升,在高频时甚至变为正值,而相对相位谱受耦合感应影响较小,整体仍为负值,且平缓变化,主要反映激电信息.由图 7bdf可见,另一测点电磁耦合效应弱而激电效应较强,电压信号在对应的供电时间段内缓慢变化,相位谱和相对相位谱差别不大,都可以反映激电信息.上述对比说明计算相对相位可以对耦合干扰严重的数据进行滤波,同时不会对不含耦合干扰的数据带来偏差.需要说明的是,本文的频域去耦方法仅适用于观测频率较低的激电数据,如果激电观测中采用的频率较高,电磁耦合干扰已完全压制了激电曲线,则目前并无较好方法进行直接去耦处理.有关学者已经开展了同时考虑激电效应和电磁耦合的正反演研究(王大勇等,2010范翠松等,2014曹金华等,2017薛园兵等,2017),可以采用其方法进行含耦合效应的激电数据反演.

2 激电三维实测数据抗干扰处理及应用分析

基于提出的相关分析、经验模态分解、稳健叠加及时频域激电参数提取等算法,编制了一套完整的全波形激电数据抗干扰自动处理程序,并应用于多个典型矿区实测的全波形激电数据处理.下面以在西藏扎西康矿集区三维分布式激电探测全波形数据为例,分析不同条件下抗干扰数据处理的效果.

2.1 分布式三维激电数据采集

扎西康铅锌多金属矿床是我国西南部以富含硫盐矿物为重要特征的大型铅锌锑银共生矿床之一,容矿岩石为含炭钙质板岩、钙质板岩、绢云母板岩、页岩和石英砂岩.矿体受近南北向和北东—南西向两组断裂控制(焦彦杰等, 2015).为了进一步圈定成矿有利区,探测浅层大地电性结构,开展了三维分布式激电探测.经过前期标本测量工作发现,在频率10-3~10-1Hz范围内进行观测可有效进行矿与非矿的区分.此次共布置100条左右的测线,每条线100个测点,点距20 m,线距40 m.观测装置包括中梯扫面装置及非常规四极测深装置.采用的仪器为中南大学自主研发的分布式多通道电法观测系统,可一次性同时布置上百道采集站并行采集,采集阵列基于GPS同步和Zigbee通讯.发射电流波形为五频组合波,采样频率64 Hz, 基频为1/256 Hz,周期长度为256 s.对于每一个测点,数据采集时间为十个周期以上,通过全波形记录及抗干扰处理,可提取出1/256(0.0039) Hz、2/256(0.0078) Hz、4/256(0.0156) Hz、8/256(0.0313) Hz、16/256(0.0625) Hz等五个频点的复电阻率幅值、相位、频散率及不同延时的极化率.受篇幅所限,下面仅通过频域复电阻率幅值和相位的计算结果来分析抗干扰数据处理效果,频散率和极化率结果与相位数据具有类似的形态.

2.2 不同频率/不同极距下抗干扰处理效果分析

由于此次激电观测采用的频率很低,而且测区处于正在开采的矿山周边,低频趋势项干扰及各类电磁噪声干扰对激电信号带来严重畸变,需要采用抗干扰技术来提高数据质量.首先以某测线为例,分别采用抗干扰方法及常规方法对每一个测点的原始数据进行处理,分析其效果.常规方法按照法国IRIS全波形激电系统的处理流程进行(Bernard, 2016),主要包括:采用手动识别方法剔除电压序列中的大尺度干扰数据段,采用线性拟合消除电压数据的趋势项漂移,采用均值叠加对观测序列进行多周期平均,然后利用电压数据和电流数据计算复电阻率和相位谱.抗干扰方法处理流程包括:采用相关分析自动剔除强干扰数据段,采用经验模态分解方法消除趋势项干扰,采用稳健统计进行多周期平均,采用相对相位代替复电阻率相位等.图 8展示了某测线不同频率复电阻率幅值、相位的剖面计算结果,该测线采用中梯装置,AB=5000 m,测线长为2000 m,点距20 m,所有测点一次性布设完成,同步观测.

图 8 某测线抗干扰处理前后不同频率视电阻率和视相位曲线(相位取负) 从上到下依次为0.0039 Hz、0.0078 Hz、0.0156 Hz、0.0313 Hz、0.0625 Hz. Fig. 8 Apparent resistivity and phase at various frequencies for a survey line (phase value is negative) From top to bottom, frequency is 0.0039 Hz, 0.0078 Hz, 0.0156 Hz, 0.0313 Hz, 0.0625 Hz, respectively.

图 8可见,电阻率数据计算结果受噪声干扰影响相对较小,采用线性拟合和均值叠加等常规方法和本文提出的抗干扰方法进行处理,结果相差不大.这是因为电阻率参数表征介质导电性,计算结果主要受总场影响,总场信号强,而且原始数据观测周期较多,通过多周期均值叠加等常规方法也能压制噪声干扰的影响.但是对于相位数据,则受噪声干扰影响较大,这是因为相位参数表征介质的激电性,电压和电流的相位差主要受二次场影响,因为本身变化范围小,易受噪声干扰影响,抗干扰处理后,相位数据质量得到极大提升.同时对比不同频率的处理结果可见,整体而言频率越低,激电数据受噪声干扰的影响越大,抗干扰处理效果越明显.

为了进一步分析抗干扰数据处理效果,分别对不同供电极距下的某测线数据进行处理.该测线采用多组极距的中梯装置,测线长2000 m,点距20 m,所有测点一次性完成布设,每次供电,所有测点同步采集.供电极距依次是3900 m、3600 m、3300 m、3000 m、2700 m.图 9显示了不同极距下均值叠加等常规方法和抗干扰方法所得0.0039 Hz对应的复电阻率幅值相位处理结果.

图 9 某测线抗干扰处理前后不同发射极距视电阻率和视相位曲线(相位取负) 从上到下电极极距AB依次为3900 m、3600 m、3300 m、3000 m、2700 m. Fig. 9 Apparent resistivity and phase of various electrode space AB for a survey line (phase value is negative) From top to bottom, electrode spacing AB is 3900 m, 3600 m, 3300 m, 3000 m, 2700 m respectively.

图 9可见依然是电阻率数据处理前后相差不大,而相位数据得到极大改善,剖面曲线更加光滑合理.整体而言,供电电极极距越大,噪声干扰越严重,抗干扰处理效果越好.而且对于大极距供电,中间部分的测点受噪声干扰严重,这是因为供电极距越大,中间部分测点信号越弱,所以需要进行抗干扰处理,以保证深部探测的数据质量.

2.3 三维激电抗干扰处理结果

经上述分析可知,通过对每一个测点进行抗干扰处理,可有效提高大极距、低频点时整体剖面的激电数据质量.将该处理方法应用于测区全体数据,得到了超过5000个测点的高品质大规模三维数据体.图 10显示了由部分中梯扫面数据计算得到的最低频视电阻率和视相位平面图.图 11显示了由部分测线测深数据计算得到的最低频视电阻率和视相位断面图,视深度由AB/4表示.

图 10 对某测区大规模中梯扫面数据进行抗干扰处理后计算得到的视电阻率和视相位平面图 (黑色实心点为测点位置). Fig. 10 Apparent resistivity and phase pseudo planes of various scanning survey lines after anti-interference processing (Black dot is measurement site)
图 11 对某测区不同测线测深数据进行抗干扰处理后计算得到的视电阻率和视相位断面图 Fig. 11 Apparent resistivity and phase pseudo sections of various sounding survey lines after anti-interference processing

图 10可见,对中梯扫面数据进行抗干扰处理,进而得到的视电阻率和视相位拟平面图可有效反映测区大地整体的导电性和激电性分布,可有效区分出低阻高极化区、高阻中极化区、高阻低极化区及低阻低极化区.由图 11可见,对不同测线的测深数据进行抗干扰处理,进而得到视电阻率和视相位拟断面图,可对垂向上低阻高极化体的埋深及规模进行初步刻画.在大规模实测数据中的应用验证了本文提出的抗干扰数据处理的效果.综合扫面及测深结果,再结合其他地质资料进行联合反演,可对测区电性结构做出进一步精细解释.

3 结论

基于近年来我国在大功率分布式电法勘探系统方面的研制突破,本文针对多周期全波形激电数据提出了一套完整的抗干扰数据处理方法.通过对西藏扎西康铅锌矿矿区激电探测数据进行处理分析,可以得到以下结论:

(1) 在激电探测中,激电参数相比于电阻率参数,更易受到噪声干扰的影响.线性拟合、均值叠加等常规方法虽然可以压制电阻率参数中的噪声干扰,但对激电参数中的噪声干扰压制能力不足.

(2) 本文提出的相关分析、经验模态分解、稳健统计、多参数提取等方法可以构成一个方法体系,对激电时间序列中的突发强噪声、离群值、低频趋势项和电磁耦合等干扰进行压制,提高激电参数的数据质量.

(3) 野外数据采集的频点越低、发射极距越大,观测数据受噪声干扰的影响越严重,抗干扰方法的处理效果越明显.

References
Bernard J. 2016. Data acquisition and processing of a distributed 3D induced polarisation imaging system. ASEG-PESA-AIG 25th Geophysical conference & exhibition, 1-57.
Cao J H, Li T L, Liu Y L, et al. 2017. Study on the influence of IP and EM effects of 3D resistivity modeling results. Progress in Geophysics (in Chinese), 32(2): 579-583. DOI:10.6038/pg20170217
Chen R J, Luo W B, He J S, et al. 2003. The data acquisition system in the high-precision multi-frequency electric method. Geophysical and Geochemical Exploration (in Chinese), 27(5): 375-378.
Chen Y K, Zhou Y T, Chen W, et al. 2017. Empirical low-rank approximation for seismic noise attenuation. IEEE Transactions on Geoscience and Remote Sensing, 55(8): 4696-4711. DOI:10.1109/TGRS.2017.2698342
Coggon J H. 1984. New three-point formulas for inductive coupling removal in induced polarization. Geophysics, 49(3): 307-309. DOI:10.1190/1.1441665
Deo R N, Cull J P. 2015. Denoising time-domain induced polarisation data using wavelet techniques. Exploration Geophysics, 47(2): 108-114.
Di Q Y, Lei D, Wang Z X, et al. 2016. An integrated test of the multi-channel transient electromagnetic system. Chinese Journal of Geophysics (in Chinese), 59(12): 4399-4407. DOI:10.6038/cjg20161201
Eaton P, Anderson B, Queen S, et al. 2010. NEWDAS—The Newmont distributed IP data acquisition system.// SEG Technical Program Expanded Abstracts. Denver, Colorado: SPE, 1768-1772.
Fan C S, Li T L, Wang D Y, et al. 2014. Research on 2.5D SIP inversion with topography. Acta Geologica Sinica (in Chinese), 88(4): 755-762.
Fu G H, Pan Z, Cheng H, et al. 2016. Decoupling model and simulation analysis from signal transmitters in frequency domain IP. Progress in Geophysics (in Chinese), 31(4): 1569-1574. DOI:10.6038/pg20160421
Gazoty A, Fiandaca G, Pedersen J, et al. 2013. Data repeatability and acquisition techniques for time-domain spectral induced polarization. Near Surface Geophysics, 11(4): 391-406. DOI:10.3997/1873-0604.2013013
Gayen A K. 1951. The frequency distribution of the product-moment correlation coefficient in random samples of any size drawn from non-normal universes. Biometrika, 38(1-2): 219-247. DOI:10.1093/biomet/38.1-2.219
Gharibi M, Killin K, McGill D, et al. 2012. Full 3D acquisition and modelling with the quantec 3D system-the hidden hill deposit case study. ASEG Extended Abstracts, 2012(1): 1-4.
Hallof P G. 1974. The IP phase measurement and inductive coupling. Geophysics, 39(5): 650-665. DOI:10.1190/1.1440455
He J S. 2006. Dual-Frequency IP Method (in Chinese). Beijing: Higher Education Press.
He J S. 2010. Wide-Field Electromagnetic Method and Pseudo-Random Signal Electrical Method (in Chinese). Beijing: Higher Education Press.
Holland P W, Welsch R E. 1977. Robust regression using iteratively reweighted least-squares. Communications in Statistics-Theory and Methods, 6(9): 813-827. DOI:10.1080/03610927708827533
Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193
Huber P J. 1964. Robust estimation of a location parameter. The Annals of Mathematical Statistics, 35(1): 73-101. DOI:10.1214/aoms/1177703732
Jiao Y J, Liang S X, Guo J, et al. 2015. Comparative research on the combinational test of geophysical methods in the Zhaxikang lead-zinc ore concentration area, Tibet. Geophysical and Geochemical Exploration (in Chinese), 39(2): 245-252.
Karaoulis M, Revil A, Tsourlos P, et al. 2013. IP4DI: A software for time-lapse 2D/3D DC-resistivity and induced polarization tomography. Computers & Geosciences, 54: 164-170.
Kemna A, Binley A, Cassiani G, et al. 2012. An overview of the spectral induced polarization method for near-surface applications. Near Surface Geophysics, 10(6): 453-468. DOI:10.3997/1873-0604.2012027
Kingman J E E, Donohue J G, Ritchie T J. 2007. Distributed acquisition in electrical geophysical systems.// Proceedings of Exploration 07: 5th Decennial International Conference on Mineral Exploration. Toronto, Canada, 425-432.
Li G, Xiao X, Tang J T, et al. 2017. Near-source noise suppression of AMT by compressive sensing and mathematical morphology filtering. Applied Geophysics, 14(4): 581-589. DOI:10.1007/s11770-017-0645-6
Li J, Zhang X, Gong J Z, et al. 2018. Signal-noise identification of magnetotelluric signals using fractal-entropy and clustering algorithm for targeted de-noising. Fractals, 26(2): 1840011. DOI:10.1142/S0218348X1840011X
Li J, Zhang X, Tang J T, et al. 2019. Audio magnetotelluric signal-noise identification and separation based on multifractal spectrum and matching pursuit. Fractals, 27(1): 1940007. DOI:10.1142/S0218348X19400073
Li J M. 2005. Geoelectric Field and Electrical Exploration (in Chinese). Beijing: Geological Press.
Li Y, Oldenburg D W. 1992. Approximate inverse mappings in DC resistivity problems. Geophysical Journal International, 109(2): 343-362. DOI:10.1111/j.1365-246X.1992.tb00101.x
Li Y E, Wu X P. 2015. Application of wavelet transform on the data processing of dual-frequency IP survey. Progress in Geophysics (in Chinese), 30(1): 460-465. DOI:10.6038/pg20150168
Li Y G, Oldenburg D W. 1994. Inversion of 3-D DC resistivity data using an approximate inverse mapping. Geophysical Journal International, 116(3): 527-537. DOI:10.1111/j.1365-246X.1994.tb03277.x
Lin P R, Zheng C J, Wu W L, et al. 2015. Techniques and systems for large-depth and multi-function electromagnetic survey. Geological Survey of China (in Chinese), 2(8): 60-66.
Liu W, Cao S Y, Chen Y K. 2016. Applications of variational mode decomposition in seismic time-frequency analysis. Geophysics, 81(5): V365-V378. DOI:10.1190/geo2015-0489.1
Liu W Q, Lin P R, Lü Q T, et al. 2017. Time domain and frequency domain induced polarization modeling for three-dimensional anisotropic medium. Journal of Environmental and Engineering Geophysics, 22(4): 435-439.
Luo W B, Li M, Chen R J. 2017. Induced polarization method based on circular cross-correlation method and field experiment. Progress in Geophysics (in Chinese), 32(3): 1393-1398. DOI:10.6038/pg20170359
Luo X Z, Li D W, Peng F P, et al. 2014. Implementation and applications of an coded electrical instrument with anti-interference ability. Progress in Geophysics (in Chinese), 29(2): 944-951. DOI:10.6038/pg20140263
Luo Y Z, Lu Z G, Sun G L, et al. 2016. Algorithms of second solution for instruments of electrical and electromagnetic prospecting using pseudo random signal. Progress in Geophysics (in Chinese), 31(6): 2604-2608. DOI:10.6038/pg20160634
Luo Y Z, Zhang G Q. 1988. Principle of Frequency Domain Induced Polarization Method (in Chinese). Beijing: Geological Press.
McGill D, Farquhar-Smith D. 2015. A comparison of 3D DCIP data acquisition methods. ASEG Extended Abstracts, 2015(1): 1-4.
Olsson P I, Dahlin T, Fiandaca G, et al. 2015. Measuring time-domain spectral induced polarization in the on-time: decreasing acquisition time and increasing signal-to-noise ratio. Journal of Applied Geophysics, 123: 316-321. DOI:10.1016/j.jappgeo.2015.08.009
Olsson P I, Fiandaca G, Larsen J J, et al. 2016. Doubling the spectrum of time-domain induced polarization by harmonic de-noising, drift correction, spike removal, tapered gating and data uncertainty estimation. Geophysical Journal International, 207(2): 774-784. DOI:10.1093/gji/ggw260
Pelton W H, Ward S H, Hallof P G, et al. 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP. Geophysics, 43(3): 588-609. DOI:10.1190/1.1440839
Ran J L, Liu J Y. 2018. The application and study of induced polarization group device. Geophysical and Geochemical Exploration (in Chinese), 42(6): 1259-1263.
Tang J T, Li G, Zhou C, et al. 2018. Power-line interference suppression of MT data based on frequency domain sparse decomposition. Journal of Central South University, 25(9): 2150-2163. DOI:10.1007/s11771-018-3904-7
Wang D Y, Li T L, Li J P, et al. 2010. Forward simulation of a 3D complex resistivity model. Progress in Geophysics (in Chinese), 25(1): 266-271. DOI:10.3969/j.issn.1004-2903.2010.01.035
Wang H, Campanyà J, Cheng J L, et al. 2017. Synthesis of natural electric and magnetic Time-series using Inter-station transfer functions and time-series from a Neighboring site (STIN): Applications for processing MT data. Journal of Geophysical Research: Solid Earth, 122(8): 5835-5851. DOI:10.1002/2017JB014190
Xiao D, Guo P, Lin P R, et al. 2016. Phase-IP experimental effects under the condition of strong interference and difficult grounding. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 38(5): 593-597.
Xu X, Li Z H, Yang Y. 2017. De-noising of IP data based on EMD-ICA. Application Research of Computers (in Chinese), 34(6): 1737-1739, 1744.
Xue Y B, Li J M, Li T L, et al. 2017. Research on forward modeling of three-dimensional complex electromagnetic method with topography. Progress in Geophysics (in Chinese), 32(2): 781-785. DOI:10.6038/pg20170244
Yang Y, Deng F H, Li D Q. 2013. The application of great depth induced polarization method based on pseudo-random signal to oil and gas exploration. Geophysical and Geochemical Exploration (in Chinese), 37(3): 438-442.
Ye J L, Luo Y C. 2015. The method for improving the grounding resistivity of large-power IP prospecting in Gobi region of Xinjiang. Geophysical and Geochemical Exploration (in Chinese), 39(6): 1156-1159.
Yuan S Y, Wang S X, Luo C M, et al. 2018. Inversion-based 3-D seismic denoising for exploring spatial edges and spatio-temporal signal redundancy. IEEE Geoscience and Remote Sensing Letters, 15(11): 1682-1686. DOI:10.1109/LGRS.2018.2854929
Zhang C H, Li J M, Jiang M, et al. 1998. The application of fractal geometry to the time-domain IP. Progress in Geophysics (in Chinese), 13(1): 79-83.
Zhang R F, Yan L J, Sun S M, et al. 2016. IP parameter extraction from TFEM data in the time domain. Oil Geophysical Prospecting (in Chinese), 51(6): 1227-1232.
Zhang W X, Zhou F D, Lin J, et al. 2012. Application of distributed electromagnetic system in deep groundwater prospecting. Journal of Jilin University (Earth Science Edition) (in Chinese), 42(4): 1207-1213.
Zhang Y W, Yan J Y, Zhang K, et al. 2015. Review of distributed 3D DC/IP method. Progress in Geophysics (in Chinese), 30(4): 1959-1970. DOI:10.6038/pg20150459
Zhao B R, Zhao J, Zhang H K, et al. 2006. The PS100 high precision earth-eletrictity instrument system (IP to IP) with controllable source-Application of CDMA technology to the measurement of earth-resistivity for the first time. Progress in Geophysics (in Chinese), 21(2): 675-682. DOI:10.3969/j.issn.1004-2903.2006.02.053
曹金华, 李桐林, 刘永亮, 等. 2017. 激电和电磁效应对三维复电阻率正演结果的影响研究. 地球物理学进展, 32(2): 579-583. DOI:10.6038/pg20170217
陈儒军, 罗维炳, 何继善, 等. 2003. 高精度多频电法数据采集系统. 物探与化探, 27(5): 375-378. DOI:10.3969/j.issn.1000-8918.2003.05.012
底青云, 雷达, 王中兴, 等. 2016. 多通道大功率电法勘探仪集成试验. 地球物理学报, 59(12): 4399-4407. DOI:10.6038/cjg20161201
范翠松, 李桐林, 王大勇, 等. 2014. 起伏地形下复电阻率法2.5维反演研究. 地质学报, 88(4): 755-762.
付国红, 潘志, 程辉, 等. 2016. 频率域激电发送端去耦建模与仿真分析. 地球物理学进展, 31(4): 1569-1574. DOI:10.6038/pg20160421
何继善. 2006. 双频激电法. 北京: 高等教育出版社.
何继善. 2010. 广域电磁法和伪随机信号电法. 北京: 高等教育出版社.
焦彦杰, 梁生贤, 郭靖, 等. 2015. 西藏扎西康铅锌矿集区的物探方法组合试验. 物探与化探, 39(2): 245-252.
李金铭. 2005. 地电场与电法勘探. 北京: 地质出版社.
李耀恩, 吴小平. 2015. 小波变换在双频激电数据处理中的应用. 地球物理学进展, 30(1): 460-465. DOI:10.6038/pg20150168
林品荣, 郑采君, 吴文鹂, 等. 2015. 大深度多功能电磁探测技术与系统集成. 中国地质调查, 2(8): 60-66.
罗维斌, 李梅, 陈儒军. 2017. 相关辨识谱激电法及野外实验. 地球物理学进展, 32(3): 1393-1398. DOI:10.6038/pg20170359
罗先中, 李达为, 彭芳苹, 等. 2014. 抗干扰编码电法仪的实现及应用. 地球物理学进展, 29(2): 944-951. DOI:10.6038/pg20140263
罗延钟, 陆占国, 孙国良, 等. 2016. 伪随机信号电法仪器的"第二方案"算法. 地球物理学进展, 31(6): 2604-2608. DOI:10.6038/pg20160634
罗延钟, 张桂青. 1988. 频率域激电法原理. 北京: 地质出版社.
冉军林, 刘俊岩. 2018. 组合激电测深装置的应用与研究. 物探与化探, 42(6): 1259-1263.
王大勇, 李桐林, 李建平, 等. 2010. 三维复电阻率模型电磁场正演模拟研究. 地球物理学进展, 25(1): 266-271. DOI:10.3969/j.issn.1004-2903.2010.01.035
肖都, 郭鹏, 林品荣, 等. 2016. 相位激电法在强干扰区的应用试验. 物探化探计算技术, 38(5): 593-597. DOI:10.3969/j.issn.1001-1749.2016.05.03
徐信, 李志华, 杨越. 2017. 基于EMD-ICA的激电数据降噪处理方法. 计算机应用研究, 34(6): 1737-1739, 1744. DOI:10.3969/j.issn.1001-3695.2017.06.030
薛园兵, 李京谋, 李桐林, 等. 2017. 带地形三维复电阻率正演研究. 地球物理学进展, 32(2): 781-785. DOI:10.6038/pg20170244
杨洋, 邓锋华, 李帝铨. 2013. 基于伪随机信号的大深度激发极化法在油气勘探中的应用. 物探与化探, 37(3): 438-442.
叶俊麟, 罗有春. 2015. 新疆戈壁地区大功率激电勘探接地电阻改善方法. 物探与化探, 39(6): 1156-1159.
张春贺, 李金铭, 姜枚, 等. 1998. 分形几何在激电时间谱中的应用. 地球物理学进展, 13(1): 79-83.
张锐锋, 严良俊, 孙社敏, 等. 2016. 时频电磁法时域激电参数提取与应用. 石油地球物理勘探, 51(6): 1227-1232.
张文秀, 周逢道, 林君, 等. 2012. 分布式电磁探测系统在深部地下水资源勘查中的应用. 吉林大学学报(地球科学版), 42(4): 1207-1213.
张亚伟, 严加永, 张昆, 等. 2015. 分布式三维电法研究进展. 地球物理学进展, 30(4): 1959-1970. DOI:10.6038/pg20150459
赵璧如, 赵健, 张洪魁, 等. 2006. PS100型IP到端可控源高精度大地电测仪系统—CDMA技术首次在地电阻率测量中的应用. 地球物理学进展, 21(2): 675-682. DOI:10.3969/j.issn.1004-2903.2006.02.053