CN114611329B - 一种基于变分模态分解的时域电磁法近场噪声压制方法 - Google Patents

Abstract

本发明属于地球物理勘探技术领域，具体公开了一种基于变分模态分解(Variational ModeDecomposition，VMD)的时域电磁法近场噪声压制方法，包括设置固定基站和至少一移动站同步采集信号，利用固定基站和移动采集站之间噪声的相关性来进行信号压制，①得到基站信号，选取多个周期并使用VMD算法处理后重构得到去噪的基信号；②将信号减去基信号得到噪声；③将移动采集站信号减去噪声得到LNC压制后的信号。本方法根据变分最优解方法对基站信号处理得到噪声信号，再利用噪声进行近场人文噪声压制(Local Noise Compensation，LNC)，进而获得更高信噪比的可控源电磁信号。

Description

一种基于变分模态分解的时域电磁法近场噪声压制方法

技术领域

本发明属于地球物理勘探技术领域，具体涉及一种基于变分模态分解的时域电磁法近场噪声压制方法。

背景技术

电磁勘探发展至今已经衍生出许多方法，主要分为天然源电磁法和可控源电磁法，不同电磁法对于放大有用信号和去除其他噪声都有着特色的方法。瞬变电磁法是典型的可控源电磁法，其由发射系统和采集系统构成。TEM在地下矿物资源和地下结构变化监测，地层识别等各方面已经得到广泛应用，可控源电磁法是目前应用较广泛的一种电磁方法，野外采集的可控源电磁信号通常包含真实信号、人文噪声及背景噪声。传统的噪声压制方法是通过设计滤波器滤除周期性噪声，通过叠加压制随机噪声(即叠后处理)。但是在强干扰地区就需要采集较多的周期来满足叠加需求，这导致工作量加大的同时也带来频谱和信号畸变的问题，这些问题都增大了去噪工作的难度，最终不能有效得到真实信号。

发明内容

本发明意在提供一种基于变分模态分解的时域电磁法近场噪声压制方法，根据变分最优解方法对基站信号处理得到噪声信号，在利用得到的噪声进行近场人文噪声压制，进而获得更高信噪比的可控源电磁信号。

为了达到上述目的，本发明的技术方案如下：

一种基于变分模态分解的时域电磁法近场噪声压制方法，包括，

设置一固定基站和至少一个移动站同步采集信号，利用固定基站和移动采集站之间噪声的相关性来进行噪声压制，具体如下：

①得到基站信号b(t)，选取多个周期并使用VMD算法处理后重构得到去噪的基信号S；

②将信号b(t)减去基信号S得到噪声n(t)：

n(t)＝b(t)-S；

③将移动采集站的信号m(t)减去噪声得到LNC压制后的信号c(t)：

c(t)＝m(t)-n(t)。

进一步的，所述周期的筛选标准如下：

对各个周期与基站所有周期信号的均值做相关性计算，挑选相关系数绝对值大于0.95的周期，若无周期满足，则挑选相关系数最大的那部分周期。

进一步的，还包括对个周期的LNC压制后的信号c(t)进行指数加权滑动平均的步骤。

进一步的，固定基站和移动采集站之间使用性能和型号一致的仪器采集数据。

进一步的，步骤①中挑选基站信号的n个周期进行叠加并取平均值，得到基叠加信号：

用于VMD算法处理。

进一步的，步骤①的VMD算法，根据以下约束变分模型和增广拉格朗日函数求解各阶窄带分量IMF：

其中f是输入的原信号，u_k是分解得到的K阶IMF，ω_k是各IMF的中心频率，α是二次惩罚因子，λ(t)是拉格朗日乘子。

进一步的，步骤①的VMD算法中，在获得的VMD分解各阶IMF_i后，使用傅立叶变换计算原信号的主频和IMF_i主频，根据频谱结果去除特征明显的噪声IMF_i去除，最后把剩余的IMF_i重构为去噪后的基叠加信号/>如下：

其中m为剔除噪声对应的IMF_i后余下的IMF_i的数量。

进一步的，步骤②中，在得到去噪的基叠加信号S_stand(t)后，将各周期的基站信号S_base(t)减去去噪的基叠加信号S_stand(t)得到各周期的噪声信号N_stand(t)如下：

{N_stand(t)}＝{S_base(t)}-S_stand(t)。

进一步的，步骤③中，将各周期的移动站信号S_move(t)减去各周期的噪声信号N_stand(t)得到去噪后的各周期的移动站信号newS_move(t)如下：

{newS_move(t)}＝{S_move(t)}-{N_stand(t)}。

进一步的，按以下方式对去噪后的各周期的移动站信号newS_move(t)进行指数加权滑动平均得到最后的信号结果S_vmd-lnc(与)：

其中β为权重。

进一步的，β取0.9。

本发明的有效果在于，使用本发明中的噪声压制方法后得到的结果与其他方法相比信噪比得到大幅提高，而且方法简便易与实现，在压制强干扰方面有着优秀的表现，通过变分模态分解(Variational Mode Decomposition，VMD)后再重构信号本身可以很好的压制高斯白噪声，通过剔除噪声对应的IMF，本方法还可以针对性压制人文噪声，这使得通过本方法中的重构过程得到的信号十分接近发射信号，从接收信号中减去该信号就可以比较准确的得到“噪声信号”，于是可以利用该“噪声信号”对移动站的接收信号中的真实噪声进行压制。相比于现有技术中的叠加后处理技术，本方法根据噪声本身特征进行了针对性压制，由于噪声压制机制的存在，信号采集的时候不再需要采集大量周期的信号，缩减了采集时间，节省了采集成本。

附图说明

图1为本发明实施例中的基于变分模态分解的时域电磁法近场噪声压制方法的流程图。

图2为本发明实施例中的仿真实验所用到的场景设置示意图。

图3为本发明实施例中的仿真实验所用的可控源电磁电场曲线示意图。

图4为本实施例中的添加噪声前后的基站及移动站的电场曲线对比图。

图5为从本发明实施例中添加了噪声的基站信号的时间序列周期中随机抽取四个周期组成的展示图。

图6为从本发明实施例中添加了噪声的移动站信号时间序列周期中随机抽取四个周期组成的展示图。

图7为本发明实施例的模拟实验中得到的IMF₁的频谱图。

图8为本发明实施例的模拟实验中得到的IMF₂的频谱图。

图9为本发明实施例的模拟实验中最终VMD-LNC-EWMA处理得到的信号结果与原信号以及添加了噪声的信号的曲线对比图。

图10为以本发明实施例中的模拟实验为基础仿真得到的最终信号的信噪比与VMD参数之间的关系示意图。

图11为本发明实施例中实测采集的部分电性源瞬变电磁信号的示意图。

图12为本发明实施例中实测采集的移动站信号的频谱图。

图13为本发明实施例中从实测采集的数据中分解得到的IMF₁的频谱图。

图14为本发明实施例中从实测采集的数据中分解得到的IMF₂的频谱图。

图15为本发明实施例中对实测采集的数据进行VMD-LNC-EWMA处理后得到的结果的频谱图。

图16为本发明实施例中对实测采集的移动站信号原始数据以及进行常规处理和VMD-LNC-EWMA处理得到的最终结果间的对比图。

具体实施方式

下面将结合本发明实施例的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例中的基于变分模态分解的时域电磁法近场噪声压制方法基本如图1所示，包括：

设置一固定基站和至少一个移动站同步采集信号，利用固定基站和移动采集站之间噪声的相关性来进行噪声压制，具体如下：

①得到基站信号b(t)，选取多个周期并使用VMD算法处理后重构得到去噪的基信号S；

②将信号b(t)减去基信号S得到噪声n(t)：

n(t)＝b(t)-S； (1)

③将移动采集站的信号m(t)减去噪声得到LNC压制后的信号c(t)：

c(t)＝m(t)-n(t)。 (2)

其中，VMD算法由EMD算法改进而来，但VMD是在变分这一数学框架下实现的，这使得其分解结果有依据且更加准确。VMD算法和EMD算法分解得到的模态都被称为本征模态函数(IMF)，不同于EMD中的定义，VMD算法中的IMF被定义为一个调幅-调频信号：

上式中A_k(t)为u_k(t)的瞬时幅值。瞬时频率ω_k(t)是对时间的导数，即：

VMD算法求解IMF是通过搜索增广拉格朗日函数最优解来进行的，人为设置输出IMF的数量K后通过拉格朗日乘数法计算每个IMF的中心频率和频带并依次迭代，最终将K个窄带IMF输出。其约束变分模型和求解引入的增广拉格朗日函数如下(式(5)、(6))：

上面二式中f是输入的原信号，u_k是分解得到的K阶IMF，ω_k是各IMF的中心频率，α是二次惩罚因子，λ(t)是拉格朗日乘子。

于是，本实施例中，首先将挑选出的基站信号的n个周期进行叠加并取平均值，得到基叠加信号：

所需周期的筛选标准是，对各个周期与基站所有周期信号的均值做相关性计算，挑选相关系数绝对值大于0.95的周期，若无周期满足，则挑选相关系数最大的那部分周期。而相关系数计算公式如下：

式中Cov(x，y)为x，y的协方差，Var[x]、Var[y]为x，y的方差。

下一步依照式(5)、(6)求解的VMD分解各阶IMF_i，需要选取的参数包括二次惩罚因子α、起始频率ω₀和分解次数K。

接下来需要将各阶IMF_i信号进行筛选以用于重构，对IMF_i的筛选依据频谱分析，使用傅立叶变换计算原信号的主频和IMF主频，根据频谱结果将特征明显的噪声IMF去除(如电网附近的50Hz工频干扰、矿场采矿机造成的类脉冲干扰、自然界随机白噪声等)，最后把剩余的IMF重构为去噪后的基叠加信号，另一方面由于VMD分解本身自带抗噪性，在重构的信号中高斯白噪声也已被压制。

在得到去噪的基叠加信号后，将各周期的基站信号减去去噪的基叠加信号得到各周期的噪声信号，这两步过程如下式(8)、(9)：

{N_stand(t)}＝{S_base(t)}-S_stand(t) (9)

接下来则是对移动站信号的噪声压制，首先将各周期的移动站信号S_move(t)减去各周期的噪声信号N_stand(t)得到去噪后的各周期的移动站信号，再进一步叠加后使用指数加权滑动平均(Exponential Weighted Moving Average，EWMA)压制随机噪声得到最终结果。

指数加权滑动平均(Exponential Weighted Moving Average，EWMA)是统计学中常用的一个方法，它在机器学习和信号分析领域常用来分析数据。EWMA给时间序列每个样本点赋予一个随时间衰减的权重，从而在减少均值计算量的同时剔除一些随机误差。对于一时间序列，其t时刻的加权平均可以表示为：

EWM_t＝(1-β)·(β⁰·e_t+β¹·e_t-1+…+β^t-1·e₁) (10)

式中β为人为设置的权重，它的大小影响加权窗的宽度，设置为0时相当于没有加权取滑动加窗，等于0.98时相当于取了近50个时间点的加权平均，一般情况我们取为0.9，即近10个时间点的加权滑动加窗。相比于传统的滑动平均，EWMA更节省计算内存和时间，同时对时间序列中的随机干扰也能起到一定的压制作用。

于是，本实施例中的移动站信号的最终压制结果S_vmd-_lnc(tj)通过以下两式求得：

{newS_move(t)}＝{S_move(t)}-{N_stand(t)} (11)

最后，根据实际地形情况结合提取的信号特征对参数进行修正，起始频率ω₀会影响初始运算速度，K和α取值不当易导致信号失去准确度，根据实际地形情况下信号的频谱分析来调整校正这几个参数，可以防止过拟合和去噪不完全，有助于提高处理结果的精确度。

VMD方法的参数将决定重构信号的准确度，不同的噪声对应的去噪参数是不同的，不同信号的主频不同对应的参数也是不同的。所以应当结合实际情况测试其参数并校正后再进行噪声提取工作，即在完成测区调研后选择合适地址布设一组仪器完成发射采集，然后根据原信号频谱和地区噪声特征选择最合适的一组参数，具体的，应当依据主频设定去噪起始频率，然后再根据噪声特征和频谱分析确定K和α取值。在完成测试并得到参数后再以这些参数为参考对后续大规模数据进行去噪工作。当原信号的频谱分析结果发现主频较少且噪声也比较单一时，K应取值小一些，此时K值略大于主频数和噪声主频之和即可，1/α值可以设置的稍大一些(如设置为5000～10000)，起始频率ω₀设置成接近主频的大小；当主频较为多且噪声组成相对复杂时，K值需要相应取大一些，但是不应该超过20，1/α值可以设置的小一点，但是不应该小于100，过小会使得信号失真，起始频率ω₀可以在主频附近测试适当调整，以最快计算值为佳。

基站应当距离发射源更近，可以压制干扰，这样就可以先提取出发射信号，然后利用基站信号减去发射信号就得到噪声信号，移动站信号减去噪声信号就得到有用信号。

实验例

以下通过具体的仿真实验说明本实施例的技术效果，实验中初始模型发射极距设置为80m，偏移距为450m，地层模型参考实验地形建立成A型地电模型：第一层为厚度300米的50Ω低阻层，第二层为厚度500m的200Ω地层，最后一层为厚度20000米的5000Ω高阻层。现在分别按图2设置发射及基站和移动采集站，基站离发射源近一些，移动站可以有一个或者多个，根据前文的发射极距和地层模型(发射距80米，A型地层)，分别按偏移距320米和偏移距480米设置基站和移动采集站，对正演得到的电场值添加一定信噪比的噪声，再对合成的模拟信号使用VMD-LNC方法验证该方法的去噪效果。噪声是使用Python3.8的numpy库中random.randn语句编写代码生成的，参考实际实验中单点的采集时长，共合成了100时长为10秒的周期，每个周期都为功率谱恒定且概率密度符合正态分布高斯白噪声，可控源电磁一维正演的电场Ex分布结果，信号相关情况在如图3～图6中示出。由于移动站距离发射源较远，因此其正演得到的电场值也会在数值上偏小一些，这就导致对于基站来说相对幅度较小的白噪声添加到移动站信号中后噪声几乎将信号完全覆盖，已经无法辨识出实际信号。

对于模拟信号正演生成并添加噪声后按周期保存即视为完成了信号采集。实验中，发射一般使用双极性TD50发射波形，这样可以有效压制仪器和电极噪声影响。采集则使用凤凰公司V8按3200点每周期采样，同时需要记录场源坐标，接收电极坐标等信息。随后开始发射并记录数据，发射电流电压信息和接收信号同时记录。

实验中，为了便于评价，采用相似度(NCC)和信噪比(SNR)来衡量去噪效果和噪声强度，NCC和SNR的计算方式如下：

上面二式中f(n)是原信号，g(n)是VMD-LNC处理得到的信号。利用这两个公式计算得基站信号添加噪声后信噪比SNR是14.150，相似度NCC是0.982，与实际添加的情况相符。移动站信号相对较小，故信噪比和相似度都会更差一些，计算得信噪比SNR为1.588，相似度为0.789，属于受到强干扰污染信号。

对于受到强干扰污染的基站信号和移动站信号，首先将基站信号叠加(实测信号一般是挑选周期后再叠加)得到基叠加信号。

下一步进行VMD分解重构得到去噪后的基叠加信号，对于模拟添加的白噪声，由于其不存在固定的频率，因此使用VMD分解时，主要调整的参数是二次惩罚因子α(在程序中为了便于计算使用了其倒数)、起始频率ω₀和分解次数K。实验发现对于只含白噪声的基站信号，在一定范围内增大K值可以使处理得到的信噪比变大，而α值的本质是约束因子，因此α过大会使得分解后的IMF失真、而过小会出现过拟合效应，同样使得信噪比下降。本实验中，选择了α为3000，K值一般情况下应在2-8之间，此处为2，起始频率ω₀固定为200。

图7和图8分别是实验模拟的信号分解后得到的IMF₁和IMF₂的频谱图，由于VMD算法自身具有抗噪性，IMF₁和IMF₂的频谱图中已经几乎看不到高斯噪声的影响，因此不去除任何IMF即可直接将信号重构。在得到去噪的基叠加信号后，将各周期的基站信号减去去噪的基叠加信号得到各周期的噪声信号；最后将各周期的移动站信号减去各周期的噪声信号得到去噪后的各周期的移动站信号，再进一步叠加后使用指数加权滑动平均(ExponentialWeighted Moving Average，EWMA)压制随机噪声得到最终结果。

最终经过VMD-LNC-EWMA处理得到的结果与原信号以及加噪后的信号的对比如图9所示，可以发现信号得到了恢复，由信噪比公式(13)和相似度公式(14)计算得其信噪比为40.117，与初始添加噪声的信噪比1.588相比有较大提升，其相似度高达0.999952，这些指标都充分说明了VMD-LNC-EWMA方法在压制强干扰方面有着优秀的表现。

另外，使用VMD-LNC-EWMA处理模拟数据后得到的结果与其他方法相比信噪比得到大幅提高，而且方法简便易与实现。同时相比于之前的叠加后处理技术，VMD-LNC-WMA技术根据噪声本身特征进行了针对性压制，由于噪声压制机制的存在，信号采集的时候不再需要采集大量周期的信号，缩减了采集时间，节省了采集成本。

图10显示了以上述实验为基础仿真得到的，经过VMD-LNC-EWMA处理后最终得到信号的信噪比与VMD参数之间的关系的示意图。从图中明显的可以看到前述的结论，在一定范围内增大K值可以使处理得到的信噪比变大，而α过大会使得分解后的IMF失真、而过小会出现过拟合效应，同样使得信噪比下降。表1中示例性的列举了本例中的多组参数组合及其带来的信号信噪比值。

表1信噪比SNR与VMD参数之间的关系对照表

实测例

在通过模拟实验验证了VMD-LNC-EWMA方法的随机噪声压制性能之后，本例中使用该方法对某实验场地采集的电性源瞬变电磁数据进行了处理实验，电性源瞬变电磁信号的部分周期时间序列形态如图11所示。观测实测电性源瞬变电磁信号可以发现由于噪声的干扰，采集到的时间序列不能观测到衰减的二次场形态，真实信号被噪声所掩盖，并且噪声表现出规律网状分布，是典型的50Hz工频干扰。

电性源瞬变电磁方法从提出到现在已经获得了较大的发展应用，野外采集的数据主要包括：自然噪声、文化干扰、系统干扰、激发极化(IP)效应等，对于这些干扰传统处理一般分为叠前处理和叠后处理。叠前处理是对信号周期回放后选择性叠加，如果条件允许也会在这一步进行时域递归滤波；叠后处理则主要用来处理工业干扰、DC漂移和IP效应等，完成这些工作后再对数据进行加窗平滑。区别与传统的方法，LNC技术中基站得到的基叠加信号是决定LNC处理结果优劣的关键点，显然传统的选择性叠加和数字滤波器压制在处理随机噪声方面还有所不足。为了解决这一问题，本例采用VMD-LNC-EWMA方法来代替传统数字滤波器进行基站信号的处理。

如图12所示，对原移动站信号进行傅氏变换可以观测到原信号含有的噪声信号大部分成分是50Hz周期噪声干扰，同时从图12中可以发现，还有一些随机噪声干扰。于是，采用VMD算法处理时，50Hz工频干扰对应的IMF需要去除。

图13、图14中示出了对本例中采集到的基站信号的进行分解后得到的IMF₁和IMF₂的频谱图，从图中可以看到，IMF₂的主频为50Hz，因此被视为应当被剔除的IMF。

如图15所示，VMD-LNC-EWMA方法处理后的结果的频谱显示不仅没有50Hz干扰，随机干扰振幅也接近于0。

为提升提取的信号质量，本发明使用变分模态分解(VMD)替代数字滤波器。VMD是由维纳滤波器演变而来，相比于数字滤波器在频率域滤波，VMD直接在时间域对时间序列进行处理，充分考虑了信号的连续性。如图16所示，与传统处理方法相对比，VMD-LNC-WEMA处理方法不存在晚期的震荡效应，中期不平稳部分也压制了随机干扰。

通过实际实验采集数据得到的结果和正演结果对比，可以发现实际数据处理结果和预估实验地形的参数的正演模型构造的地层情况大致是吻合的，从侧面说明了处理结果的准确性。对比原始数据、常规处理(EWMA加窗)和VMD—LNC-EWMA处理结果可以发现，常规处理已经成功压制了50Hz工频干扰，但是常规处理方法对于1ms和10ms之间的干扰处理不完全，同时尾部有震荡，但是VMD—LNC方法则将这些问题有效解决了。综上所述，VMD-LNC方法对噪声的压制是有效的，通过VMD-LNC方法提取的信号也更加接近准确的信号时间序列。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (8)

1.一种基于变分模态分解的时域电磁法近场噪声压制方法，其特征在于，包括，

设置一固定基站和至少一个移动站同步采集信号，利用固定基站和移动采集站之间噪声的相关性来进行噪声压制，具体如下：

①得到基站信号b(t)，选取多个周期并使用VMD算法处理后重构得到去噪的基信号S，其中包括：挑选基站信号的n个周期进行叠加并取平均值，得到基叠加信号：

用于VMD算法处理；

VMD算法中，在获得的VMD分解各阶IMF_i后，使用傅立叶变换计算原信号的主频和IMF_i主频，根据频谱结果去除随机噪声对应的IMF_i，再将特征明显的其他噪声IMF_i去除，最后把剩余的IMF_i重构为去噪后的基叠加信号S_stand(t)如下：

其中m为剔除噪声对应的IMF_i后余下的IMF_i的数量；

②将信号b(t)减去基信号S得到噪声n(t)：

n(t)＝b(t)-S；

③将移动采集站的信号m(t)减去噪声得到LNC压制后的信号c(t)：

c(t))＝m(t)-n(t)。

2.根据权利要求1所述的方法，其特征在于，所述周期的筛选标准如下：

对各个周期与基站所有周期信号的均值做相关性计算，挑选相关系数绝对值大于0.95的周期，若无周期满足，则挑选相关系数最大的那部分周期。

3.根据权利要求1所述的方法，其特征在于，还包括对各周期的LNC压制后的信号c(t)进行指数加权滑动平均的步骤。

4.根据权利要求1所述的方法，其特征在于，步骤①的VMD算法，根据以下约束变分模型和增广拉格朗日函数求解各阶窄带分量IMF：

其中f是输入的原信号，u_k是分解得到的k阶IMF，ω_k是各IMF的中心频率，α是二次惩罚因子，λ(t)是拉格朗日乘子。

5.根据权利要求1所述的方法，其特征在于，步骤②中，在得到去噪的基叠加信号S_stand(t)后，将各周期的基站信号S_base(t)减去去噪的基叠加信号S_stand(t)得到各周期的噪声信号N_stand(t)如下：

{N_stand(t)}＝{S_base(t)}-S_stand(t)。

6.根据权利要求5所述的方法，其特征在于，步骤③中，将各周期的移动站信号S_move(t)减去各周期的噪声信号N_stand(t)得到去噪后的各周期的移动站信号newS_move(t)如下：

{newS_move(t)}＝{S_move(t)}-{N_stand(t)}。

7.根据权利要求6所述的方法，其特征在于，按以下方式对去噪后的各周期的移动站信号newS_move(t)进行指数加权滑动平均得到最后的信号结果

S_vmd-lnc(t_j)：

其中β为权重。

8.根据权利要求7所述的方法，其特征在于，β取0.9。