曹小玲 刘开元 严良俊

(①长江大学油气资源与勘探技术教育部重点实验室，湖北武汉 430100； ②非常规油气湖北省协同创新中心，湖北武汉 430100；③长江大学信息与数学学院，湖北荆州 434023； ④中国石油川庆钻探工程有限公司地球物理勘探公司，四川成都 610213)

1 引言

大地电磁测深(MT)是一种被动源电磁探测技术，通过在地表测量正交电、磁场分量的扰动值，研究地下介质的地电结构[1]。MT探测深度的范围从地表到上地幔，它既可以研究大尺度地质单元的电磁感应、深部地壳地幔电性结构及壳幔中的高导层，也可以开展油气资源普查、地热资源勘查、天然地震活动前兆监测以及地下水资源勘查和环境保护等，因而得到广泛应用[2]。但大地电磁测深观测的天然电磁场具有信号弱、频带宽等特点，极易受到各种噪声的干扰，因此大地电磁数据的去噪一直是研究热点。地球物理学家们提出一系列技术，这些技术主要集中在四个方面：远参考方法、稳健阻抗估计、时频转换技术和时间域去噪[3-7]。

传统的数据处理基于傅里叶变换—功率谱估计，利用傅里叶变换将信号通过一个低通或带通滤波器，映射到整个频域内加以分析处理，但无法表达信号的时域局部性质，而这恰恰是大地电磁信号作为一种非平稳信号的最基本和最关键的特征。为此，小波分析被陆续应用于大地电磁的去噪处理[8-13]。王书明等[14,15]根据MT信号的双相干系数、三相干系数以及在z平面上的零点分布，分析研究不同地区实测MT信号的特征，得出MT信号一般具有非高斯性、非线性和非最小相位特征,这个结论启发人们用新的信号处理方法解决大地电磁的噪声问题。近年来，随着盲源分离技术的迅猛发展，独立分量分析理论[16-19]被广泛应用于生物医学信号处理、语音信号处理、图像处理、地震勘探等多个领域，在信号去噪中表现优异。独立分量分析要求源信号在统计特性上具有非高斯性、非依赖性、非白性[16-19]，大地电磁信号恰好符合这些要求；另一方面，独立分量分析不需要知道信号的先验知识，即它只要通过观测信号就能推测甚至恢复原始信号，这对大地电磁测深来说意义重大，因为由于各种客观因素的制约，导致在大地电磁测深中想要获得相关的先验知识一般来说是非常困难的。因此，使用独立分量分析技术对大地电磁信号的各种信息成分进行分析和提取，尤其是对噪声成分进行分析和提取，符合大地电磁测深实际应用的特点和需要。独立分量分析处理的对象是一组相互统计独立的信号源经线性组合而产生的混合信号,最终从混合信号中提取各独立的信号分量。就大地电磁噪声干扰问题,目标信号和噪声干扰分别由不同的信号源产生，因此彼此相互独立。由统计原理，相互独立的随机变量必然不相关； 但反过来却不一定成立[20]，即不相关的随机变量不一定相互独立。独立分量分析方法在统计独立意义下对混合信号进行分离，以达到消除噪声的目的，这显然比基于不相关性度量的传统自适应滤波方法的效果更好。因此，根据大地电磁测深数据的特点和其噪声的来源，以及大地电磁资料的处理流程和解释流程，研究基于独立分量分析的大地电磁数据噪声消除方法，并通过实际大地电磁数据验证方法的有效性。

2 独立分量分析

独立分量分析(ICA)[16-19]是20世纪90年代发展起来的一种信号处理技术，它是从多维统计数据中找出隐含因子或分量的方法。从线性变换和线性空间角度来说，源信号为相互独立的非高斯信号，可以看作线性空间的基信号，而观测信号则为源信号的线性组合。ICA就是在源信号和线性变换均未知的情况下，从观测的混合信号中估计出数据空间的基本结构或源信号。ICA的算法流程如图1所示：在源信号s(t)中各分量相互独立的假设下，由观察信号x(t)通过解混系统B把他们分离开来，使输出信号y(t)逼近s(t)。

图1 ICA的算法流程图

2.1 数据预处理

ICA通常需要进行零均值化和白化预处理。

一般情况下，观测数据是具有相关性的，因此需要对零均值化后的数据进行白化处理，以去除相关性、简化后续分量的提取过程，数据的白化处理还能有效增强算法的收敛性。

若均值为零的随机向量Z=(Z1,…,ZM)T满足E(ZZT)=I(I为单位矩阵，M为随机向量中元素的个数，E(·)是期望值算子)，则称Z为白化向量。白化的本质是去相关。对观测信号X，需寻找一个变换，使X经变换后变成白化向量，即有：Z=QX,其中Q为白化矩阵。

零均值化和白化作为ICA的预处理可以有效地降低问题的复杂度。

2.2 FastICA算法

FastICA算法又称固定点(Fixed-Point)算法，它采用批处理方式，速度比较快,是一种快速寻优迭代算法。FastICA算法有多种形式，本文采用基于负熵最大的FastICA算法。它以负熵最大作为搜寻方向，可以实现顺序提取独立源；此外，还采用了定点迭代的优化算法，使收敛更加快速、稳健。

MT时间序列信号可视为随机变量在某个时刻的采样，因此可以将时间序列信号设为随机变量。假定这里有m个观测信号，第i个观测信号xi是由n个相互独立的未知源信号s1、s2、…、sn混合而成(xi、sj均为随机变量)，即

xi=ai，1s1+ai，2s2+…+ai，nsn

(1)

上式中ai,j(j=1,2,…,n)为常数系数。用矢量X表示观测信号变量(x1,x2,…,xm)T，矢量S表示源信号变量(s1,s2,…,sn)T(m≥n)，Am×n表示混合矩阵(ai,j)，则式(1)可表示为

X=AS

(2)

FastICA的目标就是寻求分离矩阵Wm×n，以使Y=WTX=WTAS具有最大的非高斯性。

设离散随机变量y的熵定义为

(3)

式中：ki指y的可能取值；P(·)指求概率。式中对数一般取2为底，单位为比特。其负熵定义为

Ng(y)=H(yGauss)-H(y)

(4)

式中yGauss是一个与y具有相同方差的高斯随机变量。负熵Ng(y)可以作为随机变量y非高斯性的测度。由于计算离散随机变量的熵需要知道随机变量的概率分布，其不易获取，因此常采用下式计算负熵

Ng(y)={E[g(y)]-E[g(yGauss)]}2

(5)

式中g(·)为特定非线性函数。

FastICA算法原理如下。

首先，通过对E[g(WTX)]进行优化获得y的负熵的最大近似值。根据Kuhn-Tucker条件，在E[(WTX)2]=‖W‖2=1的约束下，E[g(WTX)]的最优值在满足下式的点上获得

E[Xg(WTX)]+βW=0

(6)

然后，利用牛顿迭代法求解上述方程。记P=E[Xg(WTX)]+βW，可得P的雅可比矩阵为

JP=E[XXTg′(WTX)]-βI

(7)

式中g′表示g的一次导数。由于数据被白化，有E(XXT)=I,故

E[XXTg′(WTX)] ≈E(XXT)·E[g′(WTX)]

=E[g′(WTX)]I

(8)

因而JP变成了对角阵，容易求逆。如此便得到近似牛顿迭代公式

(9)

式中：W*是W进行一次迭代后的取值；β=E[WTXg(WTX)]。该式简化后就得到FastICA算法的迭代公式

(10)

MT观测的是电场和磁场的时间序列，把获得的时间序列看作一个随机变量在特定时候的采样，即将记录的电场或磁场时间序列视为随机变量x1，x2，…，xm,则可以利用上面所述的FastICA算法进行迭代处理，便可在获得分离矩阵W后获得具有最大的非高斯性的随机变量矢量

Y=WTX=WTAS

(11)

Y既包含原始信号也包含噪声，但因为可以将大地电磁资料中的原始电磁信号和噪声看作是相互独立的信号源，因此可以利用FastICA算法进行原始信号与噪声的分离和提取。

2.3 基于小波分析的独立分量分析

经过多次试验后发现，小波变换对高信噪比的含噪信号去噪效果较好，对低信噪比的含噪信号去噪效果较差，而属于盲源分离的独立分量分析算法对低信噪比的含噪信号的去噪效果却较好。

根据盲源分离理论[21]，在不知道源信号及混合矩阵的任何信息的情况下，只需假设源信号是相互统计独立的，就能使用独立分量分析将源信号从混合信号中分离出来。在大地电磁勘探的实际应用中，并不能通过有效方法证明源信号之间是严格独立的，更不能证明源信号是严格线性瞬时混合的，但从各种源信号的不同存在方式来看，它们之间至少是独立的，也就是说它们满足相互统计独立的要求。从另一方面来说，由于大地电磁有效信号与噪声信号的来源不一致，因而它们的自身属性也不相同，所以从属性角度看它们也是相互独立的。因此，可以把大地电磁资料中的有效信号与噪声信号视为相互统计独立的信号，这样就可以使用独立分量分析算法来对它们进行分离。

但独立分量分析也有其弊端，那就是它的适用条件要求采集到的信号的数目要大于或者等于源信号的数目[16-19]，这就意味着必须同时获得至少两组在同时同地采集到的数据才能将所得数据分离为有效信号和噪声信号两组数据。而且，一般而言，使用的采集信号的组数越多，提供的噪声特征的信息量越丰富，因此在许多情况下要求使用在同一时间、同一地点采集到的至少三组大地电磁信号来进行实验和研究。而在实际野外作业中，由于地质条件、气候条件、环境条件及人力物力条件的制约，要获得这样的三组及以上的数据是非常困难的，且成本非常高。受陈艳等[22]研究的启发，首先对信号进行小波分解和去噪处理，再应用独立分量分析技术进行盲源分离，以提取源信号和噪声。为了减小观测信号的信噪比对去噪效果的影响，本文引入动态的自适应因子∂，根据观测信号的幅值情况不同而取不同的∂值；另外，由于噪声一般包含在小波高频分量中，所以重构信号由不含噪声的低频分量(近似信号)获得。根据如下公式重构信号

Y=R+∂*T

(12)

式中：R为小波近似信号；T为ICA分离后提取的有效源信号；Y为去噪后信号。

基于小波分析的ICA去噪的算法流程如图2所示。

图2 基于小波分析的ICA去噪算法流程

3 仿真试验

为研究独立分量分析算法对信号处理的效果，首先进行仿真试验。分别采用周期信号和非周期信号进行实验研究。

以t表示时间变量，周期信号取y1=20sin150t+5cos100t，非周期信号取leleccum信号，即y2=leleccum(t)，这里的leleccum信号为matlab所提供的一种非周期非规则信号。

分别运用七种方法进行去噪： ①启发式软阈值去噪法(5层分解)； ②启发式硬阈值去噪法(5层分解)； ③SURE软阈值去噪法(5层分解)； ④SURE硬阈值去噪法(5层分解)； ⑤自适应默认阈值去噪法(5层分解)； ⑥重构信号Y=R+∑T′去噪法(T′为ICA分离后未经处理的原始结果信号)(5层分解)； ⑦由式(12)所示的去噪法(5层分解)。得到的实验结果见图3。限于篇幅，这里只显示加入噪声(SNR=10)时y1的实验结果。

图3 原始信号和去噪结果 (a)原始信号和含噪信号； (b)含噪信号5层小波分解后的结果； (c)ICA处理后的结果； (d)七种方法分别处理后的结果

为了比较去噪效果，这里采用Pearson(皮尔逊)相关系数ρY,O度量去噪后信号Y与原始信号O(无噪信号)之间的相关度

(13)

显然Y、O包含的元素个数相等。

表1和表2是分别运用上述七种去噪方法对信号y1和y2进行去噪处理后的结果。可以看出，无论哪一种信号、无论信噪比大小，本文所采用的去噪方法的去噪效果都比较稳定。所采用的恢复信号既包含低频小波系数，又含有独立分量分析分离出的源信号，并且通过自适应因子动态改变独立分量分析算法的分离效果，使观测信号的信噪比对去噪结果的影响减小。

由表2可以发现，对信号y2而言，方法⑥明显比方法⑦更好，这是因为信号y2含有的高频成分较多，把ICA分离得到的原始信号进行直接叠加，减少了高频成分的损失。也就是对变化较快的高频信号而言，方法⑥会比方法⑦效果更好。但就MT信号来说，大多数情况下需要获取的是低频率信号、而不是高频信号，因此应用于大地电磁测深数据时，方法⑦整体性能仍然优于方法⑥。

图4对比了方法①与方法⑦的稳定性，可以发现方法⑦的稳定性优于方法①。其他方法(方法②～方法⑥)与方法⑦的稳定性比较结果与之类似，这里不再赘述。

表1 信号y1七种方法去噪处理后的相关系数

表2 信号y2七种方法去噪处理后的相关系数

图4 方法①与方法⑦的稳定性对比图 (a)数据y1计算结果； (b)数据y2计算结果

4 实际观测资料去噪处理

选取含有噪声的两个测点(分别记为A和B)进行研究，其中测点A采自内蒙古境内，测点B采自河南省境内，因环境影响，测点B含噪声较多。对这两个测点数据运用本文所述的方法进行去噪处理，得到图5和图6所示的一系列结果(限于篇幅，仅展示1000个采样点的结果)。

图7和图8分别是A、B两点去噪前、后的视电阻率曲线和相位曲线。可以看出，运用本文方法进行去噪处理后，无论视电阻率曲线还是相位曲线，都

图5 测点A去噪前(左)、后(右)的时间序列对比 (a)Ex； (b)Ey; (c)Hx; (d)Hy; (e)Hz

图6 测点B去噪前(左)、后(右)的时间序列对比 (a)Ex； (b)Ey; (c)Hx; (d)Hy; (e)Hz

图7 测点A(a)和测点B(b)去噪前(左)、后(右)视电阻率曲线对比

图8 测点A(a)和测点B(b)去噪前(左)、后(右)相位曲线对比

变得更加光滑、流畅，即去噪效果都比较好。虽然仍有极少数频点去噪效果并不理想，但这些频点在数据预处理时可以作为“飞点”予以删除，而并不影响整条曲线的处理结果。

5 结论

本文结合小波分析与独立分量分析的优势，提出基于小波分析的独立分量分析去噪方法并对大地电磁数据进行去噪。理论上证明了将独立分量分析应用于大地电磁数据去噪中的可行性；模拟信号仿真实验表明该方法的去噪稳定性优于传统小波阈值去噪方法；实际大地电磁观测资料去噪结果说明，此方法能比较有效地去除噪声，特别是在非极低频区域，去噪效果尤其明显。经该方法处理后的视电阻率曲线和相位曲线变得光滑、流畅，为后续数据处理提供了高质量的数据。

[1] 谢成良.大地电磁测深资料综合处理软件系统研究[学位论文].北京:中国地质大学(北京),2013.

[2] 程正璞,胡祥云,李烨等.民丰凹陷大地电磁探测研究.石油地球物理勘探,2016,51(2):391-403. Cheng Zhengpu,Hu Xiangyun,Li Ye et al.A magnetotelluric survey in Minfeng Sag.OGP,2016,51(2):391-403.

[3] 严良俊,胡文宝,陈清礼等.远参考MT方法及其在南方强干扰地区的应用.石油天然气学报,1998，18(4):34-38. Yan Liangjun,Hu Wenbao,Chen Qingli et al.Application of remote reference MT to noisy area.Journal of Oil and Gas,1998，18(4):34-38

[4] 蔡剑华,龚玉蓉,王先春.基于Hilbert-Huang变换的大地电磁测深数据处理.石油地球物理勘探,2009,44(5):617-621,629. Cai Jianhua,Gong Yurong,Wang Xianchun.Magnetotelluric data processing based on Hilbert-Huang transform in oil and gas exploration.OGP,2009,44(5):617-621，629.

[5] 张全胜,王家映.大地电磁测深资料的去噪方法.石油地球物理勘探,2004，49(增刊1):17-23. Zhang Quansheng,Wang Jiaying.The method of de-noising for magnetotelluric data.OGP,2004，49(S1):17-23.

[6] 蔡剑华,王先春.EMD在大地电磁信号分析中的问题及解决方法.石油地球物理勘探,2016,51(1):204-210. Cai Jianhua,Wang Xianchun.Problems and solutions of EMD in magnetotelluric signal analysis.OGP,2016,51(1):204-210.

[7] 汤井田,李晋,肖晓等.数学形态滤波与大地电磁噪声压制.地球物理学报,2012,55(5):1784-1793. Tang Jingtian,Li Jin,Xiao Xiao et al.Mathematical morphology filtering and noise suppression of magnetotelluric sounding data.Chinese Journal of Geophy-sics,2012,55(5):1784-1793.

[8] Trad D O,Travassos J M.Wavelet filtering of magnetotelluric data.Geophysics,2000,65(2):482-491.

[9] Manoj C,Nagarajan N.The application of artificial neural networks to magnetotelluric time-series analysis.Geophysical Journal of the Royal Astronomical Society,2010,153(2):409-423.

[10] Weckmann U,Magunia A,Ritter O.Effective noise se-paration for magnetotelluric single site data processing using a frequency domain selection scheme.Geophysical Journal International,2005,161(3):635-652.

[11] 何兰芳,王绪本,何展翔等.MT时间序列的小波去噪分析.地震地质,2001,23(2):222-226. He Lanfang,Wang Xuben,He Zhanxiang et al.Wavelet-based denoising of MT time series.Seismology & Geology,2001,23(2):222-226.

[12] 严家斌,刘贵忠.基于小波变换的脉冲类电磁噪声处理.煤田地质与勘探,2007,35(5):61-65. Yan Jiabin,Liu Guizhong.Like-impulse electromagnetic noise processing based wavelet transform.Coal Geology & Exploration,2007,35(5):61-65.

[13] 蔡剑华,肖晓.基于小波自适应阈值去噪的MT信号处理方法.地球物理学进展,2015,30(6):2433-2439. Cai Jianhua,Xiao Xiao.Method of processing magnetotelluric signal based on the adaptive threshold wavelet.Progress in Geophysics,2015,30(6):2433-2439.

[14] 王书明,王家映.大地电磁信号统计特征分析.地震学报,2004,26(6):669-674. Wang Shuming,Wang Jiaying.Analysis on statistic characteristics of magnetotelluric signal.Acta Seismologica Sinica,2004,26(6):669-674.

[15] 王书明,王家映.关于大地电磁信号非最小相位性的讨论.地球物理学进展,2004,19(2):216-221. Wang Shuming,Wang Jiaying.Discussion on the non-minimum phase of magnetotelluric signals.Progress in Geophysics,2004,19(2):216-221.

[16] Rinen A,Oja E.A fast fixed-point algorithm for independent component analysis.MIT Press,USA，1997,1483-1492.

[17] Hyvarinen A,Oja E.Independent component analysis:algorithms and applications.Neural Networks,2000,13(4):411-430.

[18] 杨福生,洪波.独立分量分析的原理与应用.北京:清华大学出版社,2006,97-101.

[19] 海韦里恩,卡尔胡恩,奥亚.独立成分分析.北京:电子工业出版社,2014,138-152.

[20] 盛骤,谢式千.概率论与数理统计及其应用.北京:高等教育出版社,2004,87.

[21] 余先川,胡丹.盲源分离理论与应用.北京:科学出版社,2011,35-42,51-104.

[22] 陈艳,何英,朱小会.基于小波变换的独立分量分析及其在图像分离中的应用.现代电子技术,2007,30(24):131-134. Chen Yan,He Ying,Zhu Xiaohui.ICA based wavelet transform and its application to image separation.Modern Electronics Technique,2007,30(24):131-134.

[23] 胡昌华.基于MATLAB的系统分析与设计—小波分析.陕西西安:西安电子科技大学出版社,1999,223-225.

[24] Donoho D L,Johnstone I M.Ideal spatial adaptation by wavelet shrinkage.Biometrika,1994,81(3):425-455.