大地电磁测深(magnetotelluric sounding, MT)是研究地球深部构造的重要方法之一, 该方法利用天然电磁场作为场源因而易受到各种干扰的影响, 其中一种重要的干扰类型为工频干扰。工频干扰是大地电磁信号采集过程中最为普遍的干扰之一, 近年大地电磁信号中的工频干扰噪声愈加严重, 极大影响了大地电磁勘探的效果。工频干扰噪声的信号能量强, 一般的工频干扰噪声能量为正常大地电磁信号的几倍, 有的工频干扰噪声能量甚至比正常大地电磁信号高几个数量级。在高压输电线附近, 工频干扰噪声信号能量尤其巨大, 常常使得大地电磁有效信号完全淹没在工频干扰噪声中, 严重影响后续的数据处理和反演解释^[1-3]。随着国内外对大地电磁信号噪声压制的重视和持续研究, 一些学者针对广泛存在的工频干扰噪声问题提出了解决方案。蔡剑华等^[4]提出先对包含工频干扰噪声的MT信号进行傅里叶变换, 得到其实部和虚部后, 再用小波阈值去噪法分别对实部及虚部进行去噪处理, 最后将去噪后的实部和虚部联合起来进行傅里叶逆变换得到去噪后的信号。TANG等^[5]提出先对采集的大地电磁信号进行傅里叶变换, 然后设计与干扰信号相匹配而对有用信号不敏感的冗余字典原子, 再结合改进的正交匹配追踪(improved orthogonal matching pursuit, IOMP)算法分离出频域信号中的工频干扰成分, 最后将处理后的频域信号进行傅里叶逆变换得到去噪后的信号。上述两种方法压制大地电磁信号工频干扰噪声效果较好, 但均为基于傅里叶变换及其反(逆)变换的方法, 故往往受限于傅里叶变换的优势和劣势。CAO等^[6]和曹小玲等^[7]提出基于盲源分离的大地电磁工频干扰噪声压制方法, 将受到工频干扰噪声的大地电磁信号进行离散小波变换(discrete wavelet transform, DWT)和集合经验模态分解(ensemble empirical mode decomposition, EEMD)处理, 而后再进行盲源分离以压制噪声, 该方法采用了基于DWT的处理模型并引入了自适应权重因子, 不仅有效降低了盲源分离算法恢复信号幅值的不确定性, 而且在处理含有工频干扰噪声且信噪比较低的MT数据时仍然效果良好。该方法采用小波变换方法构造多通道信号, 因小波变换基函数的选择和分解层次的选择过于依赖于人工经验及试验分析, 故方法的自适应性不强。我们在研究中一直试图避免采用小波变换方法, 探索不过分依赖人工经验及试验分析的信号分解方法和工频干扰噪声压制方法。

经验模态分解^[8-9](empirical mode decomposition, EMD)显著特点为: 根据信号特性通过迭代的方式自适应地获取基函数和分解层次, 即无需事先人为给定基函数和分解层次, 无需依赖人工经验及试验分析, 因其具有良好的自适应性故被广泛应用于噪声压制。盲源分离(blind source separation, BSS)是近年来发展起来的一种现代信号处理技术, 已在通信信号处理、机械故障识别、语音信号去噪、地震信号去噪等领域获得广泛应用, 前期我们将其应用于大地电磁信号去噪, 取得了较好的效果^[10]。为了结合经验模态分解和盲源分离的优势并弥补二者的不足, 本文提出一种基于单通道盲源分离的大地电磁工频干扰噪声压制方法, 主要利用工频干扰噪声的特点、EMD及盲源分离的优良特性, 首先将含有工频干扰噪声的大地电磁信号进行经验模态分解(EMD)和主成分分析(principal component analysis, PCA), 然后利用衰减率确定有效的固有模态函数(intrinsic mode function, IMF)的个数(即源信号个数, 也即有效主成分的个数), 最后利用盲源分离中的独立分量分析(independent component analysis, ICA)方法对有效主成分进行处理, 以压制工频干扰噪声。本文方法根据观测信号确定源信号的个数, 解决了在欠定盲源分离情况下源信号个数无法确定的问题。

1 大地电磁工频干扰噪声特征

大地电磁天然场源的特点^[11-12]决定了大地电磁信号非常容易受到噪声的干扰, 且这些噪声形态各异、种类繁多、强弱不一。根据噪声产生的机理划分, 大地电磁测深信号包含的噪声包括工频干扰噪声、仪器噪声、人文噪声、地质噪声等。工频干扰噪声一般指电力传输系统、人工电磁场、通信过程及有线广播等产生的噪声。这些电磁噪声呈现非平面波的形态, 且与观测点距离较近, 因此不符合大地电磁测深对场源的要求。

因工频干扰噪声很大一部分是人为导致的, 故也有学者认为应将其归为人文噪声。如文献[1]认为大地电磁测深中的噪声包括场源噪声、地质噪声和人文噪声3个类别, 其中人文噪声主要是指“人类活动所产生的干扰性电磁波以及人类的活动产生的其它电磁噪声”^[1], 这里的“人文噪声”包括了工频干扰噪声, 即认为工频干扰噪声是人文噪声中的一部分。

工频干扰噪声主要特征为: ①等频率(50 Hz)且等振幅, 通常导致原始的大地电磁有效信号几乎完全被淹没而不可见; ②形态较为规则, 一般呈现出类似于正弦曲线的形状; ③幅值通常很大, 因为与原始有效信号幅值对比强烈, 故常常使得原始有效信号无法被识别, 因此对电磁勘探等实际生产的破坏性非常大; ④谐波干扰(一般为奇次谐波干扰)严重, 影响后期的数据处理和反演工作。本文从实际大地电磁信号中挑选出了一段包含工频干扰噪声的实测大地电磁信号(图 1), 该测点地址位于31:08.487(N), 114:06.865(E)。一般而言, 工频干扰噪声可能同时出现在所有电道和磁道中, 如图 1所示E_x、E_y、H_x、H_y、H_z各道信号中都出现了工频干扰噪声, 原始有效信号几乎被完全覆盖。

2 方法原理 2.1 盲源分离基本原理

盲源分离^[13]的流程如图 2所示, 其中, s(t)=[s₁(t), s₂(t), …, s_n(t)]^T是n维未知源信号向量, A是未知混合系统, x(t)=[x₁(t), x₂(t), …, x_m(t)]^T是m维观测信号向量, 它们是源信号向量受噪声向量n(t)=[n₁(t), n₂(t), …, n_m(t)]^T干扰而形成的。x(t)被分离系统W处理后得到输出信号y(t), 它是源信号向量s(t)的估计。

盲源分离的目的是在源信号s(t)和混合系统A及噪声n(t)均未知的情况下, 利用观测信号向量x(t)调整分离系统W(寻求分离矩阵W), 使得输出y(t)是源信号s(t)的估计, 即:

$ \mathit{\boldsymbol{y}}\left( t \right) = \mathit{\boldsymbol{W}}\left( {\mathit{\boldsymbol{x}}\left( t \right)} \right) \cong \mathit{\boldsymbol{s}}\left( t \right) $

(1)

根据王书明等^[14-15]对不同地区实测大地电磁信号的性质特征分析可知, 大地电磁信号符合盲源分离对源信号的统计性质要求。

独立分量分析(ICA)方法是近年发展的一种有效盲源分离技术, 其应用范围极为广泛^[16-17]。从线性变换和线性空间角度来看, 源信号为相互独立的非高斯信号, 可以看作线性空间的基信号, 而观测信号为源信号的线性组合, ICA方法是在源信号和线性变换均不可知的情况下, 从观测的混合信号中估计数据空间的基本结构或者信源信号^[7]。ICA方法的算法流程^[7]如图 3所示, 在信源信号s(k)中各分量相互独立的假设下, 由观察信号x(k)通过解混系统B将信源信号分离, 使输出y(k)逼近s(k)。ICA方法中应用较多的是FastICA算法, 本文也采用此算法。

2.2 EMD方法基本原理

EMD方法将相对复杂的原始信号分解成有限个相对简单且具有不同特征尺度的数据序列, 即固有模态函数(IMF)分量, 反映了原始信号的本质和真实信息, 因此我们对其进行操作和处理就等同于对原始信号进行操作和处理。EMD方法的基本步骤如下。

1) 假设原始信号为x(t), 首先我们找出x(t)中所有的局部极值点(包括极大值点和极小值点), 用三次样条函数连接所有的局部极大值点构成上包络线, 用三次样条函数连接所有的局部极小值点构成下包络线。

2) 求出上包络线与下包络线的均值(设为m₁), 我们设信号x(t)与m₁的差为h₁, h₁= x(t)－m₁, 如果h₁满足IMF的条件, 那么h₁就是原始信号的第一个IMF分量。

3) 如果h₁不是IMF分量, 我们就将h₁看成原始信号进行处理, 再根据上述步骤, 得到h₁₁=h₁－m₁₁, 其中, m₁₁是h₁信号的上、下包络线均值。按照上面的操作反复处理k次得到h_1k, 如果h_1k满足IMF的条件, h_1k被称为IMF, 且有h_1k=h_1(k－1)－m_1k, c₁=h_1k为原始信号的第一个IMF分量, 代表x(t)最高频率的分量。

4) 从x(t)中分离c₁, 得到r₁= x(t)－ c₁, 将r₁视为原始信号重复上述过程, 可以得到信号x(t)的第二个IMF分量c₂。重复上述过程n次, 得到n个IMF分量, 依次为c₁, c₂, …, c_n, 且r₂= r₁－ c₂, …, r_n= r_n－1－ c_n。当r_n为一个极小的常量或单调函数时, 就无法提取IMF, 此时我们终止分解过程, 可以得到:

$ \mathit{\boldsymbol{x}}\left( t \right) = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{c}}_i} + {\mathit{\boldsymbol{r}}_n}} $

(2)

其中, r_n称为残差, 它反映了信号x(t)的集中趋势。IMF分量c₁, c₂, …, c_n包含了信号的特征尺度信息, 也包含了信号的特征频率信息, 因此利用c₁, c₂, …, c_n就可以了解原始信号完整的特征尺度信息和特征频率信息。

本文没有选择EEMD而选择EMD处理的原因如下: 一是因为EEMD是在信号中添加白噪声后再进行EMD处理, 添加白噪声可能对原有的工频干扰噪声产生影响; 二是因为大地电磁信号普遍数据量庞大, 采用EEMD处理会极大地增加计算量。

2.3 基于单通道盲源分离的大地电磁工频干扰噪声压制方法

大地电磁信号通常是由多道信号数据构成, 由于每道信号数据包含的噪声形式和噪声数量不尽相同, 因此我们需要对每道信号数据分别进行去噪处理, 即需要进行单通道大地电磁工频干扰噪声压制。盲源分离一般要求观测信号的数目不少于源信号的数目, 而单通道大地电磁工频干扰噪声压制问题属于欠定盲源分离问题, 即观测信号的数目少于源信号的数目, 所以必须构造其它观测信号, 增加观测信号的个数, 使得观测信号的数目不少于源信号的数目, 才能应用盲源分离方法。本文提出的基于经验模态分解和盲源分离的大地电磁工频干扰噪声压制方法, 可以实现单通道大地电磁工频干扰噪声压制。该方法的流程如图 4所示, 首先采用EMD将观测信号分解为若干个IMF分量(这里假设为M个), 即实现观测信号的升维, 以满足盲源分离对源信号数目的要求; 然后将IMF分量组成矩阵, 应用奇异值分解方法求取矩阵特征值, 根据特征值比求衰减率, 根据衰减率确定源信号的数目和有效主成分的个数; 再对获得的IMF分量进行PCA分析, 确定并选择起主要作用的主成分分量; 为进一步压制噪声和减少主成分信号之间的相关噪声, 最后对选取的若干个有效主成分分量进行ICA处理, 获得压制工频干扰噪声之后的信号。

该方法主要优势在于EMD模型、PCA处理及盲源分离方法的应用, 受文献[18]的启发, 我们利用衰减率确定源信号的数目和待处理的主成分的个数, 解决了欠定盲源分离中源信号个数无法确定的问题, 实现了单通道观测信号的源信号的分离和提取。

PCA处理在盲源分离处理领域是一种极其有效的处理手段, 我们对EMD分解得到的IMF分量进行PCA处理^[19], 可以有效降低信号数据的维数, 通过寻找信号中能量最大的分量得到信号的主要特征量, 即起到“化繁为简”的作用。另外, PCA处理作为一种统计分析方法, 还能实现对信号的去相关处理。因为它能够使得变换后产生的新分量正交或不相关, 故变换后的分量能量更集中、性质更稳定^[20]。PCA处理的优势在于利用某种变换将数据原有的大量特征变换为几个主要特征, 而这些特征包含了原始数据的最主要的特征信息。因此, 当我们仅仅提取前面的部分特征时, 可以既减小特征的个数又保留原有信号最主要的特征信息。同时, 由于PCA处理后的信号不相关, 所以噪声子空间和源信号子空间相互分离, 进而可以达到初步压制部分噪声的目的, 即初步压制部分工频干扰噪声。

PCA处理虽然能够去除信号的相关性, 但在高阶累积量的定义下信号可能仍然具有相关性。为使所得信号具有的统计独立性最大, 我们又进行了ICA处理, 选取有用分离向量与混合矩阵相乘重构得到分离后的独立分量信号, 进而得到消噪信号。

3 仿真模拟实验

如图 5a所示, 利用蒙特卡洛方法对仿真信号中的原始信号随机构造3组信号(幅值为400)。因为工频干扰噪声一般为50 Hz的周期信号及其谐波, 故构造工频干扰噪声信号(噪声信号)可表示为S=800 sin(2π×50t)。现实情况下强烈的工频干扰噪声主要来自电力传输系统, 其频率一般为50 Hz, 从方便分析的角度考虑, 本文设工频干扰噪声为50 Hz。因工频干扰噪声信号幅值一般是有效电磁信号幅值的数倍, 所以这里幅值应设较大值, 即为原始信号的数倍(图 5b)。图 5c为工频干扰噪声加入原始信号后得到的含噪信号。

对图 5c所示的含噪信号进行EMD分解, 结果如图 6a所示, 对该分解结果进行主成分分析(PCA), 结果如图 6b所示, 受篇幅限制, 只展示图 5中左侧信号的处理结果。因为工频干扰噪声的幅值和能量均较大, 因此它可能被分解到不同阶的IMF分量中(出现EMD分解的混频现象), 对这些IMF分量进行主成分分析后发现, 工频干扰噪声主要包含在前几个主要的主成分之中。因此如果需要提取噪声, 我们只需要对前几个主成分进行处理。正如我们在获得观测信号后无法确定其源信号数目一样, 确定有效主成分的个数非常困难。根据本文方法的思路和独立分量分析方法的前提条件, 如果我们确定了待处理的有效主成分个数, 也就确定了盲源分离问题中的源信号数目。

我们根据参考文献[18]利用衰减率确定源信号的数目, 进而确定待处理的有效主成分个数。计算由IMF分量所构成的矩阵的特征值及衰减率, 步骤如下: 假设IMF分量所构成的矩阵为X =[I_IMF1, I_IMF2, …, I_IMFN], 对其进行奇异值分解, 得到的特征值为λ₁, λ₂, …, λ_M, 衰减率为v_i=λ_i/λ_i+1(i=1, 2, …, M－1), 待处理的有效主成分的个数N计算公式如下:

$ 令\left( {{v_{max}}, k} \right) = {\rm{max}}\left\{ {{v_1}, {v_2}, \cdots , {v_{M - 1}}} \right\}, 则 N = k + 1 $

(3)

式中: v_max为衰减率最大值; k为衰减率最大值的序号。

为了验证此方法有效性, 将图 5中3组含噪信号分别按上述方法进行处理, 得到的衰减率如表 1所示。

此处k=1, 即v₁值最大, 故得到的有效主成分个数N为2, 与理论情况一致。如果观测信号主要为噪声, 则有用信号能量很弱, 通常这种情况下压制噪声几乎不可能, 因此要求工频干扰噪声幅值不得超过原始有用信号的10倍。

本文进行了多次实验, 如将50, 150, 250 Hz的正弦波分别添加或一起添加到原始信号中, 实验结果均证明利用上述方法可以确定源信号的数目。因为原始信号是随机构造的且上述方法相关理论知识较为成熟, 故采用该方法确定的源信号数目有效。

得到待处理的有效主成分个数N后, 我们可以只选取PCA处理后的前N个主成分进行ICA处理, 前述处理得到N=2, 故本文选取2个主成分(记为p1, p2)进行ICA处理。图 7a、图 7b和图 7c分别为原始信号、工频干扰信号和加噪后信号; 图 7d为对加噪后的信号进行ICA处理后的结果。可以发现, 在原始信号被干扰噪声完全淹没的情况下, 原始信号波形已被较好地提取出来(如图 7d中C₁所示), 但C₁的幅值与原始信号存在差距。根据观测信号的幅值范围确定权重因子^[10], 并将权重因子作用于C₁, 最终得到去噪后的信号。

图 8a、图 8b和图 8c分别为原始信号、加噪后信号和去噪后的信号。比较去噪后的信号与原始信号可以看出, 提取出幅值较大的工频干扰噪声后, 信号波形和幅值均得到了较好的恢复。

因为EMD算法的卓越性能和PCA分解能量的优化操作, 所以对EMD之后的IMF直接进行PCA处理, 并提取第一个主成分作为消噪后的信号, 也能达到一定的工频干扰压制效果。为了方便对比分析, 本文将这种方法称为PCA方法。表 2为噪声信号与原始信号的幅值之比R分别取5种不同值的情况下, 3种方法(小波方法、PCA方法、本文方法)的性能。为了比较提取的恢复信号与原始信号的相似程度, 比较去噪效果, 我们采用相关系数来度量去噪后信号与原始信号(无噪信号)的相似程度, 定义如下:

$ {\rho _{X, Y}} = \frac{{\sum\limits_{i = 1}^m {\left( {{X_i} - \bar X} \right)\left( {{Y_i} - \bar Y} \right)} }}{{\sqrt {\sum\limits_{i = 1}^m {{{\left( {{X_i} - \bar X} \right)}^2}} } \sqrt {\sum\limits_{i = 1}^m {{{\left( {{Y_i} - \bar Y} \right)}^2}} } }} $

(4)

式中: X_i, Y_i分别为信号X和Y中的元素; m为信号X和Y的元素数量; X, Y分别为信号X和Y的平均值。表 2中第2行为含噪信号与原始信号的相关系数, 第3、4、5行分别为采用3种方法得到的消噪信号与原始信号的相关系数。

从表 2中可以发现, 随着工频干扰噪声幅值的不断增加, 采用本文方法提取的去噪信号与原始信号相关系数始终维持在一个较高的水平(大于0.84), 而采用小波方法得到的结果始终维持在一个较低的水平(小于0.18), 这说明本文方法压制信号中的工频干扰噪声明显优于小波分析方法。PCA方法在R值较小时, 性能较优越, 但当噪声信号与原始信号的幅值之比大于1后, 其性能迅速下降, 即PCA方法稳定性欠佳, 而大地电磁中工频干扰噪声幅值一般都较大, 甚至是原始有效信号的数倍, 所以PCA方法的实用性不高。

4 实际大地电磁数据处理

选取西部某地采集的大地电磁测点信号作为研究对象, 因该地区电力线信号、通讯信号等均非常微弱, 故可以认为采集到的大地电磁测点信号是不包含工频干扰噪声的原始信号。此测点信号的视电阻率曲线和相位曲线如图 9所示, 曲线形态平滑、流畅, 也从侧面印证此测点为未被噪声污染的理想测点。

此测点原始信号的电场和磁场时间序列如图 10a所示(考虑到篇幅限制这里仅展示1 000个采样点), 对该测点的5道信号分别添加50 Hz工频干扰噪声(幅值为10 000), 结果如图 10b所示, 可以看出, 添加工频干扰噪声后, 原始信号被完全淹没, 几乎完全不可见。此时若直接进行数据处理, 视电阻率曲线和阻抗相位曲线均会在50 Hz附近出现畸变现象, 进而影响后续反演解释工作, 因此必须压制工频干扰噪声。

采用本文方法进行工频干扰噪声压制, 得到了整个测点的工频干扰噪声压制结果, 本文限于篇幅, 仅选取其中1 000个采样点展示处理结果(图 11)。图 11的各个子图中, 从上到下依次为原始信号、添加了工频干扰噪声的信号、提取的工频干扰噪声、压制工频干扰噪声后的信号(即消噪信号)。可以看出, 原始信号添加工频干扰噪声后被工频干扰噪声完全淹没, 几乎看不出任何原始信号的信息。但经过本文方法处理后, 原始信号的概貌基本得到了恢复。

表 3列出了该测点(整个时间序列)5道信号的提取噪声与添加噪声的相关系数、消噪信号与原始信号的相关系数, 可以看出, 本文方法对大地电磁中工频干扰噪声的提取效率非常高, 5道相关系数均达到了0.96, 接近于1。从表 3还可以看出: 消噪信号与原始信号的相关系数虽然也较高, 但无法达到0.93。这是因为大地电磁信号属于不稳定的随机信号, 添加大幅值的工频干扰噪声后, 对原始信号能量及幅值均有较大影响, 故目前还难以完全恢复出原始信号。

添加工频干扰噪声(幅值以万为计数单位)后以及压制工频干扰噪声后得到的视电阻率曲线和相位曲线分别如图 12所示(因MTL在添加噪声前后几乎没有变化, 故这里只展示MTH在去噪后的结果)。图 12a为含工频干扰噪声大地电磁信号的视电阻率曲线和相位曲线, 可以发现, 视电阻率曲线在50 Hz(横坐标值约为1.7)附近存在畸变, 在5 Hz(横坐标值约为0.7)附近存在较大畸变, 阻抗相位曲线在50 Hz和5 Hz附近均存在较大畸变。图 12b为本文方法处理结果, 本文方法使视电阻率曲线和相位曲线在50 Hz附近的偏移获得了较大的改善。因为添加的工频干扰噪声幅值较大, 采用本文方法能将原始测点信号的视电阻率曲线和相位曲线恢复至现在的程度, 因此说明本文方法在压制工频干扰噪声方面具有较好的效果。不足的是, 本文方法无法同时使5 Hz附近曲线平滑, 如何对5 Hz附近信号进行有效处理, 是今后的研究方向。

5 结论

本文提出基于经验模态分解和盲源分离的大地电磁工频干扰噪声压制方法, 摆脱了小波分析等经典方法对小波基函数和分解层次的限定和人工经验的制约, 充分利用了经验模态分解和盲源分离的优良特性。该方法在利用衰减率确定出源信号个数的同时, 利用PCA处理EMD分解后的IMF, 而后采用ICA处理PCA提取的若干有效主成分, 最终得到了压制工频干扰噪声之后的去噪信号。该方法在工频干扰噪声的幅值为原始信号幅值的数倍的情况下能在一定程度上压制工频干扰噪声, 鲁棒性较强。实验研究发现, 工频干扰噪声信号与原始信号的幅值比为4~7, 本文算法的去噪性能最佳, 未来将着眼于拓展算法的应用范围, 并结合人工智能技术, 将此方法进一步推向智能化。

致谢: 感谢长江大学物探重点实验室为本研究提供的野外实际大地电磁观测数据。