CN110146929A - 基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法 - Google Patents

Abstract

本发明提供一种基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法，包括：首先利用数学形态滤波从含噪的大地电磁时间序列中提取出粗糙的低频有效信号；然后利用互补集合经验模态分解对粗糙的低频有效信号进行圆滑处理，获得准确的低频有效信号，含噪的大地电磁时间序列减去提取到的低频有效信号即得到含噪的高频信号；最后通过设计合适的过完备字典，利用压缩感知重构算法对含噪的高频信号进行信噪分离，得到去噪后的高频有效信号；低频有效信号与高频有效信号相加即得到全频带的大地电磁有效信号；本发明在较好的保留大地电磁有效信号的前提下，去除低频大地电磁数据中的强人文噪声，显著提高大地电磁数据的信噪比，改善视电阻率与相位曲线。

Description

基于过完备字典与压缩感知重构算法的低频大地电磁数据去 噪方法

技术领域

本发明属于地球物理信号处理领域，涉及一种基于过完备字典与压缩感知重构算法的大地电磁数据去噪方法，尤其是10Hz以下低频段大地电磁数据的去噪方法。

背景技术

大地电磁法具有探测深度大、无需人工场源等优点，被广泛应用于地球深部探测。然而天然的大地电磁信号具有很强的随机性并且幅度微弱，极易受到噪声的影响。随着人类文明的不断进步，各类人文噪声日益增强，获得高质量的大地电磁信号越来越难。

面对严峻的人文噪声污染，通用的做法是设置远参考站并进行Robust稳健估计。然而，良好的远参点与高质量的大地电磁数据同样难以获取，远参考处理所带来的有益效果十分有限。由于人文噪声常常是持续性的，而Robust估计法要求有一定比例的高质量数据，因此Robust稳健估计也常常无法奏效。为此，国内外学者纷纷从时间域进行处理，小波变换、经验模态分解(Empirical Mode Decomposition，EMD)、数学形态滤波均被引入大地电磁信噪分离，但这些方法直接从含噪信号中提取人文噪，极易损失低频信号，难以用于10Hz以下低频段大地电磁信号的处理。

专利CN201610618524.6将经验模态分解与数学形态滤波组合，先利用经验模态分解将含噪信号分解成低阶模态分量(即高频信号成分)与高阶模态分量(即低频信号成分)；然后利用形态滤波与自适应阈值去噪方法对低阶模态分量进行去噪得到去噪后的高频信号，利用平滑滤波法消除高阶模态分量中的基线漂移，得到去噪后的低频信号。然而利用平滑滤波法对高阶模态分量(低频部分)进行去噪，极易损失10Hz以下的低频有效信号，因此该方法主要用于10Hz以上的高频大地电磁信号(或称为音频大地电磁信号)的处理。此外，经验模态分解(EMD)存在严重的模态混淆问题，导致不同频率成分的信号存在于同一个模态分量中，或者同一频率成分的信号分布在多个模态分量，无法准确区分高频与低频有效信号。

李广(也是本申请的第一发明人)和汤井田等(Li et al，2017，AppliedGeophysics；汤井田等，2017，地球物理学报)将稀疏表示引入大地电磁信号处理，通过设计与人文噪声匹配但对有效信号不敏感的冗余字典，实现仅对含噪片段进行处理；李晋等(CN201710938046.1、201810858136.4、201810858135.X)先将原始信号辨识为高质量片段和强干扰段，然后仅对强干扰段进行去噪。这两类方法能够有效的保护高质量片段中的低频信号，但是对于强干扰段的处理仍然是直接提取人文噪声，损失了强干扰段中的部分低频有效信号，因此在处理10Hz以下的低频大地电磁数据时难以取得良好的效果。

发明内容

为解决上述问题，本发明创造性地将MMF、CEEMD以及稀疏表示结合，增加低频有效信号提取的环节，先从含噪信号中分离出低频有效信号予以保护，然后再从剩余信号中提取人文噪声，从而提出一种基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法。

数学形态滤波(Mathematical Morphology Filtering，MMF)是由数学形态学发展而来的一种非线性滤波技术，它包括膨胀和腐蚀两个基本运算。腐蚀运算是收缩的过程，用来消除物体的边界点，减少峰值。膨胀运算是扩张的过程，用来填平边界不平滑的凹陷部分，增大谷值。形态滤波可以用于冲击类成分与低频有效信号的分离。在大地电磁信号处理领域，形态滤波通常被用于提取噪声(例如论文Li et al.，2017，Applied Geophysics)，本发明则相反，利用形态滤波提取大地电磁低频有效信号。

互补集合经验模态分解(Complementary Ensemble Empirical ModeDecomposition，CEEMD)是在Empirical Mode Decomposition，EMD和Ensemble EmpiricalMode Decomposition，EEMD的基础上改进而来的一种非线性、非平稳信号处理方法，它利用白噪声均匀分布的特性，通过添加成对且互补的白噪声，并多次分解取平均值，有效解决了EMD存在的模态混淆以及EEMD残留噪声的问题。CEEMD将信号分解成有限个数不同频率的固有模态分量(Intrinsic Mode Function，IMF)，能够将混合信号的高频部分和低频部分有效分离。

稀疏表示包括过完备字典设计和信号重构两部分，通过设计合适的过完备字典，利用信号重构算法，可以在稀疏度等条件的约束下，准确分离出复杂信号中与过完备字典相匹配的成分。

发明人通过不断地探索，将MMF、CEEMD以及稀疏表示有机结合，用于大地电磁数据去噪，尤其是10Hz以下低频段大地电磁数据的去噪，取得了出人意料的效果，有效解决了传统方法损失低频有效信号的问题，显著提高了数据质量。形成了完整的技术方案。

本发明技术方案，基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法，其步骤如下：

步骤1、利用数学形态滤波将含噪的大地电磁时间序列A分解成粗糙的低频有效信号B和残差B_r。

步骤2、利用CEEMD将粗糙的低频有效信号B分解成准确的低频有效信号C和残差C_r，残差B_r与C_r相加得到含噪的高频信号X。

步骤3、利用稀疏分解将X分解成高频有效信号Y_h和人文噪声Y_c，高频有效信号Y_h和低频有效信号C相加得到全频带的大地电磁有效信号Y。

进一步的，步骤1中所述数学形态滤波，其操作如下：

设A＝a(n)＝[a₁，a₂，…，a_N]^T是长度为N的待处理含噪信号，g＝g(m)＝[g₁，g₂，…，g_M]^T是长度为M的结构元素，M要求略大于大地电磁信号中人文噪声所占的采样点个数，处理电通道信号时，结构元素的长度M取值范围为80～150。处理磁通道信号时，结构元素的长度M取值范围为50～100。用数学形态滤波提取到的粗糙的低频有效信号B为：

B＝ψ_OC(CO)(g)＝[OC(A)+CO(A)]/2。(1)

上式中OC(A)和CO(A)分别表示形态开-闭滤波器和形态闭-开滤波器：

式中“ο”表示开运算，“·”表示闭运算，开运算和闭运算定义为：

式中“Θ”和分别表示腐蚀和膨胀运算：

进一步的，所述步骤2中CEEMD的分解步骤如下：

2.1给信号B分别添加白噪声W和-W，得到加噪后的信号B₁和B₂：

B₁＝B+W， (8)

B₂＝B-W， (9)

分别对B₁和B₂进行EMD分解，其步骤如下：

2.2计算加噪后的信号B₁(B₂)的全部极大值与极小值，采用三次样条插值法分别将极大值拟合成上包络线，将极小值拟合成下包络线。

2.3将极大值包络和极小值包络相加后取均值，得到均值包络m⁺(m^-)。

2.4提取临时局部震荡h⁺＝B₁-m⁺(相应的h^-＝B₁-m^-)。

2.5对h⁺(h^-)重复步骤2.2～2.6，直到均值包络接近于零则停止。此时，h⁺(h^-)即一个固有模态分量，记为I⁺(I^-)。

2.6计算残差r⁺＝B₁-I⁺(相应的r^-＝B₂-I^-。)

2.7对残差r⁺(r^-)重复执行步骤2.2～2.6，得到第二个IMF，重复执行z次(固有模态分量的个数)或一直执行直到残差成为无法再分解的单调信号时停止。得到分解结果如下：

2.8利用不同的白噪声，对步骤2.1～2.7重复执行α次(即集合次数)，并取平均值，得到最终的分解结果如下：

其中：

本发明中，固有模态分量的个数z取10，集合次数α取10，经过大量实践，将低频有效信号C确定为第5至第10个固有模态分量以及残差之和：

相应的，可得到残差C_r为：

C_r＝B-C， (16)

含噪的高频信号X为：

X＝B_r+C_r。 (17)

进一步的，所述步骤3包括过完备字典设计以及信号重构两部分，所述过完备字典，其表达式为：

其中d_γ是由τ和a定义的归一化的原子，且||d_γ||＝1，τ∈(1，N)表示持续时间，a表示幅度。

所述压缩感知重构算法为子空间追踪算法(Subspace Pursuit，SP)。设D＝{d_γ}_γ∈Γ是设计的过完备字典，d_γ是字典中的某个原子，γ是原子的序号，Γ是γ的集合，K是稀疏度，R^l是迭代l次之后的残差，Λ_l是候选原子集合，ψ_l是已选原子集合，则子空间追踪算法步骤如下：

输入：K，D，X，l，Λ_l，ψ_l，R^l。

初始化设置：l＝1，R⁰＝X，

迭代：

1)更新候选原子集合Λ_l：

Λ_l＝Λ_l-1∪{使得表达式|<R^l，D>|值最大的K个原子}。

2)计算投影系数u_l：

3)更新已选原子集合Ψ_l：

Ψ_l＝{使得u_l最大的K个原子}。

4)更新重构信号Y_l与残差R^l：

R^l＝X-Y_l。 (21)

5)如果||R^l||₂＞||R^l-1||₂，则Ψ_l＝Ψ_l-1且停止迭代。否则，l＝l+1，并返回步骤1)继续迭代。

输出：重构信号Y_c＝Y_l。

重构信号Y_c即为稀疏分解提取到的人文噪声，含噪的高频信号X减去人文噪声Y_c即得到去噪后的高频有效信号Y_h：

Y_h＝X-Y_c。 (22)

将低频有效信号C与高频有效信号Y_h相加，即得到去噪后的全频带的大地电磁有效信号Y：

Y＝C+Y_h。 (23)

本发明具有以下优点：

1、本发明在大地电磁信噪分离方法方面，与现有技术易损失低频有效信号相比，在去噪前增加了保护低频有效信号的环节，有效避免了低频有效信号的损失。在保留低频有效信号的前提下，准确去除强烈的文化噪声，在对10Hz以下低频大地电磁数据去噪时，取得远远优于直接采用MMF、小波、经验模态分解等现有方法的效果。

2、在大地电磁信号处理领域，形态滤波通常被用于提取噪声，本发明将形态滤波用于提取低频有效信号，并与有效解决了模态混淆问题的CEEMD算法结合，提供了一种基于MMF与CEEMD的低频有效信号提取方法。该方法能够将大地电磁信号中低频缓变部分准确提取出来，为后续精准的信噪分离提供了良好的基础。低频缓变部分的提取方法，不仅适用于大地电磁信号，也可以用于地震等其他信号的低频缓变部分的提取。

3、本发明在基于过完备字典与压缩感知重构算法的大地电磁信噪分离方法方面，通过设计与噪声匹配但是对有效信号不敏感的过完备字典，并利用压缩感知重构算法提高效率，能够快速准确的识别出人文噪声的位置，并将其分离出来。基于过完备字典与压缩感知重构算法的去噪方法，仅对过完备字典相匹配的成分敏感，因此能更好的保留有效信号。

4、总之，本发明能够有效的用于10Hz以下、尤其是1Hz以下低频段的电磁数据的去噪。频率越低，探测深度越大，因此本发明能够为深部地球物理探测提供更加可靠的数据，为提高深地探测的分辨率奠定了良好的基础。

附图说明

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

图2为低频有效信号提取实验结果对比图。

图3为不同方法对合成数据的信噪分离结果对比图。

图4为合成数据处理前后视电阻率-相位曲线对比图。

图5为实测数据处理前后的时间序列效果图(脉冲噪声)。

图6为实测数据处理前后的时间序列效果(方波噪声)。

图7为实测数据处理前后视电阻率-相位曲线图。

具体实施方式

本发明结合附图和实施例的描述做出详细说明。

实施例

附图1为本发明的基本流程。首先，利用数学形态滤波将含噪的大地电磁时间序列A分解成粗糙的低频有效信号B和残差B_r。然后，利用CEEMD将粗糙的低频有效信号B分解成准确的低频有效信号C和残差C_r，残差B_r与C_r相加得到含噪的高频信号X。最后，利用设计的过完备字典以及压缩感知重构算法根据稀疏分解原理将X分解成高频有效信号Y_h和人文噪声Y_c，高频有效信号Y_h和低频有效信号C相加得到全频带的大地电磁有效信号Y。

附图2所示为低频有效信号提取的实例。图2中(a)是合成的含噪信号，图2中(b)是通过步骤1和步骤2提取的大地电磁低频有效信号，图2中(c)是含噪的高频信号X，主要由文化噪声和高频有效信号组成。由图2可知，本发明可以从含噪信号中准确的提取出低频缓变部分。这样操作的目的是将低频信号保护起来，使得后续的信噪分离不会损伤低频有效信号。

附图3为不同方法对合成数据的信噪分离结果。图3中(a)为实测的高信噪比大地电磁信号。图3中(b)为合成的含噪信号。图3中(c)为直接采用形态滤波法去噪后的信号。图3中(d)为直接采用小波阈值法去噪后的信号。图3中(e)为直接使用OMP法去噪后的信号。图3中(f)为本发明所提MMF-CEEMD-SP法去噪后的信号。如图3所示，形态滤波成功地消除了方波噪声，但留下了尖峰。小波阈值方法成功地消除了方波噪声和脉冲噪声，但丢失了许多有用的信号。OMP算法的原理表明，识别噪声的基础是噪声和原子的最大内积。因此，OMP方法有时会将具有大振幅的信号误认为是噪声。另一方面，一些幅度小的噪声可能被错误地判断为有用信号。MT低频缓变部分的幅度相对较大。直接使用OMP方法降低噪声容易导致过度处理或处理不足。如图3(e)所示，添加的噪声尚未完全消除(例如205秒的位置)，并且在某些地方(例如100秒的位置)已经去除了一些有用的信号。只有本发明所提方法不仅能准确地消除所有噪声，而且没有明显的有效信号的损失。

附图4(a1，a2)所示为使用原始无噪声数据计算的视电阻率和相位曲线，除低于0.001Hz以下有个别频点有小幅度的跳变以为，整条曲线平滑且连续变化，表现出较高的数据质量。如图4(b1，b2)所示，使用加噪后的合成数据获得的视电阻率和相位曲线在1Hz和0.001Hz之间严重失真。特别是0.1Hz左右的曲线，方差非常大。如图4(c1，c2)所示为Robust稳健估计法处理后的结果，由于噪声是持久的，因此稳健的统计方案导致更差的结果。将本发明所述完整的去噪方法应用于合成噪声数据，得到了视电阻率和相位曲线如图4(d1，d2)所示。除最后一个频率外，本发明所提方法处理后得到的结果与图4(a1，a2)基本一致。这个案例证明了我们方法的有效性和准确性。

附图5和附图6中的MT数据分别受到脉冲噪声和方波噪声的污染。噪声的幅度明显大于有效信号。图5(a)和图6(a)为去噪前的信号，图5(b)和图6(b)为去噪后的信号。经过本发明所提方法处理后，消除了各种文化噪声，且保留了低频缓变部分。附图7所示为青海某实测点的处理结果，其中图7(a1，a2)为使用本发明去噪前的数据计算的视电阻率和相位曲线，显然，该曲线在5Hz和0.05Hz之间严重失真。应用本发明所提去噪方法处理后的数据计算得到如图7(b1，b2)视电阻率和相位曲线，显然，曲线整体的连续性和平滑度得到极大地改善。

综上所述，本发明可以有效去除低频大地电磁数据中的强人文噪声，显著提高数据质量，取得明显优于现有方法的效果。

本发明的描述和实施是说明性的，并非想将本发明的范围限制在上述实施例中。这里所披露的实施例的变形和改变是可能的，凡在本发明的精神和原则之内，做出的简单修改、等同替换等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于过完备字典与压缩感知重构算法的低频大地电磁(Magnetotelluric，MT)数据去噪方法，其特征在于，包括以下步骤：

步骤1、利用数学形态滤波(Mathematical Morphology Filtering，MMF)将含噪的大地电磁(Magnetotelluric，MT)时间序列A分解成粗糙的低频有效信号B和残差B_r，即A＝B+B_r；

步骤2、利用互补集合经验模态分解(Complementary Ensemble Empirical ModeDecomposition，CEEMD)将粗糙的低频有效信号B分解成准确的低频有效信号C和残差C_r，残差B_r与C_r相加得到含噪的高频信号X；

步骤3、利用稀疏分解(Sparse Decomposition，SD)将X分解成高频有效信号Y_h和人文噪声Y_c，高频有效信号Y_h和低频有效信号C相加得到全频带的大地电磁有效信号Y。

2.根据权利1所述基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法，其特征在于，所述步骤3中稀疏分解包括过完备字典设计和信号重构两部分；所述过完备字典，即解析型字典，通过数学表达式确定；所述信号重构通过以匹配追踪算法为基础的压缩感知重构算法完成。

3.根据权利1所述基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法，其特征在于，步骤1中所述数学形态滤波，其操作如下：

设A＝a(n)＝[a₁，a₂，…，a_N]^T是长度为N的待处理含噪信号，g＝g(m)＝[g₁，g₂，…，g_M]^T是长度为M的结构元素，M要求略大于大地电磁信号中人文噪声所占的采样点个数，处理电通道信号时，结构元素的长度M取值范围为80～150；处理磁通道信号时，结构元素的长度M取值范围为50～100；用数学形态滤波提取到的粗糙的低频有效信号B为：

B＝ψ_OC(CO)(g)＝[OC(A)+CO(A)]/2； (1)

上式中OC(A)和CO(A)分别表示形态开-闭滤波器和形态闭-开滤波器：

式中表示开运算，“·”表示闭运算，开运算和闭运算定义为：

式中“Θ”和分别表示腐蚀和膨胀运算：

4.根据权利1所述基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法，其特征在于，所述步骤2中CEEMD的分解步骤如下：

2.1给信号B分别添加白噪声W和-W，得到加噪后的信号B₁和B₂：

B₁＝B+W， (8)

B₂＝B-W， (9)

分别对B₁和B₂进行EMD分解，其步骤如下：

2.2计算加噪后的信号B₁(B₂)的全部极大值与极小值，采用三次样条插值法分别将极大值拟合成上包络线，将极小值拟合成下包络线；

2.3将极大值包络和极小值包络相加后取均值，得到均值包络m⁺(m^-)；

2.4提取临时局部震荡h⁺＝B₁-m⁺(相应的h^-＝B₁-m^-)；

2.5对h⁺(h^-)重复步骤2.2～2.6，直到均值包络接近于零则停止；此时，h⁺(h^-)即一个固有模态分量，记为I⁺(I^-)；

2.6计算残差r⁺＝B₁-I⁺(相应的r^-＝B₂-I^-)；

2.7对残差r⁺(r^-)重复执行步骤2.2～2.6，得到第二个IMF，重复执行z次(固有模态分量的个数)或一直执行直到残差成为无法再分解的单调信号时停止；得到分解结果如下：

2.8利用不同的白噪声，对步骤2.1～2.7重复执行α次(即集合次数)，并取平均值，得到最终的分解结果如下：

其中：

固有模态分量的个数z取10，集合次数α取10，经过大量实践，将低频有效信号C确定为第5至第10个固有模态分量以及残差之和：

相应的，可得到残差C_r为：

C_r＝B-C， (16)

含噪的高频信号X为：

X＝B_r+C_r。 (17)

5.根据权利2所述基于过完备字典与压缩感知重构算法的低频大地电磁数据去噪方法，其特征在于，所述步骤3包括过完备字典设计以及信号重构两部分，所述过完备字典，其表达式为：

其中d_γ是由τ和a定义的归一化的原子，且||d_γ||＝1，τ∈(1，N)表示持续时间，a表示幅度；

所述压缩感知重构算法为子空间追踪算法(Subspace Pursuit，SP)；设D＝{d_γ}_γ∈Γ是设计的过完备字典，d_γ是字典中的某个原子，γ是原子的序号，Γ是γ的集合，K是稀疏度，R^l是迭代l次之后的残差，Λ_l是候选原子集合，Ψ_l是已选原子集合，则子空间追踪算法步骤如下：

输入：K，D，X，l，Λ_l，Ψ_l，R^l；

初始化设置：l＝1，R⁰＝X，

迭代：

1)更新候选原子集合Λ_i：

Λ_l＝Λ_l-1∪{使得表达式|<R^l，D>|值最大的K个原子}；

2)计算投影系数u_l：

3)更新已选原子集合Ψ_l：

Ψ_l＝{使得u_l最大的K个原子}；

4)更新重构信号Y_l与残差R^l：

R^l＝X-Y_l； (21)

5)如果||R^l||₂＞||R^l-1||₂，则Ψ_l＝Ψ_l-1且停止迭代；否则，l＝l+1，并返回步骤1)继续迭代；

输出：重构信号Y_c＝Y_l；

重构信号Y_c即为稀疏分解提取到的人文噪声，含噪的高频信号X减去人文噪声Y_c即得到去噪后的高频有效信号Y_h：

Y_h＝X-Y_c； (22)

将低频有效信号C与高频有效信号Y_h相加，即得到去噪后的全频带的大地电磁有效信号Y：

Y＝C+Y_h。 (23)。