2018, Vol. 61
2. 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083
2. Key Laboratory of Metallogenic Prediction of Non-Ferrous Metals and Geological Environment Monitor, Ministry of Education, Central South University, Changsha 410083, China
可控源电磁法(Controlled Source Electromagnetic Method，简称CSEM)是在大地电磁法的基础上发展起来的一种人工可控源电磁测深法，弥补了天然场源大地电磁法的不足(李金铭, 2005；何继善, 2007).广义上，CSEM包括所有采用人工源的电磁勘探方法，其最大的特点为采用人工源控制发射信号，有效信号是周期的.通过多周期采集能够压制噪声，提高野外采集信号的信噪比，如CSAMT、广域电磁法(WFEM)及时频电磁法(TFEM)等方法.但随着工业的发展及人类活动范围的不断扩大，各种电磁干扰越来越强，噪声问题变得日益严重.在很多情况下，尤其在矿集区电磁噪声的影响复杂且没有特定规律(汤井田和李晋等，2012；汤井田和徐志敏等，2012；汤井田和任政勇等2015；汤井田和刘子杰等, 2015, 李晋和汤井田等，2017)，易导致信号信噪比降低，有效频点减少，严重影响CSEM的纵向分辨率和探测效果，急需发展新的方法提高CSEM有效信号挖掘能力.
CSEM信号几乎都是基于伪随机编码(Pseudo Random Binary Signal)发送的，周期方波信号是伪随机信号的最简单形式，众多商用频率域电磁仪采用方波信号，一般采用扫频方式进行勘探，即通过不断改变主频，获得不同频率的电磁响应，野外施工效率并不高(汤井田，2005).为了提高勘探效率和信号抗干扰能力，伪随机信号被引入电磁勘探领域，其中具有代表性的2n序列伪随机信号(何继善, 2010a, 2010b)能够同时采集多个频率信息，而且对数坐标下其频率域能量分布均匀，更符合勘探需要.伪随机信号电法的核心是发送机将不同频率的电流波形合成为伪随机电流后向地下供电，接收机同时接收这些频率经大地后的响应并将其分离，因此，一次供电，便可完成多频率的测量，具有快速、高效、电源利用率高等优点.
理论上，方波和伪随机信号均含有大量谐波成分，如果能充分利用谐波成分，将极大的提高CSEM的频率密度，进而提高勘探的分辨率.以Zonge为代表的学者将方波的奇次谐波提取出来加以利用，称为“奇次谐波法”，并在GDP系列仪器中保留该方法(Zonge, 1981).但在野外工作中，有效信号总会受到各种电磁噪声的干扰，不同频率均受到噪声的影响而导致信噪比不同程度的降低.在这样情况下，有学者认为谐波能量不大，更易受到噪声干扰等原因，倾向于只应用主频信号，或根据经验选取部分低阶谐波(席振铢, 1996).复杂的伪随机信号同样要面对这样的问题，而且其谐波成分更加复杂.
如果能够快速提取伪随机信号中的所有有效信号(主频和其谐波)，在野外施工条件理想情况下，只需发射一组伪随机信号(2n序列13频波伪随机信号)就可以获取整个频段的高频率密度的有效信号，能够极大的提高电磁法的纵向分辨率和能源利用率.从公开的文献来看，目前还没有一种标准判断谐波什么时候可以利用，什么时候不可以利用.在野外勘探数据中，噪声在频率域的影响是不均匀的，可能是对谐波影响大，对主频影响小；也可能是对主频影响大，而对谐波频率影响小；如果是后者，高阶谐波虽然在幅度小，但相比主频和低阶谐波拥有更多的叠加次数，稳定性更好，此时可以将其提取并加以利用，而不是以往简单的认为幅度大，就更可靠，而幅度小就不可靠.
本文提出一种在频率域基于小波变换和Hilbert解析包络的噪声评价方法，首先对原始时间域数据进行混合基离散傅里叶变换得到准确频谱曲线，在CSEM频率位置对频谱进行预处理，同时保持其他频率频谱不变，基于小波分解(Mallat, 1989；Daubechies, 1992；Shensa 1992)将预处理后频谱分成低频和高频部分，对高频部分进行基于希尔伯特变换(Vakman，1972；Boashash, 1992；Picinbono, 1997；Feldman, 2011)的解析包络，获得上包络线，并与低频部分重构后获得频谱整体的上包络线(上界).在CSEM频率，计算包络线数值与原始振幅的比值(噪声评价系数)，估计噪声的最大影响程度，根据阈值筛选出高信噪比的主频和谐波信号.在评价过程中，所有主频和谐波均为候选的有效信号，如果包络线数值与CSEM频率频谱之比小于某个阈值(如5%)，则将该频率信息提取并加以利用，反之则丢弃处理.本方法可应用于所有频率域电磁勘探信号的有效频率筛选.
1 基本理论人工可控源电磁勘探方法的发射信号几乎都是周期的，其有效信号频谱在频率域是离散的，只在有限个频率位置存在能量.值得注意的是：由于很多电磁勘探仪器的信号采样频率并不是2的n次方，并不能直接采用基于2n的快速傅里叶变换(Welch, 1967；Marple, 1999)计算频谱，而应该采用更具广泛性的混合基快速傅里叶变换(Mix-Radix FFT)(Singleton，1969)，只有在此情况下，可控源信号的频谱才会是离散的，否则将发生能量泄漏造成频谱不离散.
基于有效信号的周期性特点，本文将信号分成周期信号和非周期信号.周期信号包括CSEM信号和周期噪声，周期噪声在频率域很容易和CSEM信号分开，为了表述方便本文中不再涉及周期噪声的影响，后面文中所指噪声均为非周期噪声.周期信号的频率将被统称为CSEM频率(包含信号主频和其谐波)；非周期信号则全来自噪声.由于工频噪声并不只存在于50 Hz，往往是50 Hz附近的一个频带，在本文中属于非周期噪声范畴.本文所处理信号为完整时间域信号，而不仅仅是有限个CSEM频率的频率域系数.
CSEM有效信号的频谱为一系列的离散脉冲，噪声的频谱则在频率域是连续的，能量存在于整个频率域.图 1为一组仿真信号的频谱，主频为1 Hz，其幅值为60 mV, 并包含大量谐波成分(3~199 Hz)，谐波能量以
|
图 1 仿真信号的频谱 Fig. 1 Spectrum of simulated signal |
尖脉冲位置即为CSEM频率位置, 当不存在噪声干扰时，信号将只在脉冲位置存在能量，在其他频率位置幅值为零.但实测数据总是存在着各种各样的电磁干扰，CSEM频率系数受到噪声不同程度的影响.此时，在CSEM频率位置，实测数据的频率域系数等于非周期部分(噪声)和周期部分(CSEM信号)系数在复平面内的向量和.
如果能够估计CSEM频率受噪声影响的最大值，CSEM频率对应的真实频率域系数将分布于复平面中的黑色虚线框内(见图 2).定义噪声评价系数为噪声的最大影响值(虚线圈半径)与CSEM频率幅值的比值，可用噪声评价系数评价噪声的影响程度，并筛选出符合质量或规范要求的高信噪比信息.
|
图 2 在CSEM频率噪声对频率域系数的影响示意图(对幅度和相位的影响) Fig. 2 The possible biggest influence of noise with certain amplitude on CSEM frequency coefficient (influence on amplitude and phase) |
通过该算法获得的噪声评价数同时能够反映相位受噪声影响程度，对噪声评价数进行反正弦运算即可获得相位角度误差的范围，相位角误差计算见式(1):
|
(1) |
其中k为噪声评价数，即噪声最大影响与实测数据幅值的比例，φ为相位误差最大角度.
噪声评价数为3%时，相位角误差为±1.72°；噪声评价数为5%时，相位角误差为±2.86°；噪声评价数为10%时，相位角误差为±5.74°.
2 算法 2.1 噪声的频谱特征非周期噪声在频率域的频率谱是振荡的，噪声的不同出现时间对应着不同的时移因子，经傅里叶变换后在频率域表现为振荡.多个噪声在频率域叠加时，表现为锯齿状振荡.
信号的频率谱形态特征与信号在时间域的分布和形态紧密相关，即使相同噪声出现在不同时间，在频谱上的特征也将不同.为说明噪声频谱的复杂振荡特性，本文仿真了几种简单噪声及组合后的信号，通过傅里叶变换计算得出它们的频谱曲线来说明噪声频谱的复杂性.
首先，本文仿真一个高斯函数曲线，并作为一种简单的非周期噪声，对其进行傅里叶变换计算频谱，傅里叶变换及其频谱计算公式，见式(2)和(3).
|
(2) |
|
(3) |
其中，f(t)为时间域信号，t为时间，e－π(t－a)2为高斯函数曲线，a为时移因子(即噪声在时间域出现时间)，
图 3为该噪声仿真信号在时间域和频率域的曲线特征，其频谱仍然是一个高斯函数(图中只保留了正频率部分).当存在单个简单噪声时，其频谱特征依然能够保持良好的形态，这是由于频谱中振荡因子e－2πias的绝对值为1.
|
图 3 一个高斯函数信号的时间域(a)和频率域(b)特征 Fig. 3 One Gaussian function in time domain (a) and its Fourier transform (b) |
本文在此基础上，仿真了两个高斯函数曲线组成复合噪声(见图 4)，其频率域系数为两个噪声频率域系数的叠加(复数域)，其傅里叶变换和频谱计算公式如下，见式(4)和(5).
|
图 4 两个高斯函数的时间域(a)和频率域(b)特征 Fig. 4 Two Gaussian functions in time domain (a) and its Fourier transform (b) |
|
(4) |
|
(5) |
其中，e－π(t－b)2为另一高斯函数曲线，a、b为不同时移因子(即噪声在时间域出现时间)，F(s)为其频谱值(频率域系数的绝对值).
图 4为该仿真信号时间域和频率域曲线，从中能够看到频谱不再是简单的正态分布函数，而是一个振荡的频谱，此时两组振荡因子相加取绝对值不再是简单的常数，而是一个复杂的振荡因子1+e－2πi(b－a)s，其绝对值为一个振荡函数，所以其频谱出现锯齿特征.
本文另仿真了一个高斯函数曲线信号和一个方波信号的复合噪声信号(见图 5)，其频谱计算结果见式(6)和(7).
|
图 5 一个高斯函数和一个方波函数的时间域(a)和频率域(b)特征 Fig. 5 Signal with one Gaussian function and one square wave in time domain and its Fourier transform |
|
(6) |
|
(7) |
其中，Π(t－b)为方波函数曲线，a、b为时移因子(即噪声在时间域出现时间)，F(s)为其频谱值(频率域系数的绝对值).
图 5为该噪声仿真信号在时间域和频率域的曲线特征，其频率谱已不再是正态分布函数和sinπs/πs函数的简单叠加，而同样是一个振荡的频谱.
对于任何一个复杂非周期噪声的傅里叶变换结果，可以理解成多个简单噪声傅里叶变换以后的线性相加结果.每一个噪声单元可以近似认为是一个振荡因子和一个能量衰减函数的乘积，其频谱曲线(能量曲线)是取绝对值后的结果，振荡因子与噪声出现的位置直接相关.本文只假设了几个最简单的噪声，频谱已非常复杂；如噪声变得更加复杂，存在多个非周期噪声叠加时，频谱曲线将会变得异常复杂(图 1)，频谱的振荡特征很难找到特定规律，所以在频率域直接计算噪声在CSEM频率位置的频率系数几乎是不可能的.
2.2 频谱包络虽然无法获得噪声在CSEM频率的准确系数，但是如果能够估计噪声的最大影响幅度(上界)，则可以评价CSEM频率信息的可靠程度.本文通过包络方法估计噪声在不同CSEM频率位置系数绝对值的上界，并计算了不同CSEM频率信息的噪声评价系数，以其评估每个CSEM频率信息的可靠性.
本文估计对象为噪声在CSEM频率的最大影响幅度，但原始数据在CSEM位置的频谱是噪声与有效信号的向量叠加，所以在对噪声频谱估计之前，应首先对噪声在CSEM频率的幅值进行预处理，本文将CSEM频率位置频谱设为左右相邻频率幅值的均值作为初始值，从而得到一组预处理频谱曲线(Modified Spectrum)，称为预处理频谱，包络算法是针对此预处理频谱进行的，预处理频谱计算方式见式(8).
|
(8) |
其中，
图 6为预处理频谱图，相比原始频谱只在有效频率位置进行改变(为相邻两个频率频谱平均值)，其他频率位置频谱保持不变，其中红色线为预处理频谱，黑色线为原始频谱.
|
图 6 原始信号频谱与预处理频谱 Fig. 6 Raw spectrum and modified spectrum |
本文采用基于希尔伯特变换的解析包络方法获得噪声频谱的上包络线.对一个实数振荡信号进行希尔伯特变换后，以该结果作为虚部，原始信号为实部，构成一个解析信号，而该解析信号的幅度值曲线，即为该实数信号的上包络线，在很多情况下被称之为瞬时振幅(式9)：
|
(9) |
其中，f为实数信号，fA即为f的解析信号，H(f(s))为实数信号的Hilbert变换结果(Bracewell, 2000)，|fA|即为对应实数信号的上包络曲线，其中：
|
(10) |
图 7为一组仿真信号的解析包络结果，在任意位置的信号大小都小于上包络线数值.
|
图 7 基于希尔伯特变换的解析包络 Fig. 7 Example of analytic envelope based on Hilbert transform |
严格意义上，基于希尔伯特变换获得包络线的方式是针对窄带信号进行的.频率谱并不是典型的振荡曲线，但是通过小波分解(Mallat, 1989; Meyer, 1993)，频谱曲线可以被分成低频部分(L)和高频部分(H)，高频部分可以近似认为在局部符合希尔伯特变换条件.通过解析包络得到高频部分的上界，则可以获得整个曲线的上包络线(式(11))，低频部分与高频部分上界相加重构即为整个频谱曲线的上界.
|
(11) |
其中，f为预处理频谱，L为其低频部分，H为其高频部分，Hsup为高频部分上界，fsup为整体频谱上界.
对频谱曲线进行离散小波变换，将预处理频谱曲线分解成低频部分与高频部分.小波分解层数与包络概率和精度有关，分解层数越多，包络概率越大，但是包络误差也越大，包络精度越差.这实际上是一个优化问题(Boyd, 2004)，即在包络比例和包络误差之间寻找一个最优比例关系，获取最优分解方式.根据仿真试验，发现对频率谱进行三级小波分解后应用包络算法，既可以获得较高的包络比例(包络比例约为80%左右)，又可以保证较小的误差值.
所以本文对预处理频谱应用三级小波分解，图 8为三级小波分解示意图.选择具有最大消失矩的DB小波进行小波分解，低频部分为变化趋势(A3)，剩余部分为高频振荡部分(D1+D2+D3)，对仿真信号三级小波分解后，低频和高频部分见图 9b和9c.
|
图 8 三级小波分解示意图 Fig. 8 Third level wavelet decomposition |
|
图 9 对噪声频谱进行三级小波分解 (a)预处理噪声频谱曲线(S); (b)低频部分(A3); (c)高频部分(D1+D2+D3). Fig. 9 Wavelet decomposition on noise spectrum (a) Modified spectrum (S); (b) Low frequency part (A3); (c) High frequency part (D1+D2+D3). |
对高频部分(D1+D2+D3)进行基于Hilbert变换的解析包络，获得其上包络线，包络结果见图 10a.获得高频部分的上包络曲线后与低频部分重构，得到整个频谱曲线的上包络曲线，见图 10b，包括原始数据的频谱曲线，噪声预处理频谱曲线，以及噪声频谱上包络线.
|
图 10 (a) 对高频部分进行解析包络并获取上包络线；(b)高频部分包络线与低频部分重构得到预处理频谱的整体包络线 Fig. 10 (a) Analytic envelope of high frequency part; (b) Total envelope for modified spectrum by reconstruction of analytic envelope and low frequency part |
以对数坐标绘制完整的频谱曲线，见图 11a(即对图 10b对数坐标显示)，主频及谐波成分共100个频率(主频为1 Hz，谐波为3~199 Hz)，图中包括原始数据的频谱曲线，噪声预处理频谱曲线，以及噪声频谱上包络线.
|
图 11 (a) 噪声频率谱包络结果；(b)噪声评价系数曲线 Fig. 11 (a) Envelope for modified spectrum in semilog coordinate; (b) Noise rating number curve |
通过上包络线获得在CSEM频率受噪声影响的上界值，与CSEM频率实际频谱的比值即为该频率的噪声评价数，计算公式如下(噪声评价数曲线见图 11b所示)：
|
(12) |
其中，Vf为噪声在有效频率位置的影响最大值，Sf为在有效频率的总频谱值，kf即为噪声评价数.
数值越小，信噪比越高.如果以噪声评价数5%为阈值，小于该阈值为合格数据，100个频率中将有93个频率是符合要求的.对此仿真信号，与只采用主频或者几个低阶谐波信号相比，本方法可提取更多有效频率信息.
3 包络分析在进行包络时，噪声在CSEM频率的频谱影响是未知的，频谱预处理过程中，对其进行的是均值处理，包络算法的本质为根据相邻频率的振荡情况估计频谱上界，并不能保证包络值一定大于噪声在CSEM频率位置的实际幅值，但是可以根据仿真信号包络结果统计包络比例.同时为了获得更高的包络比例，可以在解析包络的基础上，再进行一次极值包络.针对图 1中仿真信号频谱进行包络估计，对包络结果进行了统计：1)应用解析包络后，估计值大于实际值的比例是82%；2)再应用一次极值包络后，包络比例达到了90%，获得了更大的包络比例.下面对包络统计情况进行详细说明.
图 12a为对高频部分进行包络算法后的包络线特征，其中蓝色线为解析包络结果，而红色线为联合解析包络和极值包络的结果，其中极值包络是针对解析包络线进行的，目的是获得更大包络比例.
|
图 12 (a) 预处理频谱高频部分包络情况示意图；(b)包络结果对比图：(b1)包络值与真实值的差值，大于0说明包络值大于真实值，小于0说明包络值小于真实值，(b2)当差值小于零时，包络值与真实噪声值比例.蓝色线为解析包络结果，红色线为联合解析包络和极值包络结果 Fig. 12 (a) Envelope for high frequency part of modified spectrum; (b) Comparative result: (b1) Difference between envelope and true value, (b2) Ratio of envelope value to true value. Blue curve is the result of analytic envelope, while red curve is the result of analytic and peak envelope |
图 12b1为在CSEM频率位置包络曲线数值与噪声实际频谱的差值，如果差值大于0，说明包络值大于真实值；如果小于0，说明包络值小于真实值，图 12b2为包络值小于噪声实际频谱时(差值小于0)，包络值与噪声在CSEM频率真实幅值的比值.图 12b中包含了只应用解析包络(蓝色)和联合解析包络和极值包络(红色)的包络结果及对比图.
对仿真信号应用解析包络后，100个CSEM频率的包络值中，18个频率包络值比噪声真实值小，另外82个频率包络值大于噪声真实值，包络率为82%；而在18个包络值小于实际频谱的频率中，14个频率的估计噪声频谱与实际频谱的比例(估计值/实际值)为70%以上.如果应用解析包络后，再进行一次极值包络，包络比例上升为90%，90个CSEM频率的包络值大于噪声真实值；10个包络值小于真实值的频率中，9个的频率估计噪声频谱与真实值的比例大于70%, 包络值能很好的反映噪声在CSEM的影响幅度.实际上，需要根据实际情况进行包络次数的选择：如果更关注上界包络的绝对性，那么可以选择更多的包络次数；而如果更关注上界大小的准确性，则选择较少的包络次数，甚至只做一次解析包络，即可以直接输出估计结果.
为了确认包络比例具有一般性，本文对仿真信号中只存在高斯噪声时的频谱进行了包络估计，对预处理频谱应用三级小波分解，其低频部分和高频部分见图 13，其高频部分包络结果见图 14a，其中，黑色线为高斯噪声的频谱曲线，蓝色线为解析包络曲线，红色线为在解析包络基础上，再进行一次极值包络的结果.
|
图 13 对高斯噪声频谱进行三级小波分解 (a)噪声预处理频谱曲线; (b)低频部分(A3); (c)高频部分(D1+D2+D3). Fig. 13 Wavelet decomposition on Gaussian noise spectrum (a) Modified spectrum (S); (b) Low frequency part (A3); (c) High frequency part (D1+D2+D3). |
|
图 14 (a) 高斯噪声频谱高频部分包络情况示意图频谱包络情况示意图；(b)针对高斯噪声包络结果：(b1)包络值与真实值的差值，大于0说明包络值大于真实值，小于0说明包络值小于真实值，(b2)差值小于零时，计算得到的包络值与真实噪声值比例.蓝色线为解析包络结果，红色线为联合解析包络和极值包络结果 Fig. 14 (a) Envelope for high frequency part of modified spectrum; (b) Comparative result: (b1) Difference between envelope and true value, (b2) Ratio of envelope value to true value. Blue curve is the result of analytic envelope, while red curve is the result of analytic and peak envelope |
图 14b1为包络曲线值与噪声实际频谱的差值，其中统计结果如下：在100个CSEM频率，只进行解析包络后84个实际频谱小于包络值，16个频率大于包络值，比例约为84%；如果进行解析包络和一次极值包络，93个频率的包络值大于噪声真实值，包络比例约为93%.图 14b2为包络值小于实际噪声的部分包络值与实际值的比例关系，可以发现虽然有些频率包络值小于真实值，但是几乎都大于真实值的60%，与真实值具有相同尺度，所以在这些频点的包络值仍然能够表征噪声的最大影响幅值.
在此基础上，本文针对同一信号分别加入50组不同的随机噪声，应用包络算法估计，对其在CSEM频率位置包络比例情况进行了统计，统计结果见表 1.
|
表 1 对频谱曲线应用基于解析包络和极值包络算法后包络比例情况 Table 1 Envelope ratio based on analytic envelope and peak envelope |
根据统计结果，只应用解析包络时，平均包络比例为80.6%，应用解析包络和一次极值包络时，平均包络比例能够达到89.24%.图 15为50次包络测试包络比例的变化情况，包括解析包络比例和解析包络+极值包络比例及其平均值.
|
图 15 针对50组高斯噪声频谱曲线应用基于解析包络和极值包络后包络比例变化曲线 Fig. 15 Envelope ratio test based on 50 sets of Gaussian noise |
综上，应用基于希尔伯特变换的解析包络算法能够很好的估计噪声对CSEM频率幅值的影响上界值.通过应用基于希尔伯特变换的解析包络方法，能够获取非周期噪声在CSEM频率的上界值，利用该上界估计值可以计算CSEM频率位置的噪声评价数，根据该评价参数可筛选信号中符合规范的有效频率成分.
4 应用实例应用该包络方法对某实测2n序列伪随机信号数据进行了有效频率成分挖掘，其中2n序列伪随机信号不同频组主频成分见表 2.
|
|
表 2 2n序列伪随机7频波各频组主频频率(Hz) Table 2 Main frequencies for 7 frequency group based on 2n pseudo random signal |
图 16为2n序列伪随机信号7-4频组对应时域波形及其频率谱.当发射电流值为100 A时，在频谱中存在大量频率(包括7个主频和其谐波)对应强电流值，其中电流值大于35 A的频率7个，大于10 A的频率19个，大于5 A的频率55个，大于1 A的频率数目235个，具有大量潜在有效频率成分.2n序列伪随机信号的其他频组均同样具有类似频谱特征.
|
图 16 2n序列伪随机信号7频波4频点时域波形(a)及频谱(b) Fig. 16 Time domain (a) and frequency domain (b) of 7-2 frequency group |
图 17为江西某地7-4频组野外实测电流数据频谱图，本文对所有频率(主频及谐波成分)的电流数值进行统计，电流大于1 A的频率个数为402个.就电流绝对大小而言，信号中大量频率(包括大量谐波)具有较大电流值，所以在原始信号中存在大量潜在可用信息.
|
图 17 发射端实测电流频谱图 Fig. 17 Spectrum of measured current from field word |
对原始接收数据应用本文提出的包络方法能够提取大量高信噪比频率信息，下面以该测区187点数据为例展开详细说明.
图 18c1为对187点7-4频组原始数据频谱应用解析包络和极值包络后的结果，其中蓝色线为原始频谱，黑色线为预处理频谱上包络曲线，获得在CSEM频率位置的包络值后，计算噪声评价数.图 18c2为所有CSEM频率的噪声评价数曲线，很多谐波的噪声评价数小于5%.对于7-4频组进行有效信号筛选，其中噪声评价数小于5%的频率个数为135个，评价数小于3%的频率个数为80个，相对于只提取7个主频成分，本方法提取了更多的有效频率信息.
|
图 18 (a) 187点7-2频组：(a1)有效信号频率谱包络结果，(a2)不同频率对应评价数；(b) 187点7-3频组：(b1)有效信号频率谱包络结果，(b2)不同频率对应评价数；(c) 187点7-4频组：(c1)有效信号频率谱包络结果，(c2)不同频率对应评价数；(d) 187点7-5频组：(d1)有效信号频率谱包络结果，(d2)不同频率对应评价数 Fig. 18 (a) Result of 187 point 7-2 frequency group: (a1) Envelope result, (a2) Noise rating curve; (b) Result of 187 point 7-3 frequency group: (b1) Envelope result; (b2) Noise rating curve; (c) Result of 187 point 7-4 frequency group: (c1) Envelope result; (c2) Noise rating curve; (d) Result of 187 point 7-5 frequency group: (d1) Envelope result, (d2) Noise rating curve |
对该测点的7-2频组、7-3频组及7-5频组数据分别应用包络算法进行有效数据提取，包络曲线和噪声评价数曲线分别见图 18a、18b及18d.其中，图(a1, b1, d1)为频谱包络结果，图(a2, b2, d2)为根据包络结果计算得到的噪声评价数.
图 18a1中，工频噪声(50 Hz)对伪随机7-2频组信号频谱造成了影响，图 18a2中在50 Hz附近几乎没有噪声评价数小于5%的频率，在图 18b中伪随机7-3频组同样存在类似的情况，因为受工频干扰在50 Hz附近的勘探频率受到的强干扰，所以其噪声评价系数很大.
汇总4个频组中噪声评价数据，并分别以3%误差和5%误差为阈值进行了有效频率提取，图 19中蓝色点为噪声评价数小于5%的频率归一化电场信号，黑色点为主频，48 Hz明显受到工频干扰，通过本文方法的筛选，受到工频干扰严重的48 Hz信号被剔除，且在32 Hz与64 Hz之间，提取了部分可用谐波成分，提高了该频率段的频率密度；图 20中红色点为3%误差阈值的提取结果后对其归一化结果，黑色点为主频.
|
图 19 以5%噪声评价数为阈值得到的结果(黑色点为主频信号分布位置，对应区间内共23个主频频率，蓝色点为提取后有效值，共257个频率) Fig. 19 Filter result based on threshold of 5% (black points are main frequencies, totally 23 frequencies; blue points are extracted frequencies, totally 257 frequencies) |
|
图 20 以3%噪声评价数为阈值得到的结果(黑色点为主频信号分布位置，区间内共21个主频频率，红色点为提取后有效频率位置，共110个频率) Fig. 20 Filter result based on threshold of 3% (black points are main frequencies, totally 21 frequencies; red points are extracted frequencies, totally 110 frequencies) |
在频率区间内，主频共有23个频率，而以5%噪声评价数阈值对数据进行有效信息提取，能够获得257个频率，以3%噪声评价数阈值提取能够获得110个频率，大大增加了可用有效频率个数，能够提高地球物理反演的纵向分辨率.在实际应用中，可根据不同的实际需求选择不同的阈值进行有效信息筛选.
5 结论(1) 通过应用本文所提出的噪声评价方法，能够在频率域快速的提取CSEM信号中有效频率成分.在不增加任何野外工作量的基础上，能够提取更多的有效频率信息.
(2) 此种方法还可以提高野外的勘探效率，由于采集低频信号中含有高频谐波成分，如果此时所含高频成分已经符合勘探的要求，不必要再额外采集高频信号，能够提高野外采集效率，节省勘探成本.
(3) 此方法可以应用于所有频率域电磁勘探信号的有效频率成分提取.
Boashash B.
1992. Estimating and Interpreting the Instantaneous Frequency of a Signal I:Fundamentals. Proceedings of the IEEE, 80(4): 519-538.
|
|
Boyd S, Vandenberghe L. 2004.
Convex Optimization. Cambridge: Cambridge University Press.
|
|
Bracewell R. 2000.
The Fourier Transform and Its Applications. (3rd ed). New York: McGraw-Hill.
|
|
Daubechies I. 1992. Ten lectures on wavelets.//CBMS-NSF Conference Series in Applied Mathematics. Philadelphia:SIAM.
|
|
Feldman M.
2011. Hilbert transform in vibration analysis. Mechanical Systems and Signal Processing, 25(3): 735-802.
DOI:10.1016/j.ymssp.2010.07.018 |
|
He J S.
2007. The new development of frequency domain electro-prospecting. Progress in Geophysics, 22(4): 1250-1254.
|
|
He J S.
2010a. Closed addition in a three-element set and 2n sequence pseudo-random signal coding. Journal of Central South University (Science and Technology), 41(2): 613-637.
|
|
He J S. 2010b.
Wide Field Electromagnetic Method and Pseudo Random Signal Electrical Prospecting. Beijing: Higher Education Press.
|
|
Li J M. 2005.
Geoelectric Field and Electrical Exploration. Beijing: Geological Publishing House.
|
|
Li J, Tang J T, Xu Z M.
2017. Magnetotelluric noise suppression base on signal-to-noise indentification in ore concentration area. Chinese Journal of Geophysics, 60(2): 722-737.
DOI:10.6038/cjg20170224 |
|
Mallat S G.
1989. A theory for multiresolution signal decomposition:the wavelet representation. IEEE Transactions on Pattern Analysis and Machine Intellignece, 11(7): 674-693.
DOI:10.1109/34.192463 |
|
Marple S L.
1999. Computing the discrete-time analytic signal via FFT. IEEE Transactions on Signal Processing, 47(9): 2600-2603.
DOI:10.1109/78.782222 |
|
Meyer Y. 1995.
Wavelets and Operators. Cambridge: Cambridge University Press.
|
|
Picinbono B.
1997. On instantaneous amplitude and phase of signals. IEEE Transactions on Signal Processing, 45(3): 552-560.
DOI:10.1109/78.558469 |
|
Shensa M J.
1992. The discrete wavelet transform:wedding the a trous and Mallat algorithms. IEEE Transactions on Signal Processing, 40(10): 2464-2482.
DOI:10.1109/78.157290 |
|
Singleton R C.
1969. An algorithm for computing the mixed radix fast Fourier transform. IEEE Transactions on Audio and Electroacoustics, 17(2): 93-103.
DOI:10.1109/TAU.1969.1162042 |
|
Tang J T, He J S. 2005.
Controllable Source Audio Magnetotelluric Method and Its Application. Changsha: Central South University Press.
|
|
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, Liu Z J, Liu F Y, et al.
2015a. The denoising of the audio magnetotelluric data set with strong interferences. Chinese Journal of Geophysics, 58(12): 4636-4647.
DOI:10.6038/cjg20151225 |
|
Tang J T, Ren Z Y, Zhou C, et al.
2015b. Frequency-domain electromagnetic methods for exploration of the shallow subsurface:A review. Chinese Journal of Geophysics, 58(8): 2682-2705.
DOI:10.6038/cjg20150807 |
|
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 |
|
Vakman D E.
1972. On the definition of concepts of amplitude, phase and instantaneous frequency of a signal. Radio Eng. Electron. Phys., 17(5): 754-759.
|
|
Welch P D.
1967. The use of fast Fourier transform for the estimation of power spectra:a method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2): 70-73.
DOI:10.1109/TAU.1967.1161901 |
|
Xi Z Z, Bao G S.
1996. A tentative discussion on hacsamy method. Geophysical and Geochemical Exploration, 20(5): 397-400.
|
|
Zonge K L, Hughes L J. 1981.
The Complex Resistivity Method, in Advances In Induced Polarization and Complex Resistivity. Tucson: University of Arizona Press.
|
|
何继善.
2007. 频率域电法的新进展. 地球物理学进展, 22(4): 1250–1254.
|
|
何继善.
2010a. 三元素集合中的自封闭加法与2n序列伪随机信号编码. 中南大学学报(自然科学版), 41(2): 613–637.
|
|
何继善. 2010b.
广域电磁法和伪随机信号电法. 北京: 高等教育出版社.
|
|
李晋, 汤井田, 徐志敏, 等.
2017. 基于信噪辨识的矿集区大地电磁噪声压制. 地球物理学报, 60(2): 722–737.
DOI:10.6038/cjg20170224 |
|
李金铭. 2005.
地电场与电法勘探. 北京: 地质出版社.
|
|
汤井田, 何继善. 2005.
可控源音频大地电磁法及其应用. 长沙: 中南大学出版社.
|
|
汤井田, 李晋, 肖晓, 等.
2012a. 数学形态滤波与大地电磁噪声压制. 地球物理学报, 55(5): 1784–1793.
|
|
汤井田, 徐志敏, 肖晓, 等.
2012b. 庐枞矿集区大地电磁测深强噪声的影响规律. 地球物理学报, 55(12): 4147–4159.
DOI:10.6038/j.issn.0001-5733.2012.12.027 |
|
汤井田, 刘子杰, 刘峰屹, 等.
2015a. 音频大地电磁法强干扰压制试验研究. 地球物理学报, 58(12): 4636–4647.
DOI:10.6038/cjg20151225 |
|
汤井田, 任政勇, 周聪, 等.
2015b. 浅部频率域电磁勘探方法综述. 地球物理学报, 58(8): 2682–2705.
DOI:10.6038/cjg20150807 |
|
席振铢, 鲍光淑.
1996. 试论HACSAMT法. 物探与化探, 20(5): 397–400.
|
|