2017, Vol. 60
2. 中国科学院矿产资源研究重点实验室, 中国科学院地质与地球物理研究所, 北京 100029;
3. 中国科学院大学, 北京 100049
2. Key Laboratory of Mineral Resources, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China
瞬变电磁方法(Transient Electromagnetic Method，TEM)是一种应用于矿产探测、地下水资源勘查等方面的常见地球物理方法(牛之琏，2007).根据噪声来源，TEM观测噪声可分为仪器内部噪声和外部噪声(McCracken et al., 1984；Olsen et al., 1988；McNeill and Labson, 1991；Barr et al., 2000).仪器内部噪声主要来源于磁场传感器，其中影响最大的是传感器前放电路噪声(刘凯等，2014；Liu et al., 2015).仪器外部噪声主要包括天然电磁场噪声、人文电磁场噪声以及感应线圈在地球静磁场中运动所产生的噪声.噪声水平的高低影响TEM观测的精度与有效探测深度(Spies, 1989)，在研发低噪声观测系统以降低仪器内部噪声的同时，开展噪声抑制方法研究以进一步降低TEM观测噪声的整体水平也尤为重要.
传统TEM采用双极性同步采样(Bipolar Synchronous Sampling)方法进行噪声抑制，对具有平稳随机特性的噪声(如远距离天电噪声、传感器电路噪声、工频噪声等)具有良好的抑制效果(蒋邦远，1998)，是目前国内外TEM仪器噪声抑制方面的标准技术(McCracken et al., 1984；嵇艳鞠等，2003).
在实际中，双极性同步采样方法的噪声抑制性能仍有不足之处：(1) 由于实际工频与理论值之间可能存在一定偏差，按照理论工频选择发射周期并不能很好实现对实际工频及其谐波的抑制；(2) 晚期信号带宽有限，现有方法晚期取样道带宽相对信号带宽余量较大，因此进入观测的高频噪声较多.为了改善这种情况，研究人员提出多种解决办法，如“滤除局部噪声观测法”、“自适应滤波法”等(Spies, 1988；李实等，2001；时圣利，2008；汤金云等，2008；李风明等，2014)，这些方法对噪声均有一定的抑制效果，但仍存在局限性(牛之琏，2007).为提升系统噪声抑制能力，国外学者开展了对常规双极性同步采样方法的优化研究，Macnae等提出在各周期进行加权叠加的思路(Macnae et al., 1984)，Nyboe和Sørensen提出对取样道加高斯窗的方法(Nyboe and Sørensen，2012).
本文在上述研究基础上，对双极性同步采样方法进行进一步优化，提升其噪声抑制性能：通过对取样道加高斯窗函数的方式优化取样道内各采样点的累加权重，提升对高频噪声的抑制能力；通过对包含所有叠加周期的完整观测序列加窗函数的形式优化各周期的叠加权重，提升对实际工频噪声的抑制能力.其中，对整观测序列所加的窗函数为高斯窗函数与50 Hz (及其倍频)梳状滤波器构成的复合窗函数.最后通过对比优化前后不同采样函数谱特性，及野外实测数据分析，结果证明本文提出的方法效果较好.
2 双极性同步采样方法的噪声抑制原理以数字TEM系统为例，其双极性同步采样过程如图 1所示：首先在配置文件中定义各取样道的中心时刻(c1, c2, …, ck)及时间宽度(Δ1, Δ2, …, Δk)；CPU按照取样道配置文件计算各取样道所含采样点的数据指针；信号经ADC采样后进入缓存区，CPU根据数据指针对取样道所包含数据进行处理，实现符号翻转、采样点累加平均等运算，得到单周期各取样道的输出值.对多个周期相同取样道的输出值再进行叠加平均，得到各取样道的最终输出值.
|
图 1 传统取样道采集方式 Fig. 1 Traditional gate function of sampling |
上述TEM同步采样过程可表示为：
|
(1) |
其中序列x表示多周期TEM观测信号，数值y表示经处理后的某单一取样道最终输出值.函数S描述至时刻t输出y值的计算方式，称作单个取样道的仪器采样函数(Instrument Sampling Function, ISF)，其定义为(Macnae et al., 1984)：
|
(2) |
其中，R (t)为参考函数，W (t)为叠加权函数.R (t)在完整TEM观测时长上定义各周期同一取样道内所含各采样点的累加权重.W (t)在完整TEM观测时长上定义各周期的叠加权重.
Macnae等(1984)证明，对于TEM原始信号中包含的平稳随机噪声，式(1) 中取样道输出值y (t)的方差σ2可以表示为：
|
(3) |
其中K (f)为噪声谱灵敏度函数(Noise Spectral Sensitivity Function)，n (f)为平稳随机噪声的单边谱密度.噪声谱灵敏度函数K (f)与S (t)的傅里叶变换S (f)之间的关系为：
|
(4) |
其中S*(f)是S (f)的共轭.
上述结果说明：可以通过优化ISF来优化噪声谱灵敏度函数，从而提升系统对噪声的抑制能力.
3 采样函数的精细优化方法为了实现对ISF优化，首先需要研究ISF的建立过程及其频谱特性.设周期为T的双极性冲激对为：
|
(5) |
周期为T的Shah函数为：
|
(6) |
(5) 式和(6) 式进行卷积计算，可以得到一个周期为T的双极性冲激序列，公式为：
|
(7) |
其傅里叶变换为：
|
(8) |
其中f0=1/T.
对(8) 式进行归一化后将之变换回时域得到得：
|
(9) |
接下来，在单周期内设计一个中心时刻为c、持续时间为Δ、幅度为1/Δ的取样道，公式为：
|
(10) |
其中Π表示矩形脉冲.将函数G (t)称作取样道函数，其傅里叶变换为：
|
(11) |
利用式(9) 和式(10) 卷积可以生成参考函数R (t)为：
|
(12) |
|
(13) |
FR (f)是R (t)的傅里叶变换.
与G (t)的定义方式类似，假定叠加周期数为N，可以定义叠加权函数W (t)为:
|
(14) |
|
(15) |
其中Q=NT表示总观测时长，w为总观测时长的时间中点，FW (f)为W (t)的傅里叶变换.
将式(12) 和式(14) 代入式(2) 即可得到传统双极性同步采样方法的ISF，记作SBB，其中S表示ISF，下标第一个字母代表取样道函数所使用的函数类型，下标第二个字母代表叠加权函数所使用的函数类型.SBB表示取样道函数与叠加权函数均使用矩形(Box)脉冲函数.
对传统ISF进行优化，实际上就是对G (t)与W (t)进行优化.Nyboe指出基于高斯窗函数良好的低通性能(Nyboe，2011)，可选择高斯窗函数作为ISF优化的基础.经理论分析与野外测试，本文认为对G (t)与W (t)的优化方式应有所区别：优化G (t)的主要功能是提升ISF的高频噪声抑制性能，降低FG (f)旁瓣的影响；优化W (t)的主要功能是提升ISF对工频噪声的抑制性能，要求对工频噪声具有更强的抑制能力，尤其对于实际工频相对其理论值有较大偏离的情况.基于上述考虑，构建优化取样道函数与叠加函数为：
|
(16) |
|
(17) |
GG(t)和WN(t)为优化后的参考函数和叠加权函数.其中，(17) 式表示优化后的叠加权函数是由函数WG(t)和WF(t)共同作用形成的复合窗函数，其中WG(t)定义为：
|
(18) |
表示对原矩形脉冲函数加一个高斯窗.WF(t)为使用全相位法(王兆华和黄翔东，2009)构建的50 Hz及其奇次谐波的梳状滤波函数.
4 优化效果比较为分别研究优化G (t)和W (t)对ISF的性能提升，拟使用式(16)、式(17) 以及式(18) 分别对应替代式(12) 和式(14) 形成三组结果，记作：SGB、SGG和SGN.SGB表示相对SBB仅使用式(16) 对取样道函数进行了优化；SGG表示使用式(16) 与式(18) 分别对取样道函数和叠加权函数进行优化；SGN表示使用式(16) 与式(17) 分别对取样道函数和叠加权函数进行优化.
表 1为计算SGB、SGG和SGN的具体参数，其描述的取样道大致等于Geonics公司ProTEM接收机25 Hz基频时的第26号取样道.
|
表 1 ISF生成算例相关的参数 Table 1 Parameters of a case of ISF generation |
首先对比SBB和SGB.图 2为SBB与SGB的时域形态比较.如图 2a所示，SBB的取样道函数为矩形脉冲函数，SGB的取样道函数为高斯函数.如图 2b反应整个观测时间内所有26号取样道的叠加权重，可见，SBB与SGB所有周期的叠加权重都相同.
|
图 2 SBB与SGB的时域形态比较 (a) SBB与SGB取样道函数对比；(b) SBB与SGB叠加权函数对比. Fig. 2 Comparison of time-domain forms between SBB and SGB (a) Sampling gate function; (b) Stacked weighting function. |
进一步对SBB与SGB进行谱域比较，如图 3a所示，由于使用高斯函数对取样道函数进行优化，从而使得SGB对于高频部分具有了更强的抑制能力.图 3b为图 3a低频局部放大，使用单对数坐标系.由该图可见，尽管使用高斯函数对取样道函数进行了优化，但对50 Hz的抑制能力两种ISF效果相同.由此表明，对取样道函数优化仅可以提升系统对高频噪声的抑制能力.
|
图 3 SBB与SGB的单边幅度谱比较 (a) SBB与SGB单边幅度谱对比；(b) SBB与SGB单边幅度谱低频细节对比. Fig. 3 Comparison of single-sided amplitude spectrum between SBB and SGB (a) Overall comparison; (b) Details at low frequency. |
下一步对比SGB和SGG.图 4为SGB与SGG的时域形态比较，如图 4a所示，SGB与SGG的取样道函数均为高斯函数.如图 4b所示，在整个观测时间内，SGB所有周期的叠加权重都相同，而SGG所有周期的叠加权重并不相同，其包络线表现为高斯函数形态.
|
图 4 SGB与SGG的时域形态比较 (a) SGB与SGG取样道函数对比；(b) SGB与SGG叠加权函数对比. Fig. 4 Comparison of time-domain forms between SGB and SGG (a) Sampling gate function; (b) Stacked weighting function. |
进一步对SGB与SGG进行谱域比较，如图 5a所示，由于使用了高斯函数进行了取样道函数优化，SGB与SGG对于高频部分均具有了较强的抑制能力.图 5b为图 5a的低频局部放大，使用单对数坐标系.由该图可见，由于SGG对各周期的叠加权重进行了优化，使其对50Hz附近的噪声抑制能力强于SGB.
|
图 5 SGB与SGG的单边幅度谱比较 (a) SGB与SGG单边幅度谱对比；(b) SGB与SGG单边幅度谱低频细节对比. Fig. 5 Comparison of single-sided amplitude spectrum between SGB and SGG (a) Overall Comparison; (b) Details at low frequency. |
最后对比SGG和SGN.如图 6a所示，由于SGN同样对取样道使用了高斯函数优化，因此其具有与SGG相同的高频噪声抑制能力.图 6b为图 6a低频局部放大，使用单对数坐标系.由该图可见，SGN在50 Hz处具有一个宽度约为1.2 Hz的陷波窗口，从而使其具有比SGG更强的工频抑制能力.
|
图 6 SGG与SGN的单边幅度谱比较 (a) SGG与SGN单边幅度谱对比；(b) SGG与SGN单边幅度谱低频细节对比. Fig. 6 Comparison of single-sided amplitude spectrum between SGG and SGN (a) Overall Comparison; (b) Details at low frequency. |
为验证本文提出的方法，中国科学院电子学研究所于2015年7月在广西自治区柳州市某矿区开展了方法验证试验.使用自主研发的CASTEM系统，布设100 m×100 m发射回线，平均发射电流14.3 A，关断时间21 μs，发射基频25 Hz，叠加周期数为75次.使用CASTEM-S2型感应式磁场传感器，线性带宽40 kHz，有效面积1000 m2.CASTEM系统接收机采样率为500 kHz，具有全时段采集、原始数据存储功能.
由于临近村落并有钻机作业，工作区域的外部电磁环境相对恶劣，如图 7所示，50 Hz磁场达到
|
图 7 工作区等效输入噪声 Fig. 7 Equivalent input noise in the work area |
以对应ProTEM第26号取样道的处理结果为例，图 8为分别使用SBB方式(图 8a)、SGG方式(图 8b)以及SGN方式(图 8c)得到的单位面积输出噪声幅度谱，此三者的对比实际为传统双极性采样方式、Nyboe提出的优化方式以及本文所提出的优化方式间的对比.
|
图 8 单位面积输出噪声幅度谱 (a)SBB方式的单位面积输出噪声幅度谱；(b)SGG方式的单位面积输出噪声幅度谱；(c)SGN方式的单位面积输出噪声幅度谱. Fig. 8 Comparison of output noise spectra normalized by area of three ISFs (a) Output noise spectrum normalized by area of SBB; (b) Output noise spectrum normalized by area of SGG; (c) Output noise spectrum normalized by area of SGN. |
图 8各图中黑线为外部噪声单位面积幅度谱，灰线为经不同方式抑制后的噪声单位面积幅度谱.通过比较可见，由于SGG及SGN方式所使用的取样道优化函数相同，因此两者对高频噪声的抑制效果相当，均显著强于SBB方式.图 9为图 8b和c中SGG与SGN方式在低频部分的细节展开，通过比较可见，SGN方式对50 Hz具有比SGG方式更强的抑制能力.
|
图 9 SGG与SGN方式工频噪声抑制效果比较 Fig. 9 Comparison of industrial noise reduction effect between SGG and SGN |
利用式(3) 对上述三种方式相同取样道的观测均方根误差进行计算，SBB方式为59 nV·m-2，SGG方式为44 nV·m-2，SGN方式为32 nV·m-2.以此为标准，SGN方式的噪声抑制能力相对传统的SBB方式提升了约46%，相对Nyboe优化方式SGG提升约27%.可见，通过本文提出方法对ISF进行优化，能够显著增强系统的噪声抑制能力.
6 结论为进一步提升瞬变电磁仪器的噪声抑制能力，本文在传统双极性同步采样方法的基础上提出优化方案：通过对取样道加高斯窗函数的方式，赋于取样道内各采样点合适的累加权重，提升系统对高频噪声的抑制能力；使用高斯窗函数与梳状滤波器构造复合窗函数，通过对包含所有叠加周期的完整观测序列加复合窗函数的形式优化各周期的叠加权重，提升系统对实际工频噪声的抑制能力.使用传统双极性采样方式SBB、一般优化方式SGG以及强化工频抑制优化方式SGN对野外实测数据处理，通过比较可见：对相同取样道使用SGN方式较SBB方式可提升系统噪声抑制能力幅度达46%，从而证明本文所提方法较传统方法具有更强的噪声抑制能力，为进一步推广打下基础.
致谢向为本研究提供可靠仪器及参与了野外测试工作的中国科学院电子所地球物理仪器研发团队成员刘丽华、李巨涛、朱万华、黄玲、吴凯、齐林、刘富波、闫彬以及毕业生刘凯博士表示衷心感谢！
| Barr R, Jones D L, Rodger C J. 2000. ELF and VLF radio waves. Journal of Atmospheric and Solar-Terrestrial Physics, 62(17-18): 1689-1718. DOI:10.1016/S1364-6826(00)00121-8 | |
| Ji Y J, Lin J, Cheng F D, Yu S B. 2003. Research on the signal processing method for the transient electromagnetic method. Geophysical & Geochemical exploration, 27(2): 142-145. | |
| Jiang B Y. 1998. Practical near-zone transient electromagnetic exploration based on magnetic source. Beijing: geological publishing house. | |
| Li F M, Liu H F, Zhang X J, et al. 2014. The influence of industrial transmission lines on transient electromagnetic's second magnetic field and the elimination of interference. Progress in Geophysics, 29(3): 1399-1405. DOI:10.6038/pg20140355 | |
| Li S, L i, C S, Zhang Y P, et al. 2001. Research on the key technique for transient electromagnetic system using on project exploration. Coal geology & exploration, 29(1): 55-58. | |
| Liu K, Mi X L, Zhu W H, et al. 2014. A design of high sensitivity induction magnetometer for TEM. Chinese Journal of Geophysics, 57(10): 3485-3492. DOI:10.6038/cjg20141034 | |
| Liu K, Zhu W H, Yan B, et al. 2015. Ultralow noise preamplifier and optimization method for induction magnetometers. IEEE Sensors Journal, 15(6): 3293-3300. DOI:10.1109/JSEN.2015.2389336 | |
| Macnae J C, Lamontagne Y, West G F. 1984. Noise processing techniques for time-domain EM systems. Geophysics, 49(7): 934-948. DOI:10.1190/1.1441739 | |
| McCracken K G, Oristaglio M L, Hohmann G W. 1986. Minimization of noise in electromagnetic exploration systems. Geophysics, 51(3): 819-832. DOI:10.1190/1.1442134 | |
| McCracken K G, Pik J P, Harris R W. 1984. Noise in EM exploration systems. Exploration Geophysics, 15(3): 169-174. DOI:10.1071/EG984169 | |
| McNeill J D, Labson V F. 1991. Geological mapping using VLF radio fields.//Nabighian M N ed. Electromagnetic Methods in Applied Geophysics-Volume 2 Application Part B, Investigations in Geophysics. Tulsa, Oklahoma:SEG, 521-640. | |
| Niu Z L. 2007. Theory of Time domain electromagnetic method. Changsha: Central south university Press. | |
| Nyboe N S. 2011. Improving instrumentation and measuring strategies for the transient electromagnetic method applied in hydrogeophysical investigations[Ph. D thesis]. Aarhus:Aarhus University. | |
| Nyboe N S, Sørensen K. 2012. Noise reduction in TEM:Presenting a bandwidth-and sensitivity-optimized parallel recording setup and methods for adaptive synchronous detection. Geophysics, 77(3): E203-E212. DOI:10.1190/geo2011-0247.1 | |
| Olsen R G, Deno D, Baishiki R S, et al. 1988. Magnetic fields from electric power lines:Theory and comparison to measurements. IEEE Transactions on Power Delivery, 3(4): 2127-2136. DOI:10.1109/61.194025 | |
| Shi S L. 2008. Study of several de-noising methods for transient electromagnetic signal[MSc thesis] (in Chinese). Changchun:Jilin University. | |
| Spies B R. 1988. Local noise prediction filtering for central induction transient electromagnetic sounding. Geophysics, 53(8): 1068-1079. DOI:10.1190/1.1442543 | |
| Spies B R. 1989. Depth of investigation in electromagnetic sounding methods. Geophysics, 54(7): 872-888. DOI:10.1190/1.1442716 | |
| Tang J Y, Yu J C, Wang Y Z, et al. 2008. Human noise influence on transient electromagnetic surveying and processing technology. Coal Mining Technology, 13(2): 32-34. | |
| Wang Z H, Huang X D. 2009. All-phase spectrum analysis and filtering technique for digital signal. Beijing: Publishing house of electronics industry. | |
| 嵇艳鞠, 林君, 程德福, 等. 2003. 瞬变电磁法中数据取样处理方法的研究. 物探与化探, 27(2): 142–145. | |
| 蒋邦远. 1998. 中华人民共和国地质矿产部地质专报-七, 第14号-普查勘探技术与方法实用近区磁源瞬变电磁法勘探. 北京: 地质出版社. | |
| 李实, 李创社, 张彦鹏, 等. 2001. 工程勘探瞬变电磁仪关键技术研究. 煤田地质与勘探, 29(1): 55–58. | |
| 李风明, 刘鸿福, 张新军, 等. 2014. 工业输电线对瞬变电磁二次场的影响及干扰的消除. 地球物理学进展, 29(3): 1399–1405. DOI:10.6038/pg20140355 | |
| 刘凯, 米晓利, 朱万华, 等. 2014. 一种用于TEM高灵敏度感应式磁场传感器设计. 地球物理学报, 57(10): 3485–3492. DOI:10.6038/cjg20141034 | |
| 牛之琏. 2007. 时间域电磁法原理. 长沙: 中南大学出版社. | |
| 时圣利. 2008. 瞬变电磁信号的几种去噪方法研究[硕士论文]. 长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-2008060743.htm | |
| 汤金云, 于景邨, 王扬州, 等. 2008. 矿井瞬变电磁人文噪声干扰及处理技术研究. 煤矿开采, 13(2): 32–34. | |
| 王兆华, 黄翔东. 2009. 数字信号全相位谱分析与滤波技术. 北京: 电子工业出版社. | |