CN108399368B - 一种人工源电磁法观测信号去噪方法 - Google Patents

Abstract

本发明公开了一种人工源电磁法观测信号去噪方法，所述方法包括：首先对含噪信号进行快速傅里叶变换，在频域信号中将工频干扰去除，进行傅里叶逆变换后回到时间域；然后对去除工频干扰后的时间域信号进行EEMD分解，衰减信号基线漂移；最后基于移不变稀疏编码进一步矫正信号的基线漂移，并去除其它人文噪声。本发明可以在不损失有用信号的前提下，将幅度远远大于有用信号的工频干扰噪声、基线漂移噪声以及脉冲等人文噪声有效去除，且允许强噪声存在于从观测开始到结束的整个过程，因此使得人工源电磁法可以在持续性、强干扰环境下观测。

Description

一种人工源电磁法观测信号去噪方法

技术领域

本发明属于勘探地球物理信号处理领域，涉及一种人工源电磁法观测信号去噪方法，尤其涉及周期信号、包括伪随机信号的去噪方法。

背景技术

人工源电磁法勘探时，发送机通过接地电极或不接地回线发送可控的电磁信号，然后在一定距离外通过接收机观测信号，据此获得地下的电性结构信息。相对于天然场源电磁法，人工源电磁法其信噪比有明显提高，但实际勘探时，接收机采集到的信号仍然不同程度的受到各种噪声的影响，因此对采集到的信号进行去噪仍然是一项必不可少的任务。

尽管现有的数据段挑选、稳健估计等数据处理方法能够取得一定效果，但如果噪声的持续时间较长，甚至存在于从开始到结束整个观测时间段时，这些方法则效果不佳。小波变换、形态滤波等方法则由于损失有用信号效果也不理想。

此外，对于普遍存在的工频干扰噪声，现有的仪器通常采用硬件陷波器或者使用DSP等器件构建数字陷波器，如果采用硬件陷波器则会因为需要的陷波器数量过多造成电路过于臃肿，且由于陷波器的性能不能理想化，不同程度的存在损伤有用信号的问题。采用数字陷波器尽管简化了电路规模，但对DSP等器件的性能有较高要求，且同样存在损伤有用信号的问题。

发明内容

为解决上述问题，本发明采用的技术方案为一种用于人工源电磁法观测信号去噪的方法。该方法包括基于快速傅里叶变换(Fast Fourier Transform,FFT)去除工频干扰噪声、基于总体平均经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)减弱基线漂移以及基于移不变稀疏编码(Shift-Invariant Sparse Coding,SISC)方法矫正基线漂移并去除其它人文噪声。

具体操作步骤如下：

(1)对观测信号进行快速傅里叶变换，得到进行快速傅里叶变换后的信号频谱序列，将观测信号频谱序列中属于工频干扰基频及其谐波所对应的频点的幅度置零，从而去除工频干扰，得到去除工频干扰的去噪信号；

(2)对去除工频干扰的去噪信号进行傅里叶逆变换，回到时间域；

(3)对回到时间域的去噪信号进行EEMD分解，得到不同频率的固有模态分量，将其中代表趋势项的低频分量去除，以减弱信号的基线漂移；

(4)采用移不变稀疏编码方法进一步矫正信号的基线漂移，并去除人文噪声。

基于快速傅里叶变换的工频干扰去除，基频为50Hz时，谐波主要包括150Hz、250Hz、350Hz等奇次谐波，有时也存在于100Hz、200Hz、300Hz等偶次谐波，因此，所述工频干扰噪声包括工频干扰的基频、奇次谐波以及偶次谐波。此外，工频干扰其频率并非是固定的，而是随电网负载的变化有所波动，因此需要将50Hz及其谐波波动范围内的频点的幅度也置零，具体的频率区间根据实际情况以及有效信号的频率进行调整。

基于快速傅里叶变换的工频干扰去除，假设观测信号x₁的采样点数为N，采样频率为f_s，观测信号x₁经快速傅里叶变换后得到的频谱序列为X₁。依照式(1)计算频率分辨率Δf：

依照式(2)确定观测信号x₁某一频率分量f在频谱序列X₁中的位置L_f：

依照式(3)确定工频干扰及其谐波在频谱序列X₁中的位置L_50n：

假设工频干扰的频率波动幅度为频率分辨率的μ倍，依照式(4)和式(5)计算工频干扰的下界d₁和上界d₂：

依据式(6)将频谱序列X₁中的工频干扰及其谐波的幅度置零：

基于EEMD减弱基线漂移的方法，EEMD将长度为N的观测信号x₂分解成w个长度仍然为N的固有模态分量即IMF，且w＝fix(log₂(N))-1，fix表示取整数部分。w个IMF从高频到低频依次排列，其中频率最低的少数几个IMF即为信号的趋势项，即基线漂移噪声。实际处理时，默认代表趋势项的IMF为4个，在移不变稀疏编码去噪后，对处理结果进行评价，如果接受处理结果，则进入下一步，否则，返回EEMD去噪，调整趋势项IMF的个数。将代表趋势项的IMF丢弃，剩余的IMF重新合并，获得矫正后的观测信号x₃。尽管EEMD相对于EMD能够显著改善模态混叠的情况，但仍无法完全避免，因此EEMD只是减弱基线漂移，并非将基线漂移完全消除。

基于移不变稀疏编码的噪声去除方法，移不变稀疏编码是一种基于数据驱动的机器学习算法，它能够从给定样本中自主学习到信号的特征结构，即特征原子，然后将学习到的特征结构组建成字典对含噪信号进行稀疏表示，在稀疏表示的过程中，与特征结构相匹配的周期性有用信号被保留下来，而与特征结构不匹配的人文噪声则被丢弃，从而实现信噪分离，将基线漂移校正并去除其他噪声。

移不变稀疏编码进一步矫正观测信号x₃的基线漂移，并去除人文噪声，对于K个长度为N的离散信号组成的集合Y＝[y₁,y₂,…y_K]^T，上标T表示矩阵的转置，y_k＝[y₁,y₂,…y_N]^T，移不变稀疏编码使用一个包含M个基的过完备字典D＝[d₁,d₂,…d_M]∈R^Q×M，Q为矩阵D的行数，将y_k表示为基与系数的卷积之和：

式(7)中*表示卷积，编码系数s_m,k∈R^P是稀疏的，即只有少量非零系数，ε为高斯白噪声，Q<N，P<N且Q+P-1＝N。

在移不变稀疏编码模型中，字典原子d_m和编码系数s_m,k均是未知的，同时求取则优化问题为非凸，很难得到稳定解，因此通过交替更新，将其转化为凸优化问题，反复迭代求取。当固定字典原子d_m不变时，基于凸优化方法求解稀疏表示系数s_m,k；当固定s_m,k不变时，基于凸优化方法求解字典原子d_m。稀疏性是两个优化问题的共同目标，评价y_k稀疏程度的代价函数为：

式(8)中||·||_F表示lF范数，β为一个约束，用于平衡重构误差与稀疏度，d_m通常需要归一化，即

所述基于移不变稀疏编码的噪声去除方法，包括如下流程：

输入：待去噪信号Y即需要进一步去噪的x₃，根据实际需求输入一个或者多个，特征原子的种类a，特征原子的宽度W＝[W₁,W₂,…W_a]，W_a为第a类特征原子在时域信号中所占采样点数，特征原子总数即稀疏度L，迭代次数z迭代次数越多，重构精度越高，z通常选择10左右即可满足要求，收敛目标e＝||ε||₁，根据精度要求设置，e越小则要求重构精度越高，取远小于1的正数；

初始化：对字典D、稀疏表示稀疏编码系数s随机赋予初始值；

重复z次以下内容：

直到达到最大的迭代次数或实现收敛目标e；

输出：过完备字典D，稀疏编码系数s，以及重构信号。

所述移不变稀疏编码，字典的学习通过以下步骤实现：保持稀疏编码系数不变，仅对基进行更新，该优化函数简化为：

式(9)中矩阵E_i,k表示第k个信号除去第i个基外的所有误差。将对第i个基的更新转化为求解关于d_i的方程组。由于d_i*s_i,k＝s_i,k*d_i，首先仅分析只有第k个信号时的优化问题，相当于求解式(10)的解：

将式(10)左边矩阵视为系数s_i,k的一种特殊的Toeplitz矩阵，则式(10)写作Toep(s_i,k)·d_i＝E_i,k。由于系数s_i,k是稀疏的，矩阵Toep(s_i,k)中的许多行为0向量，这些0向量对结果没有影响。将这些行从矩阵Toep(s_i,k)中剔除，同时将向量E_i中的对应行剔除，表示为

当同时考虑K个信号时，优化函数ψ(θ)表示为：

将式(11)简写为S·d_i＝E，根据最小二乘法可得d_i＝(S^TS)^-1(S^TE)，矩阵(S^TS)∈R^{Q ×Q}，由于Q<<N，因此最后就转化为求解一个小规模的线性方程组的解。这一求解过程通过Cholesky分解直接求取最优值。

按照随机的顺序依次对每个基进行更新，能够求取基在本轮循环中的最优解，然后依次对基d_i标准化处理

移不变稀疏编码，求取稀疏表示系数时采用OMP算法实现。设y_k为待处理信号，g_i,u为特征原子d_i平移u点补零后的原子，其长度与y_k相同，且||g_i,u||＝1，最大循环次数对应稀疏度)L_max，第L次循环后的残差为r^L，第L次循环后已选原子集合ψ_L，则OMP步骤如下：

步骤1)参数初始化，r⁰＝y_k，

L＝1；

步骤2)选出与当前残差最匹配的一个原子g_i,u，原子g_i,u满足：

步骤3)更新已选原子集合，

步骤4)根据最小二乘法计算投影系数

残差r^L＝y_k-s_Lψ_L，重构信号

步骤5)判断L是否已经达到最大值L_max，若没有，L＝L+1，并转到步骤2)循环执行；否则，输出重构信号与残差。

所述移不变稀疏编码，特征原子的种类根据实际情况设置，如果有用信号为周期信号且基频为f₀，则特征原子种类a＝1，特征原子的宽度W₁为一常数，具体通过式(13)确定：

稀疏度通过式(14)确定：

一种人工源电磁法观测信号去噪方法，适用于周期信号、含有周期特征的信号以及适用于具有相似结构的信号。对于周期信号、含有周期特征的信号、具有相似结构的信号，均属于本发明权利要求之内。如果有用信号为非周期信号但含有周期特征或者具有相似结构，则特征原子的种类a等于周期特征或者相似结构的种类，每一类周期特征或者相似结构在时域信号中所占的采样点数即为原子宽度。根据“宁多勿少”的原则从实际信号中估计出每一类周期特征或者相似结构的数量。

移不变稀疏编码，假设有用信号为非周期信号，特征原子即信号中相同或者相似的结构的种类为a，不同类型的特征原子的个数分别为L₁，L₂，…，L_a，则稀疏度通过下式确定：

L＝L₁+L₂+…L_a (15)

所述人工源电磁法信号去噪方法，基于快速傅里叶变换的工频干扰去除和基于EEMD的基线漂移校正两个步骤的先后顺序是能够进行调整的，先去除基线漂移，或先去除工频干扰，这两个步骤能够提高后续基于移不变稀疏编码的去噪效果。

基于移不变稀疏编码的噪声去除方法，适用于周期信号、含有周期特征的信号以及适用于具有相似结构的信号。对于前述周期信号、含有周期特征的信号、具有相似结构的信号，均属于本发明权利要求之内。

根据本发明实施例，本发明具有以下优点：

(1)随着城市化程度的提高，电磁法勘探受人文噪声的影响日益严重，现有的数据处理方法对于持续时间长、甚至贯穿整个观测时间段的强人文噪声通常收效甚微，使用本发明能够在不损失有用信号的前体下，从受持续性强干扰污染的观测信号中，准确分离出有用的成分，使得人工源电磁法可以在持续性、强干扰环境下进行观测，减少人工源电磁法的施工障碍；

(2)工频干扰是电磁法勘探时最普遍的人文干扰，尽管它可能并不影响最终视电阻率的计算，但是工频干扰的存在使得后续的信号处理变得更加困难，此外，它的存在也严重影响观测人员对于观测数据质量的判断。为压制工频干扰，常常采用硬件陷波器或者使用DSP等器件设计数字陷波器，但工频干扰不仅仅局限于某一个频率，而是具有很多的谐波，如果使用硬件陷波器，将需要很多个陷波器，造成硬件电路过于臃肿，并且由于硬件陷波器不可能是理想的滤波器，因此会不同程度的损伤工频干扰附近频段的有用信号；使用数字陷波器尽管能够简化电路规模，但对DSP等处理器有较高的性能要求，且同样存在损伤工频干扰附近频段的有用信号。使用本发明后，信号采集设备可以不设计抑制工频干扰的硬件陷波器或者数字陷波器，显著简化电磁法仪器接收端设备的电路规模，降低成本，还能在不损伤有用信号的前提下，将幅度远远大于有用信号的工频干扰精准的分离出来，提高工频干扰压制的效果；

(3)当信号采集设备受到温度影响或者振动等噪声的影响时，信号的整体基线容易造成明显的漂移，严重影响接收信号的低频部分，对温度产生的漂移进行标定不仅麻烦，并且标定准确度难以保证，而振动等因素造成的漂移则由于噪声源是未知的，通常无法标定。使用本发明后，不需要设计复杂的温度标定电路，即可精准的去除温度变化或者振动等因素造成的漂移。

附图说明

本发明上述的优点结合一下附图对实施例的描述做出详细说明，其中：

图1为本发明的基本流程；

图2为在实测数据中加噪声的示意图；

图3为将图2所示噪声添加到原始信号之后得到的含噪信号；

图4为使用本发明所述快速傅里叶变换法处理之后得到的信号；

图5为使用本发明所述EEMD法对图4所示信号分解后得到的固有模态分量(只显示频率最低的8个分量)；

图6为将附图5中含有明显趋势项的5个低频分量去除之后得到的信号；

图7合成数据去噪前后的对比；

图8为有效频点去噪前后的误差统计；

图9为在某地进行广域电磁法勘探时，某观测点实测的时间序列信号，该信号受到了持续性而强烈的噪声干扰，产生了明显的畸变，完全失去了伪随机信号的特征；

图10为使用本发明对图9所示信号去噪后得到的信号，显然，去噪后的信号呈现出很好的伪随机方波特征，接近于理想的伪随机方波信号，说明本发明对于实测数据同样效果明显；

图11某地进行广域电磁法勘探时，某观测点去噪前后电场值与视电阻率的曲线。

具体实施方式

图1为本发明的基本流程；

所述基于快速傅里叶变换的工频干扰去除，以中国电网基频为50Hz为例，其谐波主要包括150Hz、250Hz、350Hz等奇次谐波，有时也存在于100Hz、200Hz、300Hz等偶次谐波，因此，所述工频干扰噪声包括工频干扰的基频、奇次谐波以及偶次谐波。此外，工频干扰其频率并非是固定的，而是随电网负载的变化有所波动，因此需要将50Hz及其谐波波动范围内的频点的幅度也置零，具体的频率区间根据实际情况以及有效信号的频率进行调整。

图2-8为本发明的具体实施示意图，其中：

图2为在实测数据中加噪声的示意图，图2(a)为进行广域电磁法勘探时，实测的发送机输出信号(伪随机七频波)，图2(b)为待添加的谐波噪声，图2(c)为待添加的基线漂移噪声，图2(d)为待添加的脉冲噪声，每一种噪声其幅度均明显大于原始的实测发送机输出信号(以下简称原始信号)；

图3为将图2所示噪声添加到原始信号之后得到的含噪信号，显然，加噪后原始信号产生了严重畸变，完全失去了伪随机七频波的特征；

图4为使用本发明所述快速傅里叶变换法处理之后得到的信号，伪随机七频波信号的特征得到恢复，但存在明显的基线漂移，且含有脉冲噪声；

图5为使用本发明所述EEMD法对图4所示信号分解后得到的固有模态分量(只显示频率最低的8个分量)；

图6为将附图5中含有明显趋势项的5个低频分量去除之后得到的信号，显然基线漂移得到的明显的好转，但没有完全消除，且脉冲干扰仍然存在；

图7(b)为使用本发明所述移不变稀疏编码方法进一步处理之后得到的信号，与原始信号图7(a)相比可知，除信号的两端与原始信号有较明显的差异外，其余部分均得到准确恢复，加入的工频干扰被精准的移除，基线漂移也得到很好的校正，脉冲干扰被完全消除。经过计算，去噪后信号与原始信号的归一化互相关度(或称为相似度)为0.9926，如果将信号的第一周期与最后一个周期去除，则归一化互相关度提高到0.9990，因此可以认为本例中本发明从含有持续性强噪声的混合信号中精确的恢复了真实的原始信号；

图8为有效频点去噪前后的误差统计，去噪前，含噪信号与真实值的相对误差最大达到112.70％，而使用本发明处理后，最大误差下降到2.30％，将信号的第一个周期与最后一个周期去除后，最大误差下降到0.80％。图2～图8充分说明了本发明对数据处理结果的可靠性与有效性；

图9为在某地进行广域电磁法勘探时，某观测点实测的时间序列信号，该信号受到了持续性而强烈的噪声干扰，产生了明显的畸变，并且有严重的基线漂移，完全失去了伪随机信号的特征；

图10为使用本发明对图9所示信号去噪后得到的信号，显然，去噪后的信号呈现出很强的伪随机方波特征，具有很好的周期性，与实测发送机输出信号相似度很高，说明本发明对于实测数据同样效果明显；

图11为某地进行广域电磁法勘探时，某观测点去噪前后电场值与视电阻率的曲线，图中上部分为电场值曲线，下部分为视电阻率曲线，左边为去噪前的曲线，右边为使用本发明去噪后的曲线。从图中可知，去噪前，电场值与视电阻率均有明显的跳变，严重影响后续的反演解释，去噪后，电场值与视电阻率连续性与光滑度得到明显改善。再次说明本发明的有益效果。

综上所述，本发明提出的人工源电磁法观测信号去噪方法具有适应性强、效果好、可靠性高等优点，使得人工源电磁法可以在更加恶劣的电磁环境下施工，显著减少电磁法施工障碍。

本发明的描述和实施是说明性的，并非想将本发明的范围限制在上述实施例中。这里所披露的实施例的变形和改变是可能的，对于本领域的技术人员来说实施例的替换和等效的各种部件是公知的。在不脱离本发明的精神或本质特征的情况下，可以对这里所披露的实施例进行其它变形和改变。

Claims (4)

1.一种人工源电磁法观测信号去噪方法，其特征在于：该方法包括基于快速傅里叶变换去除工频干扰噪声、基于EEMD减弱基线漂移以及基于移不变稀疏编码方法矫正基线漂移并去除其它人文噪声；

具体操作步骤如下：

(1)对观测信号进行快速傅里叶变换，得到进行快速傅里叶变换后的信号频谱序列，将观测信号频谱序列中属于工频干扰基频及其谐波所对应的频点的幅度置零，从而去除工频干扰，得到去除工频干扰的去噪信号；

(2)对去除工频干扰的去噪信号进行傅里叶逆变换，回到时间域；

(3)对回到时间域的去噪信号进行EEMD分解，得到不同频率的固有模态分量，将其中代表趋势项的低频分量去除，以减弱信号的基线漂移；

(4)采用移不变稀疏编码方法进一步矫正信号的基线漂移，并去除人文噪声；

基于快速傅里叶变换的工频干扰去除，依照式(1)计算频率分辨率Δf：

N为观测信号x₁的采样点数，f_s为采样频率；

依照式(2)确定观测信号x₁某一频率分量f在频谱序列X₁中的位置L_f：

依照式(3)确定工频干扰及其谐波在频谱序列X₁中的位置L_50n：

设置工频干扰的频率波动幅度为频率分辨率的μ倍，依照式(4)和式(5)计算工频干扰的下界d₁和上界d₂：

依据式(6)将频谱序列X₁中的工频干扰及其谐波的幅度置零：

基于EEMD减弱基线漂移的方法，EEMD将长度为N的观测信号x₂分解成w个长度仍然为N的固有模态分量即IMF，w＝fix(log₂(N))-1，fix表示取整数部分；w个IMF从高频到低频依次排列，其中频率最低的少数几个IMF即为信号的趋势项，即基线漂移噪声；实际处理时，默认代表趋势项的IMF为4个，在移不变稀疏编码去噪后，对处理结果进行评价，如果接受处理结果，则进入下一步，否则，返回EEMD去噪，调整趋势项IMF的个数；将代表趋势项的IMF丢弃，剩余的IMF重新合并，获得矫正观测信号x₃；

移不变稀疏编码进一步矫正观测信号x₃的基线漂移，并去除人文噪声，对于K个长度为N的离散信号组成的集合Y＝[y₁,y₂,…y_K]^T，上标T表示矩阵的转置，y_K＝[y₁,y₂...y_N]^T，移不变稀疏编码使用一个包含M个基的过完备字典D＝[d₁,d₂,…d_M]∈R^Q×M，Q为矩阵D的行数，将y_k表示为基与系数的卷积之和：

式(7)中*表示卷积，编码系数s_m,k∈R^P是稀疏的，即只有少量非零系数，ε为高斯白噪声，Q<N，P<N且Q+P-1＝N；

当固定字典原子d_m不变时，基于凸优化方法求解稀疏表示系数s_m,k；当固定s_m,k不变时，基于凸优化方法求解字典原子d_m；稀疏性是两个优化问题的共同目标，评价y_k稀疏程度的代价函数为：

式(8)中||·||_F表示l^F范数，β为一个约束，用于平衡重构误差与稀疏度，d_m需要归一化，即

所述基于移不变稀疏编码的噪声去除方法，包括如下流程：

输入：待去噪信号Y即需要进一步去噪的x₃，根据实际需求输入一个或者多个，特征原子的种类a，特征原子的宽度W＝[W₁,W₂,…W_a]，W_a为第a类特征原子在时域信号中所占采样点数，特征原子总数即稀疏度L，迭代次数z迭代次数越多，重构精度越高，z选择10左右即可满足要求，收敛目标e＝||ε||₁，根据精度要求设置，e越小则要求重构精度越高，取远小于1的正数；

初始化：对字典D、稀疏表示稀疏编码系数s随机赋予初始值；

重复z次以下内容：

直到达到最大的迭代次数或实现收敛目标e；

输出：过完备字典D，稀疏编码系数s，以及重构信号；

所述移不变稀疏编码，字典的学习通过以下步骤实现：保持稀疏编码系数不变，仅对基进行更新，该代价函数简化为：

式(9)中矩阵E_i,k表示第k个信号除去第i个基外的所有误差；将对第i个基的更新转化为求解关于d_i的方程组；OMP步骤如下：

步骤1)参数初始化，r⁰＝y_k，

L＝1；

步骤2)选出与当前残差最匹配的一个原子g_i,u，原子g_i,u满足：

步骤3)更新已选原子集合，

步骤4)根据最小二乘法计算投影系数

残差r^L＝y_k-s_Lψ_L，重构信号

步骤5)判断L是否已经达到最大值L_max，若没有，L＝L+1，并转到步骤2)循环执行；否则，输出重构信号与残差；

如果有用信号为周期信号且基频为f₀，则特征原子种类a＝1，特征原子的宽度W₁为一常数，具体通过式(13)确定：

稀疏度通过式(14)确定：

2.根据权利要求1所述的一种人工源电磁法观测信号去噪方法，其特征在于：本方法适用于周期信号、含有周期特征的信号以及适用于具有相似结构的信号。

3.根据权利要求1所述的一种人工源电磁法观测信号去噪方法，其特征在于：移不变稀疏编码，假设有用信号为非周期信号，特征原子即信号中相同或者相似的结构的种类为a，不同类型的特征原子的个数分别为L₁，L₂，…，L_a，则稀疏度通过下式确定：

L＝L₁+L₂+…L_a (15)。

4.根据权利要求1所述的一种人工源电磁法观测信号去噪方法，其特征在于：人工源电磁法信号去噪方法，基于快速傅里叶变换的工频干扰去除和基于EEMD的基线漂移校正两个步骤的先后顺序是能够进行调整的，先去除基线漂移，或先去除工频干扰，这两个步骤能够提高后续基于移不变稀疏编码的去噪效果。

Images