大地电磁测深(MT)是一种被动源电磁探测技术, 通过在地表测量正交电、磁场分量的扰动值, 得到地下介质的电性结构信息^[1]。MT探测范围较宽, 应用领域较广^[2], 但观测的天然电磁场信号较弱、频带较宽, 极易受到各种噪声的干扰, 给后续的数据处理、反演和地质解释结果带来误差。为此, 地球物理研究人员提出了一系列去除大地电磁噪声的技术^[3-7]。

传统的电磁数据去噪技术基于傅里叶变换——功率谱估计, 因其无法表达信号的时域局部性质, 人们先后应用小波分析方法消除大地电磁的噪声^[8-13]。2004年, 王书明等通过研究不同地区实测MT信号的特征, 指出MT信号一般具有非高斯性、非线性和非最小相位性质^[14-15], 需要利用新的信号处理方法来解决大地电磁的噪声问题。近年来, 独立分量分析(ICA)技术^[16-19]被广泛应用于地震勘探、语音信号处理、图像处理等领域的信号分析和信号处理, 在这些领域的信号去噪中表现尤其突出。该方法不需要源信号的相关信息, 也不需要观测信号中源信号的混合方式, 通过对观测信号的分析和处理推测甚至恢复出源信号。这个特点对于地球物理勘探来说意义重大, 因为源信号的相关信息(先验知识)不太容易获得, 即使获得也成本高昂。独立分量分析要求待处理信号为非依赖性、非高斯性和非白性, 而大地电磁观测信号恰好符合这些要求。另外, 独立分量分析一般要求处理的对象是由一组在统计意义上相互独立的源信号构成的混合信号, 而大地电磁观测信号中噪声干扰和信号分别由不同的信源产生, 彼此相互独立, 符合独立分量分析对处理对象的这一要求。因此, 独立分量分析技术可以用来去除大地电磁观测信号中的噪声干扰。该方法是在统计独立意义下对混合信号进行去噪, 显然比仅仅依赖于不相关性度量的其它去噪方法效果更好。本文根据大地电磁测深数据量大的特点以及独立分量分析中M-FastICA算法的优良性能, 结合小波分析和盲源分离的相关理论, 研究并提出了一种改进的独立分量分析算法, 利用实际大地电磁资料验证了算法的有效性。

独立分量分析是20世纪90年代发展起来的一种新的统计方法。它假定第i个观测信号x_i是由n个相互独立的未知源信号s_j混合而成:

假定每个观测信号x_i和未知源信号s_j都是随机变量, 矢量x表示观测信号变量(x₁, x₂, …, x_m)^T, 矢量s表示源信号变量(s₁, s₂, …, s_n)^T(m≥n), 矩阵A_m×n表示混合矩阵(a_ij)_m×n, 则公式(1)可用矢量-矩阵形式表示为:

公式(2)也称为独立分量分析模型^[17]。

独立分量分析的目标是寻求分离矩阵W_m×n, 得到独立分量s的估计y=W^Tx=W^TAs, 使y具有最大的非高斯性。FastICA算法^[16-19]是独立分量分析的一种快速算法, 该算法采用牛顿迭代法对观测变量x的大量采样点进行批处理, 每次从观测信号中分离出一个独立分量, 采用负熵(Negentropy)^[17]作为非高斯性度量的标准, 应用固定点迭代理论寻找W^Tx的非高斯性最大值。

记$\mathit{\boldsymbol{W}} = [{\mathit{\boldsymbol{w}}_1}\;{\mathit{\boldsymbol{w}}_2}\; \ldots \;{\mathit{\boldsymbol{w}}_n}]$, 则${\mathit{\boldsymbol{W}}^{\rm{T}}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{w}}_1^{\rm{T}}}\\ {\mathit{\boldsymbol{w}}_2^{\rm{T}}}\\ \vdots \\ {\mathit{\boldsymbol{w}}_n^{\rm{T}}} \end{array}} \right]$, 故$\mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{w}}_1^Tx}\\ {\mathit{\boldsymbol{w}}_2^Tx}\\ \vdots \\ {\mathit{\boldsymbol{w}}_n^Tx} \end{array}} \right]$。设y_i是y的第i个分量, 为了得到y, 需要确定每个y_i=w_i^Tx中的w_i。由文献[17]知, FastICA算法中求w_i的迭代公式为:

式中:g(·)为某个特定非二次函数的导数, E{·}为求数学期望运算, β=E{w₀^Txg(w₀^Tx)}。

实际FastICA算法运行过程中, 由于采用牛顿迭代法, 收敛速度至多是三阶, 可能会出现收敛效果不理想或者收敛速度慢的情况, 特别是对于海量数据的处理更会如此。因此, 本文采用独立分量分析的M-FastICA算法^[20]来对数据量庞大的大地电磁观测信号进行处理。迭代过程为:

式中, F(w)=E{xg(w^Tx)}-βw, 其雅可比矩阵函数为J(w)=E{xx^Tg′(w^Tx)}－βI。

可以证明, (4)式的收敛阶数为m+1, 且每m次迭代只需计算一次雅可比矩阵J(w), 因此提高了独立分量分析的速度。一般取m=2时的(4)式为M-FastICA算法的迭代公式, 即:

笔者在利用小波分析方法去除大地电磁噪声的过程中, 参照文献[21]并进行多次试验后发现, 小波变换方法对于较高信噪比的含噪声信号去噪效果较好, 但对较低信噪比的含噪声信号去噪效果较差。大地电磁信号本身较弱, 受到的噪声干扰相对来说可能比较强, 若单独采用小波分析方法去噪可能得不到理想的效果。而属于盲源分离的独立分量分析方法对于低信噪比的含噪声信号去噪效果较好, 故本文将M-FastICA算法与小波分析方法结合起来进行去噪处理。这里进行小波分析的目的有两个:一是在去噪前期消除部分噪声, 二是解决将单个时间序列视为观测信号时观测信号数目不足的问题。

根据盲源分离理论^[22], 如果所获得的信号可以近似地认为是由各个独立的源信号经过混合之后所得的结果, 那么就可以运用独立分量分析提取出源信号。就大地电磁测深而言, 虽然我们无法证明源信号和噪声信号之间是严格独立的, 但是因为这两种信号源的产生方式不同, 因而它们之间至少是独立的。同时, 由于大地电磁源信号与噪声信号的来源不相同, 导致它们的属性也不相同, 因此, 我们可以运用独立分量分析来对大地电磁观测信号中的源信号和噪声进行分离。

运用独立分量分析进行信噪分离时要求采集到的信号的数目大于或者等于源信号的数目^[16-19], 也就是说, 必须同时获取至少两组数据, 才能将所得数据分离为信号和噪声两组数据。而在实际野外作业中, 要同时获得两组及两组以上的数据需要耗费大量的人力和物力。受文献[23-24]启发, 我们对独立分量分析算法进行了改进, 采用单通道模式对大地电磁时间序列信号进行去噪处理。首先对观测信号进行多尺度小波分解, 使信号从单道变成多道, 以满足独立分量分析对观测信号的数目需求; 然后对小波分解提取的多层高频分量进行独立分量分析, 以提取源信号和噪声。为了减小独立分量分析算法分离信号后幅值的不确定性, 本文引入动态自适应因子α来对特定源信号的权重加以限制。α的取值根据观测信号的幅值确定, 即α=max(A)/β, β=5× round { lgmax (A)}。这里A表示观测信号所有幅值的集合, round (·)表示取整数。为了尽可能地保留原始信号的基本形态，因此将小波低频分量信号加入到恢复信号中。

改进的独立分量分析算法流程如图 1所示。具体实现步骤如下:

图 1(Fig. 1)

图 1 改进的独立分量分析算法流程

1) 将观测信号进行多层小波分解, 并分别提取多层小波高频系数和某一特定层的小波低频分量, 同时观察信号的幅值, 确定自适应因子α的取值;

2) 运用M-FastICA算法对多层小波高频系数进行独立分量分析, 提取出有效源信号;

3) 按照S_de-noised=S_lW⊕0.5αS_sICA⊕0.5S_eICA的形式重构恢复信号, 其中S_de-noised为恢复信号, S_lW为特定层小波低频分量信号, ⊕为求信号的合成信号, S_sICA为特定独立分量信号, S_eICA为有效独立分量信号(有效独立分量指M-FastICA算法分离出来的独立分量, 特定独立分量指运用聚类规则从有效独立分量中挑选出来的某一个特定的独立分量)。

通过仿真实验对改进的独立分量分析算法进行了测试。仿真实验分别采用周期和非周期两种信号, 原始周期信号S1为10×sin(5t)+50×cos(3t), 原始非周期信号S2为cum(t)。噪声信号为由蒙特卡洛法构造的随机信号。

在原始信号中加入不同量的噪声, 分别采用如下7种方法进行去噪:①强制去噪法; ②中间阈值法; ③mini-max法; ④自适应默认阈值去噪法; ⑤史坦无偏似然估计去噪法; ⑥重构信号为S_de-noised=S_lW⊕S_eICA的去噪法(公式中各符号含义同前); ⑦本文改进的独立分量分析法(为了对比, 这里采用的小波基和小波分解层数分别与方法①~⑤相同)。方法①~⑤为传统的小波阈值去噪方法^[25-26]。

1) 强制去噪法。将小波分解结构中的高频系数全部置为0, 即将高频部分全部滤除。

2) 中间阈值法。设小波分解得到的各层小波系数中最大的系数为F, 最小的系数为f, 则阈值为λ=(F+f)/2。

3) mini-max法, 即最小极大方差阈值法。它是一种固定阈值的方法, 设N为小波系数的长度, 则阈值为:

4) 自适应默认阈值去噪法^[26]。这是小波阈值去噪的经典方法, 设小波系数模的中值为X, 信号的长度为N, 则阈值为$\lambda = X\sqrt {2{\rm{ln}}N} /0.6745$。

5) 史坦无偏似然估计去噪法。采用史坦无偏似然估计原理进行自适应阈值选择。

图 2和图 3分别给出了信号S1和S2加入噪声后信噪比为R_SN=15时的实验结果。

图 2(Fig. 2)

图 2 信号S1(这里S1为含噪信号)在信噪比为RSN=15时的仿真实验去噪结果 a信号S1的构成; b对信号S1进行5层小波分解后的结果; c对高频小波分量进行M-FastICA处理的结果; d对信号S1分别运用7种方法去噪后的结果

图 3(Fig. 3)

图 3 信号S2(这里S2为含噪信号)在RSN=15时的仿真实验去噪结果 a信号S2的构成; b对信号S2进行5层小波分解后的结果; c对高频小波分量进行M-FastICA处理的结果; d对信号S2分别运用7种方法去噪后的结果

为了衡量重构的恢复信号与原始信号的近似程度, 这里根据系统聚类法的基本思想, 采用相似系数来对两者的近似程度进行度量:

式中:x_ki, x_kj(k=1, 2, …, l)分别为恢复信号和原始信号; x_i和x_j分别为x_ki和x_kj(k=1, 2, …, l)的均值; c_ij为恢复信号和原始信号的相似系数; l为信号的元素个数。

表 1和表 2分别是运用7种去噪方法对信号S1和S2去噪后的结果(即恢复信号)与对应的原始信号的相似系数。考虑到大地电磁信号微弱, 获得的观测信号可能包含较多的噪声, 故这里展示了较低信噪比情况下的去噪结果。可以看出, 无论是信号S1还是信号S2, 在观测信号的信噪比较低或极低的情况下, 本文改进的独立分量分析方法去噪效果都比其它方法更为显著, 也更为稳定。这是因为本文方法结合了小波分析和M-FastICA算法各自的优势, 且所采用的恢复信号既包含低频小波系数, 又包含M-FastICA算法分离出的有效源信号和特定源信号, 同时通过自适应因子动态改变M-FastICA算法的分离效果, 使噪声消除结果受观测信号信噪比的影响大大减小。另外, 方法⑥在处理非周期信号时与方法⑦效果近似, 但在处理周期信号时明显不如方法⑦, 原因是方法⑥的恢复信号舍弃了特定源信号。

表 1(Tab. 1)

表 1 信号S1运用7种方法去噪后的相似系数 方法 RSN=-10 RSN=-5 RSN=-1 RSN=5 RSN=10 RSN=15 方法① -0.0004 -0.0001 0.0003 0.0010 0.0019 0.0027 方法② 0.3185 0.4561 0.6898 0.8792 0.9055 0.9625 方法③ 0.3130 0.4698 0.6289 0.8538 0.9382 0.9729 方法④ -0.0030 0.0015 -0.0026 0.0046 0.0034 -0.0012 方法⑤ 0.0012 0.0025 0.0040 0.0051 0.0121 0.0133 方法⑥ -0.0262 -0.0508 0.0081 0.2133 -0.3614 0.3294 方法⑦ 0.3533 0.4840 0.7093 0.9074 0.9427 0.9848

表 1 信号S1运用7种方法去噪后的相似系数

表 2(Tab. 2)

表 2 信号S2运用7种方法去噪后的相似系数 方法 RSN=-10 RSN=-5 RSN=-1 RSN=5 RSN=10 RSN=15 方法① 0.1506 0.2558 0.3594 0.5885 0.8212 0.9318 方法② 0.0763 0.1623 0.1749 0.3596 0.7143 0.7805 方法③ 0.0786 0.1346 0.1953 0.3573 0.6145 0.8147 方法④ 0.1410 0.2524 0.3613 0.5873 0.8140 0.9259 方法⑤ 0.2442 0.4980 0.6726 0.8860 0.9582 0.9841 方法⑥ 0.3890 0.6488 0.7394 0.9130 0.9678 0.9899 方法⑦ 0.3892 0.6489 0.7395 0.9131 0.9679 0.9901

表 2 信号S2运用7种方法去噪后的相似系数

将本文改进的独立分量分析去噪方法应用于实际大地电磁观测信号的噪声消除, 绘制视电阻率曲线和相位曲线, 选取含噪声较多的两个测点进行对比分析(其中测点2中的噪声含量明显高于测点1), 如图 4和图 5所示。测点1应用本文方法去噪后, 无论是视电阻率曲线还是相位曲线, 都比去噪前更加光滑、稳定, 说明本文改进的独立分量分析方法去噪效果较好。测点2因为含有的噪声多, 原始信号几乎被淹没, 但应用本文方法去噪后, 视电阻率曲线和相位曲线中的误差棒(图中红点和黑点上的竖线)相对去噪前明显变短, 甚至被消除很多, 说明本文方法能比较有效地去除信噪比极低情况下的信号噪声。

图 4(Fig. 4)

图 4 实际观测信号去噪前、后视电阻率曲线对比 a去噪前的测点1; b去噪后的测点1; c去噪前的测点2; d去噪后的测点2

图 5(Fig. 5)

图 5 实际观测信号去噪前后相位曲线对比 a去噪前的测点1; b去噪后的测点1; c去噪前的测点2; d去噪后的测点2

本文根据大地电磁测深数据量大的特点以及独立分量分析M-FastICA算法的优良性能, 结合小波分析和盲源分离的相关理论, 提出了一种改进的独立分量分析去噪方法。为了减小观测信号信噪比对去噪效果的影响, 引入动态自适应因子来限制特定独立分量的权重, 恢复信号由小波低频分量和M-FastICA算法提取的两种独立分量共同构成。模拟信号仿真实验结果表明, 本文方法的去噪性能优于传统的小波阈值去噪方法。将本文方法应用于实际大地电磁观测信号去噪处理, 无论是噪声较少还是噪声较多的情况, 都能比较有效地去除噪声, 去噪后视电阻率曲线和相位曲线比去噪前更加光滑、稳定, 误差棒相对去噪前明显变短、甚至消除。