CN110147766B - 基于移不变稀疏编码的低频大地电磁信号去噪方法 - Google Patents

Abstract

本发明提供一种基于移不变稀疏编码的低频大地电磁信号去噪方法。所提方法包括：首先利用数学形态滤波将含噪的大地电磁时间序列分解成低频有效信号和含噪的高频信号；然后将含噪的高频信号分成若干片段，利用移不变稀疏编码从含噪的高频信号中自主学习到人文噪声的特征结构构建学习型冗余字典，并利用学习到的冗余字典对含噪的高频有效信号进行信噪分离，得到去噪后的高频有效信号。低频有效信号与高频有效信号相加即得到全频带的大地电磁有效信号。本发明可以有效去除低频大地电磁数据中的强人文噪声，显著提高大地电磁数据质量，改善大地电磁法深部探测效果。

Description

基于移不变稀疏编码的低频大地电磁信号去噪方法

技术领域

本发明属于勘探地球物理信号处理领域，涉及一种基于移不变稀疏编码机器学习算法的大地电磁信噪分离方法，尤其是针对10^-4～10¹Hz范围内的低频大地电磁数据的处理方法。

背景技术

大地电磁法具有探测深度大、无需人工场源等优点，被广泛应用于矿产资源勘查以及深部电性结构探测。但天然的大地电磁信号随机性强，幅度微弱，且频率范围很宽(10^-4～10⁵Hz)，极易受人文噪声的污染。

数十年来，国内外学者提出了许多方法用于压制或者分离大地电磁信号中的噪声。其中远参考法和Robust法获得广泛认可，已成为通用方法。然而随着工业文明覆盖范围的增大，远参考点的选取变得愈发困难；Robust法要求观测到的数据大部分属于高质量的数据，因此对于持续性噪声效果不佳；数据段挑选法也是一种有效的方法(Weckmann etal.,2005)，但同样仅适用于非持续性噪声。为应对持续性人文噪声，小波滤波、形态滤波、经验模态分解(Empirical Mode Decomposition，EMD)等手段被用于MT数据处理，这些方法在10Hz以上中高频段大地电磁数据(即音频大地电磁信号)的处理中取得了较好的效果，但由于存在损失低频有效信号的问题，难以用于10Hz以下低频数据的处理。此外，EMD算法存在严重的模态混淆问题，导致不同频率成分的信号存在于同一个模态分量中，或者同一频率成分的信号分布在多个模态分量，无法准确区分高频与低频有效信号。

专利CN201610618524.6将经验模态分解与数学形态滤波结合，先利用经验模态分解将含噪信号分解成高阶模态分量(低频部分)与低阶模态分量(高频部分)；然后利用平滑滤波法对高阶模态分量进行去噪，消除基线漂移，得到去噪后的低频信号；利用形态滤波与自适应阈值去噪方法对低阶模态分量进行去噪，得到去噪后的高频信号；最后将去噪后的高频信号与低频信号合并，得到去噪后的信号。然而利用平滑滤波法对高阶模态分量(低频部分)进行去噪时，极易损失10Hz以下的低频信号，因此该方法仍然主要用于10Hz以上的高频大地电磁信号的处理。

为了更好地保留低频有效信号，汤井田和李广等(Li et al,2017,AppliedGeophysics；汤井田等,2017,地球物理学报)将稀疏表示引入大地电磁信号处理，通过设计与人文噪声匹配但对有效信号不敏感的冗余字典，实现仅对含噪片段进行处理；李晋等(CN201710938046.1、CN201810858136.4、CN201810858135.X)先将原始信号辨识为高质量片段和强干扰段，然后仅对强干扰段进行去噪。这两类方法能够有效的保护高质量片段中的低频信号，但是对于含噪片段的处理仍然是直接提取人文噪声，损失了强干扰段中的低频有效信号。

总的来说，由于大地电磁信号的频率很宽(10^-4～10⁵Hz)，现有的大地电磁数据处理方法难以满足各个频段数据的处理，尤其是对10Hz以下低频大地电磁数据处理时，现有的信噪分离方法有明显的局限性。对于1Hz以下低频大地电磁数据的处理，仍然只能采用通用的参考道法或者Robust法处理，鲜有专门针对低频大地电磁数据的信噪分离方法。

发明内容

为解决上述问题，本发明公开了一种专门针对10Hz以下的基于移不变稀疏编码的低频大地电磁信号去噪方法，创造性地利用数学形态滤波技术和移不变稀疏编码技术相结合，先通过数学形态滤波技术手段从含噪信号中提取出低频有效信号，然后再用移不变稀疏编码对剩余信号去噪。

数学形态学源于Matheron和Serra对二元多孔介质的几何形态的研究。数学形态滤波是由数学形态学发展而来的一种非线性滤波技术，被广泛用于图像与信号分析处理。它包括膨胀和腐蚀两个基本运算。腐蚀运算是收缩的过程，用来消除物体的边界点，减少峰值。膨胀运算是扩张的过程，用来填平边界不平滑的凹陷部分，增大谷值。当形态滤波结构元素尺寸大于人文噪声所占的采样点时，可以用形态滤波提取大地电磁信号中的低频缓变部分。在大地电磁信号处理领域，形态滤波通常被用于提取噪声(例如论文Li et al.,2017,Applied Geophysics)，本发明则相反，利用形态滤波提取大地电磁低频有效信号。

移不变稀疏编码(SISC)是一种基于数据驱动的机器学习算法。它引入卷积作为移位算子，使得特征原子可以自由的平移、翻转、缩放，从而通过一个特征原子方便的表示时间序列中位置不同但特征相同的信号成分。它也是一种字典学习方法，从给定样本中自主学习到同一来源信号的特征结构，然后每一个特征结构重构出一个独立的信号分量。SISC适用于各类成分复杂信号的分解，因此特别适合来源复杂的大地电磁时间序列的处理。

发明人通过不断地探索，将数学形态滤波技术和移不变稀疏编码技术相结合，用于大地电磁数据处理中，特别是针对10^-4～10¹Hz的低频大地电磁信号，取得了出人意料的效果，显著提高数据质量。

本发明技术方案，基于移不变稀疏编码的低频大地电磁信号去噪方法，包括以下步骤：

步骤1、利用数学形态滤波(Mathematical Morphological Filtering,MMF)将含噪的大地电磁时间序列A分解成低频有效信号B和含噪的高频信号X，即A＝B+X；

步骤2、利用移不变稀疏编码(Shift-Invariant Sparse coding,SISC)将X分解成高频有效信号Y_h和人文噪声Y_c，高频有效信号Y_h和低频有效信号B相加得到全频带的大地电磁有效信号Y。

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

设A＝a(c)＝[a₁,a₂,…,a_C]^T是长度为C的待处理含噪信号，g＝g(j)＝[g₁,g₂,…,g_J]^T是长度为J的结构元素，处理电通道信号时，结构元素的长度J取值为150；处理磁通道信号时，结构元素的长度J取值为100，根据实际情况，J取值可适当调整。用数学形态滤波提取到的低频有效信号B为：

B＝Φ_OC(CO)(g)＝[OC(A)+CO(A)]/2， (1)

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

式中

表示开运算，“·”表示闭运算，开运算和闭运算定义为:

式中“Θ”和

分别表示腐蚀和膨胀运算：

进一步的，所述步骤2使用移不变稀疏编码处理前，先将X分成K段，每段所占采样点数为27000(采样率为15Hz时，每段时长为30分钟)，即X＝[x₁,x₂,…,x_K]^T，x_k＝[x₁,x₂,…,x_N]^T，N＝27000。SISC使用一个包含M个基(或称为原子)的过完备字典D＝[d₁,d₂,…,d_M]∈R^{Q ×M}，将x_k表示为基与系数的卷积之和：

式中“*”表示卷积，编码系数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稀疏程度的代价函数为：

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

移不变稀疏编基本码流程如下：

输入：待处理信号X，特征原子的种类α，特征原子的宽度W＝[W₁,W₂,…W_a]，特征原子总数(即稀疏度)L，最大迭代次数z(本发明取值为15)，收敛目标e<10^-6；

初始化：对D、s随机赋予初值；

循环执行z次以下内容：

当达到最大循环次数或达到收敛目标时停止循环；

输出：过完备字典D以及重构信号

即为提取到的人文噪声(记为Y_c)，含噪的高频信号X减去人文噪声Y_c即得到去噪后的高频有效信号Y_h:

Y_h＝X-Y_c。 (10)

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

Y＝B+Y_h。 (11)

进一步的，设g_i,u为特征原子d_i平移u点补零后的原子，其长度与x_k相同，且||g_i,u||＝1，最大循环次数L_max，第L次循环后的残差为r^L，第L次循环后已选原子集合ψ_L，则系数求解步骤如下：

1)参数初始化，r⁰＝x_k，

L＝1；

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

3)更新已选原子集合，

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

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

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

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

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

将上述公式左边矩阵视为系数s_i,k的一种特殊的Toeplitz矩阵，则上式可以写作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个信号时，优化函数ψ(θ)可以表示为：

将上式简写为S·d_i＝E，根据最小二乘法可得d_i＝(S^TS)^-1(S^TE),矩阵(S^TS)∈R^Q×Q，一般情况下Q≤N，因此最后就转化为求解一个小规模的线性方程组的解。这一求解过程通过Cholesky分解直接求取最优值。按照随机的顺序依次对每个基进行更新，能够求取基在本轮循环中的最优解，然后依次对基d_i标准化处理

即可得到完整的字典D。

本发明具有以下优点：

1、本发明提供了一种针对10^-4～10¹Hz范围内的低频大地电磁数据的信噪分离方法，与现有技术易损失低频有效信号相比，在去噪前增加了低频有效信号提取的环节，有效避免了10Hz以下低频有效信号的损失。不仅能够保留高质量片段中的低频信号，还能有效保留强干扰片段中的低频有效信号，并准确去除强人文噪声，在对10Hz以下低频大地电磁数据去噪时，相对于直接提取人文噪声的方法，取得明显更好的效果。

2、本发明基于数学形态滤波的低频缓变部分提取，它能够自适应的从大地电磁时间序列中提取到代表深部信息的低频缓变部分，确保在后续的信噪分离中不会损失低频有效信号。该方法不仅可以用于大地电磁信号，也可以用于地震等其他信号的低频缓变部分的提取。

3、本发明基于移不变稀疏编码机器学习的人文噪声提取，它能够从含噪信号中自主学习到人文噪声的特征结构，构建学习型冗余字典。与解析型冗余字典仅适用于某一类或者某几类人文噪声相比，学习型字典对于各类形态的人文噪声均具有良好的相效果，不受人文噪声种类与形态的限制，灵活性大大提高。

附图说明

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

图2为合成数据处理实验结果对比图；

图3为青海实测时间序列去噪前后处理结果对比图；

图4为安徽铜陵矿集区实测时间序列去噪前后处理结果对比图；

图5为青海某实测点处理前后视电阻率-相位曲线对比图。

具体实施方式

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

附图1为本发明的基本流程，包括如下两个步骤：

步骤1，利用数学形态滤波将含噪的大地电磁时间序列A分解成低频有效信号B和含噪的高频信号X；

步骤2，对X分段，并利用移不变稀疏编码将X分解成高频有效信号Y_h和人文噪声Y_c，高频有效信号Y_h和低频有效信号B相加得到全频带的大地电磁有效信号Y。

实施例1

将本发明所述信噪分离方法应用于合成数据的处理，得到如图2所示实验结果。图2中(a)是青海柴达木盆地采集的时间序列片段，采样频率为15Hz，显然，该时间序列中没有明显的强人文噪声，信号变化平缓，属于高质量数据；图2中(b)是在(a)中添加仿真的强人文噪声后得到的含噪信号，图2中(c)是直接使用数学形态滤波法去噪后的结果，不难发现形态滤波法在去噪的同时，去除了几乎全部的低频缓变部分，损失了大量的低频有效信号，且残留有少量的脉冲噪声；图2中(d)是使用本发明所述的方法去噪后得到的信号，显然，含噪信号中的强人文噪声被准确去除，且没有明显的损失低频有效信号，处理后结果与原始高质量信号具有很高的相似度。

实施例2

与实施例1类似，将本发明所述信噪分离方法应用于青海实测时间序列的处理，得到如图3所示实验结果。不难发现，去噪前图3中(a)，时间序列中含有大量的脉冲噪声，且其幅度明显大于有效信号。去噪后图3中(b)，脉冲噪声被压制，且低频缓变信号得到很好的保留。

实施例3

与实施例1以及实施例2类似，将本发明所述信噪分离方法应用于安徽铜陵矿集区实测大地电磁数据。得到如图4所示实验结果。如图4中(a)所示，处理前时间序列中含有大量的脉冲和方波噪声，且其幅度明显大于有效信号。如图4中(b)所示，使用本发明去噪后，方波噪声和脉冲噪声被准确去除。同样的，低频信号未有明显的损失。

通过对比使用本发明所述信噪分离方法处理前与处理后的视电阻率-相位曲线，来说明本发明的有益效果。如图5所示，为青海柴达木盆地某实测点处理前后视电阻率-相位曲线。如图5中(a1,a2)，处理前，10Hz到0.01Hz范围内的频点跳变严重，尤其是1Hz附近，受到了严重的干扰，0.01Hz以下也有个别频点畸变严重，且误差棒较大。如图5中(b1,b2)，经过本发明所提方法处理后，10Hz到0.01Hz范围内曲线连续性与光滑度得到显著改善，误差棒减小，0.01Hz以下畸变的频点也有所改善。总体说来，视电阻率-相位曲线的质量得到显著提高。

综上所述，本发明可以在不损失低频有效信号的前提下，有效去除10^-4～10¹Hz范围内的低频大地电磁数据中的强人文噪声，显著提高数据质量，获得更加可靠的深部信息，改善大地电磁法在恶劣电磁环境下的应用效果。

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

Claims (1)

1.一种基于移不变稀疏编码的低频大地电磁Magnetotelluric,MT信号去噪方法，包括以下步骤：

步骤1、利用数学形态滤波Mathematical Morphology Filtering,MMF将含噪的大地电磁时间序列A分解成低频有效信号B和含噪的高频信号X，即A＝B+X；其操作如下：

设A＝a(c)＝[a₁,a₂,…,a_C]^T是长度为C的待处理含噪信号，g＝g(j)＝[g₁,g₂,…,g_J]^T是长度为J的结构元素，处理电通道信号时，结构元素的长度J取值为150；处理磁通道信号时，结构元素的长度J取值为100，根据实际情况，用数学形态滤波提取到的低频有效信号B为：

B＝Φ_OC(CO)(g)＝[OC(A)+CO(A)]/2， (1)

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

式中

表示开运算，“·”表示闭运算，开运算和闭运算定义为:

式中“Θ”和

分别表示腐蚀和膨胀运算：

步骤2、利用移不变稀疏编码Shift-Invariant Sparse Coding,SISC将X分解成高频有效信号Y_h和人文噪声Y_c，高频有效信号Y_h和低频有效信号B相加得到全频带的大地电磁有效信号Y；

使用移不变稀疏编码处理前，先将X分成K段，在采样率为15Hz，每段时长为30分钟的情况下，每段所占采样点数为27000，即X＝[x₁,x₂,…,x_K]^T，x_k＝[x₁,x₂,…,x_N]^T，N＝27000，SISC将x_k表示为基与系数的卷积之和：

式中“*”表示卷积，编码系数向量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稀疏程度的代价函数为：

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

移不变稀疏编码基本流程如下：

输入：待处理信号X，特征原子的种类α，特征原子的宽度W＝[W₁，W₂，…W_α]，特征原子总数L即稀疏度，最大迭代次数z，收敛目标e<10^-6；

初始化：对D、s随机赋予初值；

循环执行z次以下内容：

当达到最大循环次数或达到收敛目标时停止循环；

输出：过完备字典D以及重构信号

以处理强干扰环境下MT数据为例，α可设置为10，原子宽度可设置为W＝[1,3,10,30,30,50,50,100,100,150]，

即为提取到的人文噪声Y_c，含噪的高频信号X减去人文噪声Y_c即得到去噪后的高频有效信号Y_h:

Y_h＝X-Y_c， (10)

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

Y＝B+Y_h； (11)

设g_i,u为特征原子d_i平移u点补零后的原子，其长度与x_k相同，且||g_i,u||＝1，最大循环次数L_max，第L次循环后的残差为r^L，第L次循环后已选原子集合ψ_L，则系数求解步骤如下：

1)参数初始化，r⁰＝x_k，

L＝1；

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

3)更新已选原子集合，

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

残差r^L＝x_k-S_Lψ_L，重构信号

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

字典更新通过以下步骤实现：保持系数不变，仅对基进行更新，代价函数可以简化为：

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

将上式左边矩阵视为系数S_i，k的一种特殊的Toeplitz矩阵，则上式可以写作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个信号时，优化函数ψ(θ)可以表示为：

将上式简写为S·d_i＝E，根据最小二乘法可得d_i＝(S^TS)^-1(S^TE),矩阵(S^TS)∈R^Q×Q，Q≤N，因此最后就转化为求解一个小规模的线性方程组的解；这一求解过程通过Cholesky分解直接求取最优值，按照随机的顺序依次对每个基进行更新，能够求取基在本轮循环中的最优解，然后依次对基d_i标准化处理

即可得到完整的字典D。

Images