基于DWT-EEMD的盲源分离算法在MT工频干扰消除中的应用_曹小玲.pdf

第 46 卷 第 2 期 煤田地质与勘探 Vol. 46 No.2 2018 年 4 月 COAL GEOLOGY 2. Hubei Cooperation Innovation Center of Unconventional Oil and Gas, Wuhan 430100, China; 3. School of Ination and Mathematics, Yangtze University, Jingzhou 434023, China Abstract The blind source separation algorithm based on DWT-EEMD used in noise elimination of power line in- terference for magnetotelluric is put forward in this paper. This takes advantage of the excellent character- istics of DWT, EEMD and blind source separation, and adopts blind source separation to eliminate the noise after DWT and EEMD processing. The main advantage of this lies in the DWT-EEMD model adopted and the introduction of the adaptive weighting factor, in reducing uncertainty of restored signal amplitude by independent component analysis algorithm at the same time, it has made the original signals still be well separated in condition of power frequency interference noise amplitude much higher than original signal amplitude. By the real magneto- telluric signal processing, it is showed that this makes apparent resistivity curve and phase curve both be- come smooth and steady, and it better eliminates the power line interference noise in magnetotelluric signal. Keywords DWT-EEMD; MT; blind source separation; power line interference noise 大地电磁测深Magnetotelluric, MT是一种被动 源的电磁探测技术，它通过在地表测量正交电、磁 场分量的扰动值，计算得到地下介质的电性结构信 息[1-2]。该方法具有多种技术特点，但大地电磁测深 观测的天然电磁场，具有信号弱、频带宽等优势， 易受到噪声的污染。尤其是随着近些年来国民经济 的建设和现代化工业的发展，电磁信号中的工频干 扰噪声日益增加，电力传输系统、通信过程、有线 广播以及人工电磁场等产生的噪声，逐渐成为大地 电磁信号检测中重要的噪声组成部分。这些电磁噪 ChaoXing 第 2 期 曹小玲等 基于 DWT-EEMD 的盲源分离算法在 MT 工频干扰消除中的应用 165 声信号能量可能为正常信号的数倍甚至几个数量 级，尤其是在高压输电线附近，由于干扰信号强， 几乎完全淹没了有效信号，形成似等振幅的 50 Hz 的强干扰[3-5]。远参考方法、物理去噪、小波去噪、 形态滤波等方法在去除大地电磁噪声的同时为消除 工频干扰提供了有利的手段[6-10]。近年来，为了分 析非线性和非稳态信号，N. E. Huang 等[11]、孙晖[12] 提出了一种新的时频分解算法经验模态分解 Empirical Mode Decomposition, EMD。EMD 中的 基函数和分解层次不需要事先给定，而是根据信号 特性通过迭代的方式自适应地获取，基底和分解层 次会随信号的不同而改变，因此被广泛应用于信号 去噪中。但 EMD 的模态混叠效应给其发展带来严 重阻碍。为此，集合经验模态分解Ensemble Em- pirical Mode Decomposition, EEMD[13]被提出并被 应用于信号处理。 为了进一步提高EEMD对大地电磁中的工频干 扰的去噪能力，笔者提出一种基于 DWT-EEMD 的 盲源分离算法的大地电磁工频干扰去噪方法，其主 要思想是根据工频干扰噪声的特点和 DWTDiscrete Wavelet Transation、 EEMD 及 ICAIndependent Component Analysis的分解特性， 首先利用 DWT 进 行初步的噪声消除，再用 EEMD 对大地电磁信号进 行多尺度分解，并利用各阶固有模态函数Intrinsic Mode Function, IMF和原信号的相关程度找出有用 信号的 IMF， 然后利用 ICA 对各阶 IMF 进行独立分 量分析，提取出原始信号。此方法消除了 EMD 处 理噪声信号时产生的模态混叠效应，同时克服了 ICA 方法对传感器数目必须大于或等于源信号数目 的限制，并且在工频干扰的幅值高于原始信号很多 的情况下依然能较好地分离出原始信号。 1 DWT -EEMD 和盲源分离的基本原理 鉴于在 MT 中某时刻得到的时间序列信号值可 视为随机变量在某个时刻的采样，因此我们可以将 每个时间序列信号设为随机变量信号，来进行如下 一系列的变换和处理。 1.1 DWT基本原理 在数学上，小波定义为对给定函数局部化的新 领域，小波可由一个定义在有限区域的函数 x来 构造，我们通常称 x为母小波或基本小波。一 组小波基函数  , a b x，可以通过缩放和平移基本 小波来生成 , 1 a b xb x aa      1 式中 a 为缩放参数，b 为平移参数。 信号 f x以 , a b x为基的小波变换定义为    ,, 1 ,d a ba b xb Wxff xx aa         2 当 a 和 b 分别作离散化处理后对应的式2称为 对信号 fx的离散小波变换DWT。 我们在信号处理 中使用的一般都是 DWT。在实际应用中，比较常用 的母小波函数是 dbN 小波Daubechies 小波、coifN 小波Coiflet 小波、symN 小波Symlets 小波，这里 N 为消失矩。一般来说消失矩 N 越大则滤波器滤波 效果越平滑， 但是小波分解中高频分量的零点也会相 应增加， 因此也会增加去除的有用信号， 降低去噪效 果。故消失矩的取值应根据实际需要进行选取[14]。 本文结合大地电磁信号的特点和工频干扰的特性， 选择 db3 小波N3 时的 dbN 小波作为 DWT 的小波 基函数。 1.2 EEMD基本原理 EEMD[15]方法是建立在EMD基础之上的。 EMD 是一种高效的自适应信号分解方法，它根据信号自 身特征将信号分解为一组固有模态函数IMF。 当观 测信号中某个频段的分量不连续时，就会出现模态 混叠现象从而影响 IMF 的正常分解。由于白噪声信 号具有在各个频段能量一致的特点以及均值为零的 特性，因此 EEMD 通过先在原始信号中混入白噪声 后再进行 EMD 的分解方式，以保证每个固有模态 函数时域的连续性。它在原始信号中混入足够多次 的白噪声后再对全部 EMD 分解得到的各 IMF 分量 求总体均值，以消除模态混叠现象。 1.3 盲源分离基本原理 所谓盲源分离[16]，是指仅从若干观测到的混合 信号中提取、恢复分离出无法直接观测的各个原 始信号的过程，这里的“盲”是指①源信号未知或 者不可观测；②混合系统特性事先未知或只知其少 量先验知识如非高斯性、统计独立性、循环平稳性 等。在科学研究和工程实践应用中，很多观测信号 都可以看成是多个源信号的混合，也就是说，被观 测的混合信号其实是一系列传感器的输出，而每一 个传感器接收到的是源信号的不同组合。盲源分离 的主要任务就是从获得的观测数据中恢复出我们感 兴趣的信号。 盲源分离的一般表述如下已知在生产实际中获 得的观测信号为T 12 , ,, n x kx kxkxk，要求 找 到 一 个 逆 系 统 ， 以 重 构 估 计 原 始 的 源 信 号  T 12 , ,, m s ks ksksk。源信号 sk未知，源 ChaoXing 166 煤田地质与勘探 第 46 卷 信号如何混合得到观测信号 xk也未知，这体现了 求解问题的“盲”。输出可由式3表达。 y kx ks ks kWWAC 3 式中 A 为混合矩阵，W 为分离矩阵，CWA 为混合 分离矩阵。 独立分量分析ICA是近年来发展起来的一种有 效的盲源分离技术，其关键是根据一定的优化准则建 立描述输出信号独立程度的目标函数，并设计优化算 法寻求分离矩阵，使输出信号中各分量尽可能相互独 立。从线性变换和线性空间角度，源信号为相互独立 的非高斯信号，可以看作线性空间的基信号，而观测 信号则为源信号的线性组合，ICA 就是在源信号和线 性变换均不可知的情况下，从观测的混合信号中估计 出数据空间的基本结构或者信源信号。ICA 的算法流 程如图 1 所示在信源信号 sk中各分量相互独立的 假设下，由观察信号 xk通过解混系统 B 把它们分离 开来，使输出 yk逼近 sk。ICA 中应用较多的是 FastICA 算法，本文这里亦采用此方法。 图 1 ICA 的算法流程 Fig.1 The algorithm flowchart of ICA 2 基于 DWT-EEMD 的盲源分离去噪方法 针对单独每道 MT 时间序列测量信号在进行盲 源分离方面的不足即一般要求观测信号的组数大 于或等于源信号的组数，而单独每道 MT 时间序列 一般是由多个信号混合而成的，以及实际应用方面 的需要一般来说，对同一测点同时运用多台仪器进 行测量是较难实现的，笔者提出运用 DWT-EEMD 方法对单独每道 MT 时间序列测量信号进行分析研 究。该方法主要优势在于 DWT-EEMD 模型的采用 和自适应权重因子的引入，在降低独立分量分析算 法对恢复信号的幅值的不确定性的同时，实现了单 独每道 MT 时间序列测量信号的源信号和噪声分 离，即实现了测量信号数目小于源信号数目情况下 的盲源分离。该方法的实现流程如图 2 所示。 图 2 MT 去噪方法流程图 Fig.2 The flowchart of de-noising for MT 首先，对观测信号进行 DWT 去噪。因为高频 噪声是常见噪声，故将观测信号进行多层小波分解 这里根据文献[14]所述，分解层数定为 3 层。对最 高层的高频系数强制置 0，其他层的高频系数和低 频系数均不变，以实现对高频噪声的消除。接着， 利用 EEMD 对初步去噪的信号进行分解，得到一系 列的 IMF 分量。此时，EMD 的模态混叠效应得到 消除。再对得到的 IMF 分量不包括剩余项利用主 成分分析进行分析处理，利用主成分分析帮助我们 选取出能量最大的 3 个 IMF 分量。然后，把提取的 3 个有效 IMF 分量作为输入源信号，应用 FastICA 算法恢复出各独立分量。因为考虑到工频干扰的幅 值较大能量较大，所以在对恢复出的独立分量进行 选择的时候应充分考虑其与原始信号的相关程度。 另外，由式3可知源信号的近似估计值是经过与混 合分离矩阵的线性变换得到的，其能量信息在混合 矩阵和分离矩阵相乘的过程中被改变了，表现在分 离后信号的幅值较源信号幅值会发生改变。所以对 于恢复出的信号，必须对其幅值有一定的约束。因 此，为了降低独立分量分析算法对恢复信号的幅值 的不确定性，我们考虑引入一个自适应的权重因子 来对分离后的信号幅值加以约束。考虑到在分离后 幅值的取值范围，这里权重因子 β 的取值为 tmax /max,1,2,, ii RABiN，这里 Ai为经 ChaoXing 第 2 期 曹小玲等 基于 DWT-EEMD 的盲源分离算法在 MT 工频干扰消除中的应用 167 独立分量分析分离后提取信号的幅值， Bi为观测信号 的幅值，Rt为取整操作，N 为信号长度。在对大地电 磁测点 5 道信号进行相同形式的幅值处理后， 所得到 的误差在进行功率谱及阻抗计算时会相互抵消。 3 仿真模拟 在实际问题中，当我们无法区分在某区间内取 值的随机变量取不同值的可能性有何不同时，就可 以假定该随机变量服从这一区间上的均匀分布，因 此这里仿真信号中的原始信号通过蒙特卡洛方法构 造为如图 3a 所示幅值为 20。 因为工频干扰一般为 50 Hz 及其谐波，故构造工频干扰信号噪声信号 为sin2π 50 NBt，其中 B 为幅值这里为便于研 究，将 B 依次取为不同的值，工频干扰信号噪声 信号如图 3b 所示限于篇幅， 图中显示的仅是 B80 时的情况。将工频干扰噪声加入到原始信号中得到 的含噪信号如图 3c 所示。含噪信号经 DWT 初步去 噪处理后所得到的信号如图 3d 所示。 图 3 DWT 处理前后的信号对比图 Fig.3 The comparison of signal before and after DWT 对图3d信号进行EEMD分解的结果如图4a所示， 对该分解结果进行主成分分析的结果如图 4b 所示。图 4a 中信号 i1i11 是各个固有模态函数IMF，r 为残量 信号。图 4b 中 p1p11 为对图 4a 中 i1i11 进行主成分 分析PCA得到的各个主成分分量。因为工频干扰幅值 和能量都较大，因此它被分解到不同阶的 IMF 中，当 对这些 IMF 进行主成分分析之后，它会被大量的包含 于前几个主要的主成分分量之中。 提取主成分的前 3 个主要分量进行 ICA 处理的 结果如图 5 所示C1、C2、C3 均是分解出的不同独 图 4 EEMD 分解结果和 PCA 分析结果 Fig.4 The results of EEMD and PCA processing 立分量。从图中可以发现，原始信号已被提取出来 图 5a 中 C1 分量，C2，C3 为工频干扰对应的分量。 但 C1 的幅值与原始信号有差距。再根据观测信号的 幅值范围来确定权重因子 β，并将权重因子作用于 C1 得到最后的恢复信号如下图 6c 所示。 与原始信号 图 6a比较可以看出，信号幅值得到较好的恢复，且 去噪效果好于 DWT 去噪图 6b 为其去噪结果。 为了衡量提取的恢复信号与原始信号的近似程 度，根据系统聚类法的基本思想，这里采用相关系 数来对两者的近似程度进行度量。 ChaoXing 168 煤田地质与勘探 第 46 卷 图 5 ICA 分解结果 Fig.5 The result of decomposing by ICA 图 6 去噪结果比较 Fig.6 The comparison of de-noising results   1 1/2 2 2 11 n ij kikj k ij nn ij kikj kk xxxx c xxxx               4 式中 1,2,, ki xkn，1,2,, kj xkn分别为时间 序 列 ， ix为1,2,, ki xkn的 均 值 ， jx为 1,2,, kj xkn的均值， ij c为这两个时间序列的相 关系数。表 1 为当工频干扰的幅值 B 分别取不同的 值时，采用 3 种方法去噪所得的恢复信号与原始信 号的相关系数情况。其中去噪方法 1 是指强制性小 波去噪法；去噪方法 2 是指自适应性小波去噪法； 方法 3 为笔者提出的基于 DWT-EEMD 的盲源分离 算法。 表 1 3 种方法的去噪性能比较 Table 1 The de-noising perance comparison of the three s 相关系数 幅值 B 方法1 方法2 方法3 60 0.241 4 0.063 4 0.825 4 70 0.232 8 0.058 6 0.822 3 80 0.203 8 0.054 6 0.821 1 90 0.173 3 0.053 2 0.818 7 100 0.140 1 0.051 5 0.816 4 110 0.125 8 0.050 1 0.814 2 120 0.123 3 0.048 2 0.813 9 130 0.122 0 0.045 6 0.813 0 140 0.116 7 0.041 6 0.812 2 150 0.110 2 0.039 8 0.810 7 160 0.106 9 0.032 5 0.809 8 170 0.101 1 0.031 1 0.806 9 180 0.098 8 0.027 1 0.805 1 190 0.090 9 0.025 8 0.804 5 200 0.067 3 0.006 4 0.802 1 从表 1 中可以发现，随着工频干扰的幅值的不 断增加，3 种方法提取的恢复信号与原始信号的相 似程度都在降低，但下降速度不快。两种小波去噪 法作为信号去噪的常用方法对低信噪比的信号中的 工频干扰噪声的去噪能力不高，去噪效果都不甚理 想，最后得到的恢复信号与原始信号的相似程度都 很低80，而且非常稳定，不会随着噪声 幅值的增加而显著下滑。所以，该方法对于去除工 频干扰噪声是比较理想和有效的。 4 实际大地电磁信号去噪处理 选取在某地测得的大地电磁测点进行研究，研 究区地理位置位于北纬 31、东经 114，附近有一 个大型工业区及若干电力线，测点受到较大的工频 干扰影响。截取其中 3 段时间序列，显示如图 7 所 示，图 7a图 7c 的时间起点分别是052500， 102500，201500。从图 7a 可以看出，在 052500 及之后 2 s 内工频干扰对 Ex、Ey、Hy影响很大，使 曲线形态形成了工频干扰的正弦曲线，原始信号几 乎被其淹没，特别是对于 Ex，其工频干扰的幅值很 大，且几乎不受其他信号影响；在 102500 及之后 2 s 内工频干扰对 Ex、Ey、Hy仍有影响，但影响明 显减小，此时随着各种脉冲干扰的产生，削弱了工 ChaoXing 第 2 期 曹小玲等 基于 DWT-EEMD 的盲源分离算法在 MT 工频干扰消除中的应用 169 图 7 原始大地电磁时间序列 Fig.7 The original time series of magnetotelluric 频干扰的正弦曲线形态，也同时使噪声形态更加复 杂图 7b；在 201500 及之后的 2 s 内工频干扰对 Ex和 Hy有较大影响，对 Ey的影响明显减弱图 7c。 将该测点进行基于 DWT-EEMD 的盲源分离算 ChaoXing 170 煤田地质与勘探 第 46 卷 法工频干扰去噪处理后所得到的时间序列如下图 8 所示。从图 8 可以看出，经过该方法去噪后，工频 干扰噪声得到了很好的消除，图中几乎看不出工频 干扰正弦曲线波形的存在。 图 8 去噪后的大地电磁时间序列依次与图 7 中时间序列对应 Fig.8 The time series of magnetotelluric after de-noising in sequence corresponding to the time series in Fig.7 ChaoXing 第 2 期 曹小玲等 基于 DWT-EEMD 的盲源分离算法在 MT 工频干扰消除中的应用 171 将去噪前后的时间序列及其对应的参数文件进 行处理得到去噪前后的大地电磁响应曲线如图 9 所 示。图 7 中时间序列最后得到的大地电磁响应曲线 是图 9a 和图 9b即去噪前， 图 8 中时间序列最后得 到的大地电磁响应曲线是图 9c 和图 9d即去噪后。 从图 9a 和图 9b 可以看出，因为工频干扰的影响， 使视电阻率曲线在 50 Hz 处发生了断裂，且使整个 曲线多处不连续，也使相位曲线变得杂乱无章。通 过本文方法的去噪处理后，不但消除了断裂带，而 且使整个视电阻率曲线变得光滑、连续，也使相位 曲线不再杂乱无章。从图中可以看出，基于 DWT- 图 9 DWT-EEMD 法去噪前后的大地电磁响应曲线 Fig.9 The magnetotelluric response curve before and after DWT-EEMD de-noising EEMD 的盲源分离算法对工频干扰进行了有效的消 除，使因为工频干扰导致的大地电磁响应曲线的断 裂及不连续都得到了很好的修复，使大地电磁响应 曲线变得光滑流畅，为之后的反演处理及地质解释 提供了良好的数据基础。 5 结 论 a. 基于 DWT-EEMD 的盲源分离算法在大地电 磁工频干扰去噪中的应用时，充分利用了 DWT、 EEMD 以及盲源分离的优良特性，在消除了 EMD 处理噪声信号时产生的模态混叠效应的同时，克服 了 ICA 方法对传感器数目必须大于或等于源信号数 目的限制。 b. 该方法通过引入一个自适应的权重因子来 对盲源分离后的信号幅值加以约束，从而降低了独 立分量分析算法对恢复信号的幅值的不确定性。 c. 该方法在工频干扰的幅值高于原始信号幅 值很多的情况下依然能较好地分离出原始信号，即 它在处理低信噪比信号中的工频干扰噪声时亦能取 得较好的效果。 d. 本次提出基于 DWT-EEMD 的盲源分离算法 在大地电磁工频干扰去噪中的应用，旨在为大地电 磁信号中的工频干扰噪声的消除提供一种新的方法 和思路。 参考文献 [1] 石应俊，刘国栋，吴广耀，等. 大地电磁测深法教程[M]. 北 京地震出版社，1985. [2] 陈乐寿，刘任，王天生. 大地电磁测深资料处理与解释[M]. 北京石油工业出版社，1989. [3] 孙洁，晋光文. 大地电磁测深资料的噪声干扰[J]. 物探与化 探，2000，242119–127. SUN Jie，JIN Guangwen. Noise interference of magnetotelluric data[J]. Geophysical and Geochemical Exploration， 2000， 242 119–127. [4] 梁生贤，张胜业，黄理善，等. 大地电磁勘探中的电力线工频 干扰[J]. 物探与化探，2012，365813–816. LIANG Shengxian，ZHANG Shengye，HUANG Lishan，et al. Power frequency interference of electric power lines in mag- netotelluric prospecting[J]. Geophysical and Geochemical Ex- ploration，2012，365813–816. [5] 徐志敏， 汤井田， 强建科. 矿集区大地电磁强干扰类型分析[J]. 物探与化探，2012，362214–219. XU Zhimin，TANG Jingtian，QIANG Jianke. Analysis of the type of magnetotelluric strong interference in ore-concentrated area[J]. Geophysics and Geochemical Exploration，2012， 362214–219. [6] 严良俊，胡文宝，陈清礼，等. 远参考 MT 方法及其在南方强 ChaoXing 172 煤田地质与勘探 第 46 卷 干扰地区的应用[J]. 石油天然气学报，1998，20434–38. YAN Lingjun，HU Wenbao，CHEN Qingli，et al. Application of remote reference MT to noisy area[J]. Journal of Oil and Gas， 1998，20434–38. [7] 胡家华，陈清礼，严良俊，等. MT 资料的噪声源分析及减小 观测噪声的措施[J]. 石油天然气学报，1999，21469–71. HU Jiahua， CHEN Qingli， YAN Liangjun， et al. Analyzing noise sources of MT data and minimizing measurement noise[J]. Jour- nal of Oil and Gas，1999，21469–71. [8] 陈清礼，胡文宝，苏朱刘，等. 长距离远参考大地电磁测深试 验研究[J]. 石油地球物理勘探，2002，372145–148. CHEN Qingli，HU Wenbao，SU Zhuliu，et al. Study for long-distant and far-referential MT[J]. Journal of Petroleum Geophysical Prospecting，2002，372145–148. [9] 苏朱刘，胡文宝，张翔. 电磁资料中的物理去噪法[J]. 工程地 球物理学报，2004，12110–115. SU Zhuliu，HU Wenbao，ZHANG Xiang. Noise reduction by means of physical properties of electromagnetic data[J]. Journal of Engineering Geophysics，2004，12110–115. [10] 汤井田，李晋，肖晓，等. 数学形态滤波与大地电磁噪声压制[J]. 地球物理学报，2012，5551784–1793. TANG Jintian，LI Jin，XIAO Xiao，et al. Mathematical mor- phology filtering and noise suppression of magnetotelluric sounding data[J]. Chinese Journal of Geophysics，2012，555 1784–1793. [11] HUANG N E，SHEN Zheng，LONG S R，et al. The empirical mode decomposition and the Hilbert spectrum for non-linear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London A，1998，454903–995. [12] 孙晖. 经验模态分解理论与应用研究[D]. 杭州浙江大学， 2005 [13] WU Zhaohua，HUANG N E. Ensemble empirical mode de- composition A noise-assisted data analysis [J]. Advances in Adaptive Data Analysis，2009，111–41. [14] 凌振宝，王沛元，万云霞，等. 强人文干扰环境的电磁数据小 波去噪方法研究[J]. 地球物理学报，2016，5993436–3447. LING Zhenbao，WANG Peiyuan，WAN Yunxia，et al. The research of wavelet denoising of EM data in strong hu- manistic interference environment[J]. Chinese Journal of Geo- physics，2016，5993436–3447. [15] 毕凤荣，陆地，邵康，等. 基于 EEMD-ICA-CWT 的装载机室 内噪声盲源分离和识别[J]. 天津大学学报，2015，489 804–810. BI Fengrong， LU Di， SHAO Kang， et al. Blind source separation and identification of loader indoor noise based on EEMD-ICA-CWT[J]. Journal of Tianjin university， 2015， 489 804–810. [16] 余先川，胡丹. 盲源分离理论与应用[M]. 北京科学出版社， 2011. 责任编辑 聂爱兰 ChaoXing