CN110850482A - 一种基于变分模态分解原理的瞬变电磁信噪分离方法 - Google Patents

Abstract

本发明涉及地球物理勘探中瞬变电磁信号处理技术领域，公开了一种基于变分模态分解原理的瞬变电磁信噪分离方法，适用于时间域的电磁数据滤波，该方法包括：将含噪的瞬变电磁数据s(t)分解为K个本征模态函数u_k(t)；对每一个本征模态函数进行希尔伯特变换，得到单边频谱；根据单边频谱将各本征模态函数u_k(t)的频谱解调到相对应的基频带，计算解调后的L2范数，并计算分解出的本征模态函数的带宽，得到约束变分模型；引入二次惩罚因子α和拉格朗日乘法算子λ，将约束变分模型转换为非约束性变分形式；计算非约束性变分形式的最优解，得到由VMD分解后的一系列优化的本征模态分量。解决地质解释中采集的瞬变电磁数据被噪声干扰的问题，该方法简便且具有优异的噪声分离、幅值保持性能，还可以将噪声图像提取出来，为后续对噪声的分析提供条件。

Description

一种基于变分模态分解原理的瞬变电磁信噪分离方法

技术领域

本发明涉及地球物理勘探中瞬变电磁信号处理技术领域，公开了一种基于变分模态分解原理的瞬变电磁信噪分离方法，适用于时间域的电磁数据滤波。

背景技术

瞬变电磁法是一种在电磁法勘探领域中广泛使用的方法，它使用一个发射机发射电流，断开电流时产生一个时变场，在地下传播过程中产生二次场，通过接收机接收到的二次场来判断地下情况。然而，接收到的信号较弱且易受工频，天电，人文及地质干扰影响。因此，对采集到的信号进行消噪是后续的成像及地质解释的前提。

目前，关于瞬变电磁信号的消噪已经提出了很多种方法，主要有自适应滤波法、双极性采样法、近似抽道方法、小波分析法和经验模态分解方法等。其中小波方法效果很好，但是需要根据瞬变电磁信号选择合适的小波基，不够智能。2017年提出经验模态分解(Empirical Mode Decomposition，EMD)方法应用于将信号分解为未知但不同频率带的不同模式,是信号分解领域一个热门方法，广泛用于处理非平稳、非线性的信号。经验模态分解方法本质是将需要的频率成分筛选出来，存在模态混叠现象，且耗时长，在处理多尺度的数据时存在一定限制。

发明内容

本发明所要解决的技术问题在于提供一种基于变分模态分解原理的瞬变电磁信噪分离方法，解决地质解释中采集的瞬变电磁数据被噪声干扰的问题，该方法简便且具有优异的噪声分离、幅值保持性能，还可以将噪声图像提取出来，为后续对噪声的分析提供条件。

本发明是这样实现的，一种基于变分模态分解原理的瞬变电磁信噪分离方法，包括:

步骤1：将含噪的瞬变电磁数据s(t)分解为K个本征模态函数u_k(t)；

步骤2：对每一个本征模态函数进行希尔伯特变换，得到单边频谱；

步骤3：根据步骤2得到的单边频谱将各本征模态函数u_k(t)的频谱解调到相对应的基频带，计算解调后的L2范数，并计算分解出的本征模态函数的带宽，得到约束变分模型；

步骤4：引入二次惩罚因子α和拉格朗日乘法算子λ，将步骤3得到的约束变分模型转换为非约束性变分形式；

步骤5：计算步骤4得到的非约束性变分形式的最优解，得到由VMD分解后的一系列优化的本征模态分量。

进一步地，

步骤4中，根据瞬变电磁信号e指数衰减的特征，将早期与晚期的α分段。

进一步地，

步骤1中，将瞬变电磁实测数据s(t)分解为K个本征模态函数；将各个本征模态函数记为：

式中，A_k(t)是u_k(t)的瞬时幅值，为u_k(t)的相位，记t时刻的瞬时频率为ω_k(t),

在邻域范围内，u_k(t)为一个谐波信号。

进一步地，

步骤5中包含以下步骤：

I：选取VMD分解个数K值与数值小的惩罚因子α处理早期信号；

II：初始化

其中

为初始循环中第K个本证模态函数，

为初始循环中的中心频率，

为初始循环中的拉格朗日算子，设置n＝n+1，k＝k+1，交替更新{u_k(t)},{ω_k(t)}直到

k＝K结束循环:

其中，ω_k(t)为中心频率，u_k(t)为模态分量，s(ω)为s(ω)为含噪瞬变电磁信号的傅里叶变换，

是在第n+1次循环时第k个模态函数，

是在第n+1次循环时第k个中心频率；

III：重复步骤II，直到满足判定精度ε如式(3)所示，结束循环得到K个由变分模态分解后的本征模态分量；

IV：重新选取α值处理晚期信号，重复步骤I至III，得到完整的消噪后的瞬变电磁信号。

本发明与现有技术相比，有益效果在于：

本发明提出基于变分模态分解(Variational Modal Decomposition,VMD)的瞬变电磁信号噪声压制处理方法。变分模态分解的过程是求解变分模型的最优解，将每个分量的估计带宽之和最小作为约束，用交替方向乘子法不停更新各个模态及其中心频率，将每个分量的带宽和中心频率在频域内进行自适应分离，从而适当地平衡它们之间的误差，实现瞬变电磁信号与噪声分离。根据瞬变电磁信号本身的特点来选取VMD算法的参数。

附图说明

图1是本发明基于变分模态分解原理的瞬变电磁信噪分离方法的流程图；

图2是仿真的不含噪的瞬变电磁信号(a)及其频谱(b)；

图3是加入高斯噪声后的瞬变电磁信号(a)及其频谱(b)；

图4是对高斯噪声去噪结果对比图；

图5是加入工频噪声后的瞬变电磁信号(a)及其频谱(b)；

图6是对工频噪声去噪结果对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参见图1所示，一种基于变分模态分解原理的瞬变电磁信噪分离方法，将地空电磁数据分解得到K个模态函数，将每个分量的估计带宽之和最小作为约束，用交替方向乘子法不停更新各个模态及其中心频率，进而将每个分量的带宽和中心频率在频域内进行自适应分离，从而适当地平衡它们之间的误差，根据地空数据呈现e指数衰减的特点，在选取α的时候进行分段，早期信号幅值大，噪声影响较小，选择数值较小的α可以达到消噪的目的，并且幅值下降的影响不大，晚期信号幅值较小，噪声影响大，选择数值较大的α才能取的良好的消噪效果。

本发明实施例中方法包括：

1)将瞬变电磁仿真数据如图2(a)所示加入10％高斯白噪声得到图3(a)的含噪瞬变电磁数据s(t)，图2(b)为与图2(a)对应的瞬变电磁仿真数据的频谱，图3(b)为与图3(a)对应的含高斯白噪声的瞬变电磁信号频谱。s(t)分解为K个本征模态函数，K可以根据当地主要包含噪声情况选择；将各个本征模态函数记为：

式中，A_k(t)是u_k(t)的瞬时幅值，

为u_k(t)的相位，记t时刻的瞬时频率为ω_k(t),

可以在邻域范围内，认为u_k(t)是一个谐波信号。

2)对每一个本征模态函数即IMF分量进行希尔伯特变换，得到单边频谱如下：

式中，u_k(t)是第k个IMF分量，u_k(t)可看作谐波信号；δ(t)为冲激函数，*为卷积运算；

3)将步骤2)各本征模态函数u_k(t)的频谱解调到相对应的基频带，即：

式中{ω_k}＝{ω₁，…，ω_k}为各u_k(t)的中心频率。

4)解调后计算其L2范数，并用范数表达式估计分解出的IMF分量带宽，各IMF分量对应的约束变分模型，表达式如下：

5)为了将上述约束性变分问题转换为非约束性变分问题，通过引入二次惩罚因子α和拉格朗日惩罚算子λ，α的选择主要根据瞬变电磁信号频谱确定，因为瞬变电磁信号频谱范围很广，而中心频率幅值较周围大很多，因此根据瞬变电磁信号的特点，将地空电磁数据早期的α选用较小的数值，晚期选用较大的数值可以取得良好的效果。用增广lagrange函数来计算其最优解。其中，二次惩罚因子可保证信号重构精度；Lagrange乘子λ(x)可以保证模型约束条件的严格性。通过采用交替方向乘子算法,实现的交替更新，模态分量{u_k(t)}与相应的中心频率{ω_k(t)}如下所示：

其中，s(ω)为含噪瞬变电磁信号的傅里叶变换，

是在第n+1次循环时第k个模态函数，

是在第n+1次循环时第k个中心频率。

初始化

其中

为初始循环中第K个本证模态函数，为初始循环中的中心频率，

为初始循环中的拉格朗日算子设置n＝n+1，k＝k+1，根据公式(6)、(7)更新{u_k(t)},{ω_k(t)}直到k＝K结束循环。对于所有ω_k＞0，更新λ，即：

其中，λⁿ⁺¹(ω)表示第n+1次循环中的拉格朗日算子。

式中γ代表噪声容限参数,给定判定精度ε，直到满足迭代停止条件：

即得到由变分模态分解后的一系列优化的本证模态分量。

重新选取α值处理晚期信号，重复上述步骤，得到完整的消噪后的瞬变电磁信号。消噪效果对比图如图4所示。

另一实施例中与上述的实施例不同之处在于，对将瞬变电磁仿真数据如图2(a)所示加入工频谐波噪声得到图5(a)的含噪瞬变电磁数据，图5(b)为与图5(a)对应的含工频谐波噪声的瞬变电磁信号频谱。得到消噪后的对比图如图6所示。实现瞬变电磁信号与噪声分离。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于变分模态分解原理的瞬变电磁信噪分离方法，其特征在于，该方法包括：

步骤1：将含噪的瞬变电磁数据s(t)分解为K个本征模态函数u_k(t)；

步骤2：对每一个本征模态函数进行希尔伯特变换，得到单边频谱；

步骤3：根据步骤2得到的单边频谱将各本征模态函数u_k(t)的频谱解调到相对应的基频带，计算解调后的L2范数，并计算分解出的本征模态函数的带宽，得到约束变分模型；

步骤4：引入二次惩罚因子α和拉格朗日乘法算子λ，将步骤3得到的约束变分模型转换为非约束性变分形式；

步骤5：计算步骤4得到的非约束性变分形式的最优解，得到由VMD分解后的一系列优化的本征模态分量。

2.按照权利要求1所述的方法，其特征在于，步骤4中，根据瞬变电磁信号e指数衰减的特征，将早期与晚期的α分段。

3.按照权利要求1所述的方法，其特征在于，

步骤1中，将瞬变电磁实测数据s(t)分解为K个本征模态函数；将各个本征模态函数记为：

式中，A_k(t)是u_k(t)的瞬时幅值，

为u_k(t)的相位，记t时刻的瞬时频率为ω_k(t),

在邻域范围内，u_k(t)为一个谐波信号。

4.按照权利要求3所述的方法，其特征在于，

步骤5中包含以下步骤：

I：选取VMD分解个数K值与数值小的惩罚因子α处理早期信号；

II：初始化

其中

为初始循环中第K个本证模态函数，

为初始循环中的中心频率，

为初始循环中的拉格朗日算子，设置n＝n+1，k＝k+1，交替更新{u_k(t)},{ω_k(t)}直到

k＝K结束循环:

其中，ω_k(t)为中心频率，u_k(t)为模态分量，s(ω)为s(ω)为含噪瞬变电磁信号的傅里叶变换，

是在第n+1次循环时第k个模态函数，

是在第n+1次循环时第k个中心频率；

III：重复步骤II，直到满足判定精度ε如式(3)所示，结束循环得到K个由变分模态分解后的本征模态分量；

IV：重新选取α值处理晚期信号，重复步骤I至III，得到完整的消噪后的瞬变电磁信号。

Images