2014, Vol. 57
2. 中国国土资源航空物探与遥感中心, 北京 100083;
3. 成都理工大学地球探测与信息技术教育部重点实验室, 成都 610059
2. China Aero Geophysical Survey & Remote Sensing Center for Land and Resources, Beijing 100083, China;
3. Key Lab of Earth Exploration & Information Techniques of Ministry of Education, Chengdu University of Technology, Chengdu 610059, China
航空电磁方法具有勘探面积大、勘探效率高、勘探成本相对较低等优点，近几十年来，已经获得了快速的发展(Palacky et al., 1991; 张昌达, 2006).相对频率域航空电磁法而言，时间域航空电磁法勘探深度大、能获得更为精细的地下介质的电阻率分布，因此时间域航空电磁法广泛应用于地质填图、地下 水调查、矿产勘探等(Jorgensen et al., 2003; Smith, 2010; 朴化荣等, 1980; 黄皓平等, 1990).数据质量是获得好的应用效果的基本保证(Sorensen et al., 2004)，但是由于时间域航空电磁法是一种宽带电磁勘探方法，为此，其信号往往受到各种噪声的污染和干扰，所以，在不能大幅度地提高仪器硬件系统的信噪比之前，对去噪处理方法的研究就显得极为重要.
在航空电磁噪声的来源及特性方面，已有许多学者进行了详细研究并对电磁噪声的频谱进行了详细分析(Macnae et al., 1984; McCracken et al., 1984; Szarka, 1988).时间域电磁数据的频率范围大约在0.01 Hz至几十千赫兹的范围，其中低于1Hz的电磁噪声主要来自太阳散射的等离子和地球电场之间交互作用的结果，并且这些噪声的幅度和频率近似成反比；频率高于1Hz的电磁噪声主要来源大气中的雷暴活动，常常看作是天电噪声(Spies, 1988).实际上，航空时间域电磁资料的噪声不仅仅包含天然噪声，而且还包含工业电流干扰、建筑、管线等人工设施引起的人文噪声(Szarka, 1988).
噪声抑制不仅可以提高数据的信噪比，为提高 反演和解释精度奠定良好的基础，同时还可以在 不提高激发功率的条件下增大勘探深度(Spies, 1988). 为此，许多研究人员提出了多种不同的噪声抑制方法.到目前为止，所提出的噪声抑制方法可以分为三类：抑制或剔除尖脉冲或振荡的滤波器 (Macnae et al., 1984; Bouchedda, 2010)；基于预 测方法的去噪方法(Spies, 1988; Buselli, 1996; Buselli, 1998)；通用的叠加方法(Macnae et al., 1984; McCracken et al., 1984; Cull, 1991). 前两类去噪方法只能抑制或剔除某一种特定类型的噪声，而第三类去噪方法则只能抑制掉部分的天电噪声，但是对于人文干扰等噪声则不能获得较好的抑制效果(Reninger et al., 2011).此外，还有最小二乘拟合去噪法(程德福, 2003)、独立主成分分析法(刘祥平等, 2011)、小波去噪法(邵敏等, 2008)等.由于航空瞬变电磁噪声幅值较大，所以最小二乘拟合法往往很难收敛；刘祥平等(2011)所给出的独立主成分分析法主要用于去掉工频干扰；从邵敏等(2008)的处理结果看，尖脉冲干扰尚未获得很好的抑制.
为了能获得一种较为通用的去噪方法，Reninger等(2011)提出了一种基于奇异值分解的时间域航空电磁噪声抑制方法.该方法在抑制人文设施的耦合噪声和天电噪声方面获得了很好的效果.实际上，奇异值分解就是统计学中的线性主成分分析方法.从统计角度分析，该方法只考虑了数据的一阶和二阶特性，尚未考虑数据的高阶特性，只能描述平稳数据，而瞬变电磁数据往往是非平稳的，因此该方法的重构精度需要进一步的改善和提高.
为了能描述非平稳数据，许多学者提出了不同的非线性主成分分析方法，使得非线性主成分分析方法得到了极大的发展(Jolliffe, 2002).在非线性主成分分析方法中，以核主成分分析方法(Kernel Principal Component Analysis，KPCA)(Schölkopf et al.,1998)的思想最为新颖.该方法的提出极大丰富了非线性主成分分析和统计分析方法的内容.具体来说，核主成分分析方法利用了核方法的性质，即通过一个非线性映射函数，把在原始空间的线性主成分分析方法推广到了高维特征空间(即再生核希尔伯特空间)中.相对原始空间来说，它是非线性的.本文的工作就是对核主成分分析方法进行研究并设计出基于核主成分分析的电磁去噪方法，然后使用野外实测数据检验该方法的有效性.
设x={xk,k=1,2…,M}为一组随机向量，其中xk∈RN，M为随机向量个数，N表示随机向量的维数，并且
，则x的协方差矩阵可表示为
然后，利用如下等式
对C进行特征值分解和排序可以得到随机向量x的特征值λi和特征向量wi,其中p表示特征值的个数.于是可以得到随机向量x的p个主方向wi.
现在通过一个非线性函数φ把输入空间映射到特征空间：
不失一般性，假定在特征空间中
，得到在特征空间F中的协方差矩阵为
对该协方差矩阵进行特征值分解，即
式中，所有对应于特征值 λ F≠0的特征向量WF都处于φ(x1),…,φ(xM)所张成的空间中.因此，得到如下等式：
其中，
把(7)式代入(6)式，得到：
定义一个M×M矩阵K(xi,xj)为
则(8)式可以表示为
式中，α为α1,α2,…,αM的列向量.实际上要解(10)式，也就是要求下面等式的特征值和特征向量的问题.
对于(11)式，如果用λ1≥λ2≥…≥λM表示矩阵K的特征值，那么相应的α1,α2,…,αM就是矩阵的特征向量.要在高维特征空间归一化特征向量WF，即假设λkF≠0，k=1,2,…,l.利用(7)、(9)和(11)式，(WkF·WkF)=1就等价于
如果要提取某个点的主成分，必须计算该点在特征向量WkF上的投影(在特征空间F中).现在假定x为输入空间中的一个测试点，它在特征空间F中的像为φ(x)，得到
这就是所要求的主成分.
在核主成分分析方法中，核函数的选择较为关键，因为其决定数据的相似性度量.根据Hilbert-Schmidt理论，K(x,y)可以是满足Mercer定理的任意对称函数(Aron szajn,1950).目前，常用的核函数主要有多项式核函数K(x,y)=(x·y+1)d，高斯核函数
，神经网络核函数
等.在实际应用中，到底选择什么样的核函数呢？到目前为止，尚未有核函数选择的指导性理论，只能通过试验进行选择.幸运的是，针对某一类较为宽泛的数据，选定的核函数一般都具有较好的普适性.
设xi为第i个测点的衰减曲线，那么基于核主成分分析的去噪算法可以描述如下：
(1)选择一个合适的核函数计算矩阵K(xi,xj)；
(2)利用式(11)计算特征空间中的特征值λkF和特征向量WkF；
(3)在特征空间中对特征向量WkF进行归一化；
(4)对特征向量进行分析，选择反映地下介质的成分进行重构.
线性主成分分析方法，只利用了二阶统计信息，而不能获得数据的高阶特征，忽略掉了数据的非线性信息.而核主成分分析方法作为一种非线性分析工具，能获得数据的高阶统计特征，从而在分析过程中获得较为丰富的信息.航空瞬变电磁数据往往是非线性的，如果采用线性主成分分析方法进行分析，那么将会忽略掉许多丰富的细节，导致重构精度降低.反之，如果采用核主成分分析方法，因为该方法能充分获得数据的非线性信息，从而对数据分解更为精细，并获得较高的重构精度和更好的去噪效果.
为了说明核主成分分析方法的计算结果，模拟包含噪声的航空瞬变电磁二次场数据，其抽道后的归一化单点衰减曲线如图1所示，其中横坐标为时间，纵坐标为感生电动势.
 | 图1 含噪的模拟衰减曲线 Fig.1 A simulating decay curve normalized with noise |
在应用核主成分分析方法对含噪正演模拟数据进行分析时，为了简化计算，只模拟了14个测点的衰减数据或曲线.使用式(11)对这些含噪声数据进行分析，获得14个主成分如图2所示.
 | 图2 对含噪的模拟数据进行分析后获得的14个成分 Fig.2 Plots of 14 components of simulating data |
图2所显示的14个主成分是根据特征值或者成分的幅度值按降序排列，实际上这些成分也描述了其在重构时的贡献程度.从图2所示的实验结果，发现第一成分表示曲线的衰减趋势，而对应较小特征值的成分则描述了曲线的细节变化.因此，振荡和 尖脉冲可以用高阶分量解释，即高阶分量看作是噪声信号，而低阶分量可看作是反映地下介质的有效信号.
在“十一五”期间国家863计划课题组完成了直升机时间域电磁系统(CHTEM)的飞行测量试验， 并在国内某重点勘探区获得了实测数据.图3显示了野外实测但是尚未进行叠加和抽道的原始数据，它不仅包含了天电噪声，而且还包含人文噪声、随机噪声、运动噪声等.为了能获得较好的数据质量、更好的反演精度和解释结果，前期的去噪处理方法研究和应用非常关键.
 | 图3 直升机时间域电磁系统野外实测单点衰减曲线 Fig.3 A raw curve decay measured in site |
为了获得较好的处理效果，应用核主成分分析方法对时间域航空电磁数据进行去噪处理的步骤如下：(1)根据系统的采样频率特点，为了避免更高频率的噪声干扰，对原始未叠加数据进行低通滤波，滤除非有效信号频率带的噪声；(2)对低通滤波后的电磁数据进行叠加；(3)对叠加后数据进行抽道；(4)应用核主成分分析方法对抽道后的数据进行分析并重构，从而获得抑制噪声后的时间域航空电磁剖面数据；(5)对抑制噪声后的时间域航空电磁剖面数据进行漂移校正和基线校正，突出异常信息.在这些处理步骤中，叠加和抽道可以抑制大部分随机噪声；而核主成分滤波方法则主要抑制天电、人文以及部分运动噪声.
在使用核主成分分析方法进行滤波时，哪些成分是反映地下介质的电磁信号？哪些成分代表噪声？选择多少成分进行重构？对这些问题的回答成为处理结果好坏的关键.根据上述仿真数据的成分分析结果，可以认为主分量代表有效电磁信号，而微小分量代表噪声，那么从数学上看，可以采用能量占比方法，即所选择的主成分的能量占总能量的多少来确定主成分数量，例如采用如下公式进行计算：
其中，E表示能量占比，一般可选择能量占比为90%~95%；λi表示按降序排列后的第i个特征值；p为能量占比达到要求后的特征值数；N为所有的特征值数.当然，不同的勘探区需要经过多次试验才能获得合理的参数.Reninger等(2011)采用特征值加权和飞行高度相结合的方法来确定噪声分量和地质分量，但是这种方法较为繁琐.
实际上，在勘探区进行飞行测量时，往往为了校正漂移量，一般都需要进行高空基线飞行,在测线上飞行测量的电磁数据往往包括两部分：勘查数据和背景场数据.在应用核主成分分析方法分析时间域航空电磁数据时，为了判别哪些成分代表噪声以及哪些成分代表地质响应，可以分别对勘查数据和背景场数据进行分析，然后对这些成分的相似性进行判别，最后选择相似性大的成分进行重构.不同地区的外部电磁环境是有一定差异的，因而电磁场的噪声信息也不尽相同.但是，在范围不很大的勘探区，其噪声具有一定的趋同性，这种相似性度量要根据该区域的实际测量数据来确定.
采用本文所推荐的核主成分分析方法对实测时间域航空电磁数据进行了处理，获得了较为满意的效果，其中所用的核函数为多项式核.图4为在某地勘探区2070测线的瞬变剖面数据，该数据已经经过低通滤波、叠加、抽道和基线校正.从图4可以看到，数据包含较多的尖脉冲和振荡，这些噪声可能是雷电或其他外部干扰引起的结果.
 | 图4 某勘探区2070测线瞬变剖面数据 Fig.4 Electromagnetic data section in line 2070 |
对图4所示的电磁数据使用核主成分分析方法滤波后，获得了良好的处理效果(参见图5和图6).为了说明本文所推荐滤波方法的处理效果，我们对基于核主成分分析的滤波处理结果和AeroTEM商用软件的处理结果进行了对比，其中滤波处理后的异常总体形态和幅值比较见图5.从图5所示的核主成分滤波方法和AeroTEM 软件处理结果的总体形态和幅度比较可以看到，采用核主成分滤波方法不仅天电等噪声得到了较好的剔除，而且人文噪声也得到了很好的抑制，在异常总体形态上，核主成分滤波器和AeroTEM软件在波形和保幅上都获得较好的效果，两种处理方法结果相当，应用核主成分滤波方法已达到国外先进航空瞬变电磁数据处理商用软件的结果，说明采用核主成分滤波方法是有效的.
 | 图5 滤波处理异常总体形态和幅值对比结果 (a)使用AeroTEM软件对2070测线数据的滤波处理结果;(b)使用核主成分滤波方法对2070测线数据的滤波处理结果. Fig.5 Comparison of filtering results in line 2070 (a) The filtering results using AeroTEM software; (b) The filtering results using the filtering method based on KPCA. |
 | 滤波处理结果对比(最后两道数据，可以看作噪声) (a)使用AeroTEM软件对2070测线数据的滤波处理结果(最后两道数据)；(b)使用核主成分滤波方法对2070测线数据的滤波处理结果(最后两道数据)． Fig.6 Comparison of filtering result in line 2070(last two data channels) (a)The filtering result using AeroTEM software;(b)The filtering result using filtering method on KPCA. |
此外，为了比较核主成分滤波方法和AeroTEM 软件在细节处理上的优劣，我们抽取滤波后的晚期道数据进行比较.图6为从图5所示的瞬变剖面中抽取最后两道数据(晚期道，可以看作噪声)的对比图.从对比结果，如图6中的五个圆圈所圈定的对比结果来看，核主成分滤波方法在抑制尖脉冲和振荡方面获得较好的效果，而AeroTEM软件处理后还剩余部分尖脉冲或振荡.因此，所推荐的滤波方法在精细处理上要优于AeroTEM软件.
在时间域航空电磁探测中，由于实测数据往往包含多种不同类型的噪声，仅仅依赖单一的频率域滤波器或者时间域滤波器都很难获得较好的处理质量，所以应该从统计分析的角度设计出更加实用的滤波器才能对这些电磁噪声进行描述和剔除.本文对核主成分分析方法进行了研究，并设计出基于核主成分分析方法的时间域航空电磁去噪方法，即核主成分滤波器.该滤波器对电磁数据进行了从低阶到高阶的统计分析，获得了不同阶的统计信息，对地质响应信号和噪声进行了有效分离，从而可以抑制背景噪声、天电噪声等，有效地提高了数据处理的可 靠性和质量.应用该滤波器对“十一五”期间实测的 数据进行处理，获得了很好的处理效果.
通过处理仿真数据和实际资料，获得如下结论：第一、在处理某区的实测数据时，必须对该数据进行特性分析；第二、和传统的线性去噪方法相比，核主成分滤波器由于能进行高阶统计分析，因此可以实现精细的数据分解和获得较高的重构精度，从而对时间域航空电磁噪声具有较好的抑制作用；第三、世界上没有万能滤波器，因此必须组合不同滤波器，才能抑制时间域航空电磁数据中所包含的多种不同类型噪声；最后，尽管核主成分滤波器在抑制天电和人文噪声方面获得较好的效果，但是尚不能有效地抑制运动噪声.为此，寻求运动噪声的特性或规律并研究有效的去噪方法成为未来的研究课题.
| [1] | Aron szajn N. 1950. Theory of reproducing kernels. Transactions of the American Mathematical Society, 68(3): 337-404, doi: 10.2307/1990404. |
| [2] | Bouchedda A, Chouteau M, keating P, et al. 2010. Sferics noise reduction in time-domain electromagnetic systems: application to megaTEMII signal enhancement. Exploration Geophysics, 41(4): 225-239, doi: 10.1071/EG09007. |
| [3] | Buselli G, Cameron M. 1996. Robust statistical methods for reducing sferics noise contaminating transient electromagnetic measurement. Geophysics, 61(6): 1633-1644, doi: 10. 1190/1.1444082. |
| [4] | Buselli G, Pi k J P, Hwang H S. 1998. AEM noise reduction with remote referencing. Exploration Geophysics, 29(2): 71-76, doi: 10.1071/EG998071. |
| [5] | Cheng D F. 2003. Study on signal detection of near zone magnetic source transient electromagnetic method [Ph. D. thesis]. Changchun: Jilin University. |
| [6] | Cull J P. 1991. Signal processing concepts for airborne SIROTEM data. Exploration Geophysics, 22(1): 97-100, doi: 10.1071/EG991097. |
| [7] | Huang H P, Wang W ZH. 1990. Inversion of time-domain airborne electromagnetic data. Chinese J. Geophys. , (in Chinese), 33(1): 87-97. |
| [8] | Jolliffe I T. 2002. Principal Component Analysis (2nd ed). New York: Springer. |
| [9] | Jorgensen F, Sandersen P B E, Auken E. 2003. Imaging buried quaternary valley using the transient electromagnetic method. Journal of Applied Geophysics, 53(4): 199-213, doi: 10.1016/j. jappgeo. 2003. 08.016. |
| [10] | Liu X P, Wang J M. 2011. Improved ICA denoising method and application in transient electromagnetic signal processing. Journal of Beijing Normal University (Natural Science), 47(1): 34-39. |
| [11] | 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. |
| [12] | 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. |
| [13] | McCr acken K G, Oristaglio M L, Hohmann G W. 1986. Minimization of noise in electromagnetic explorations. Geophysics, 51(3): 819-832, doi: 10.1190/1.1442134. |
| [14] | Palacky C J, West G F. 1991. Airborne electromagnetic methods. // Nabighian M N ed. Electromagnetic methods in applied geophysics (Vol. 02), Tulsa: Society of Exploration Geophysicists, 811-879. |
| [15] | Piao H R, Sha S Q, Wang Y L. 1980. Time-domain electromagnetic response above the surface of a homogeneous earth. Chinese J. Geophys. (in Chinese), 23(2): 207-218. |
| [16] | Reninger P-A, Martelet G, Deparis J, et al. 2011. Singular value decomposition as denoising tool for airborne time domain electromagnetic data. Journal of Applied Geophysics,75(2): 264-276, doi: 10.1016/j. jappgeo. 2011. 06.034. |
| [17] | Schölkopf B, Smola A, Müller K R. 1998. Nonlinear component analysis as kernel eigenvalue problem. Neural Computation, 10(5): 1299-1319, doi: 10.1162/089976698300017467. |
| [18] | Shao M, Qiu N, He Z X. 2008. Effect of wavelet threshold de-noising on long-offset transient electromagnetic signal. Chinese Journal of Engineering Geophysics, 5(1): 70-74. |
| [19] | Smith R. 2010. Airborne electromagnetic methods: applications to minerals, water and hydrocarbon exploration. CSEG 2010 distinguished lecture, 7-10. |
| [20] |
Sorensen K I, Auken E. 2004. SkyTEM-A new high-resolution helicopter transient electromagnetic system. Exploration Geophysics, 35(3): 191-199, doi: 10.1071/EG04194. |
| [21] | Spies B R. 1988. Local noise prediction filtering for central induction transient electromagnetic sounding. Geophysics, 53(8): 1068-1079, doi: 10.1190/1.1442543. |
| [22] | Szarka L. 1988. Geophysical aspects of man-made electromagnetic noise in the earth-a review. Surveys in Geophysics, 9(3-4): 287-318, doi: 10. 1007/BF01901627. |
| [23] | Zhang C D. 2006. Airborne time domain electromagnetic system: Look back and ahead. Chinese Journal of Engineering Geophysics (in Chinese), 3(4): 265-273. |
| [24] | 程德福. 2003. 近区磁源瞬变电磁法信号检测技术研究[博士论文]. 长春: 吉林大学. |
| [25] | 黄皓平, 王维中. 1990. 时间域航空电磁数据的反演. 地球物理学报, 33(1): 87-97. |
| [26] | 刘祥平, 王建明. 2011. 改进ICA去噪方法在瞬变电磁信号处理中的应用. 北京师范大学学报(自然科学版), 47(1): 34-39. |
| [27] | 朴化荣, 沙树琴, 王延良. 1980. 均匀大地上空的时间域电磁响应. 地球物理学报, 23(2): 207-218. |
| [28] | 邵敏, 邱宁, 何展翔. 2008. 长偏移距瞬变电磁信号小波阈值去噪效果分析. 工程地球物理学报, 5(1): 70-74. |
| [29] | 张昌达. 2006. 航空时间域电磁法测量系统: 回顾与前瞻. 工程地球物理学报, 3(4): 264-273. |