CN117932317B - 一种基于vmd-mad的长周期大地电磁信号降噪方法 - Google Patents

Abstract

本发明提供了一种基于VMD‑MAD的长周期大地电磁信号降噪方法，涉及信号处理领域，包括：对获取的原始长周期大地电磁信号使用变分模态分解法进行分解，得到模态分量信号；使用相关系数法自适应筛选分类模态分量，去除其中的噪声分量信号，保留剩下的模态分量信号；预先设置区间长度，利用中位数绝对偏差法对保留的模态分量信号自适应识别离群点，并使用上下限阈值替换该离群点；对已完成离群点替换的模态分量信号进行信号重构，得到重构降噪信号。本发明的有益效果是：使用相关系数法去除大尺度干扰噪声，使用MAD法消除尖峰干扰，有效减小了这两类噪声干扰对长周期大地电磁信号的影响，提高了信号质量和长周期大地电磁法的勘探结果的准确性。

Description

一种基于VMD-MAD的长周期大地电磁信号降噪方法

技术领域

本发明涉及信号处理技术领域，尤其涉及一种基于VMD-MAD的长周期大地电磁信号降噪方法。

背景技术

上世纪50年代，Tikhonov和Cagniard提出了用于地球物理勘探的大地电磁(MT)方法，该方法以自然电磁场作为场源，通过长时间测量地表相互正交的电场分信号和磁场分量信号，根据采集的大地电磁信号可以计算得到地下电性结构。因此，此方法被广泛应用于矿产资源勘探、油气田普查、岩石圈电性结构探测等领域。由于长周期大地电磁法具有更低的频率，相比于其他物探方法具有更大的勘探深度，因此，大力发展长周期大地电磁法等勘探技术研究，有利于促进我国在深部地球勘探中的实际应用，有利于人类对地球深部的了解。

然而天然的大地电磁信号十分微弱，极容易受到各种噪声的干扰，尤其是长时间测量过程中存在的大尺度噪声干扰和尖峰噪声干扰，这些噪声信号在时域上的幅值会比大地电磁信号高出几个数量级，这使得计算的曲视电阻率-相位曲线会出现严重畸变和失真，大大影响测量结果的可靠性。国内外许多学者对如何去除有效大地电磁中噪声信号的研究十分深入，比如有最早期的用于去除同源相关噪声的远程参考大地电磁测深法，但该方法实施中选择一个合适距离的参考点是很难实现的；Vozoff和Hermance将最小二乘估计法用于阻抗估计以去除非相关噪声，该方法在处理过程中会出现额外飞点导致偏离真值；用于去除非高斯噪声的稳健估计法，该方法通过自适应加权的方式能有效减少飞点和跳点的出现，但是当噪声干扰过大时，该方法的处理会给出错误的阻抗结果。因此，结合现代数据处理技术对实测长周期大地电磁信号进行降噪处理已经成为电磁勘探领域的热门研究课题，高质量的长周期大地电磁信号是准确获得真实地下电性结构的必要前提。

发明内容

为了去除长周期大地电磁信号中的大尺度噪声和尖峰噪声，本发明提供了一种基于VMD-MAD的长周期大地电磁信号降噪方法，能够有效去除实际采集的大地电磁信号中存在的大尺度噪声干扰和尖峰噪声，极大地提高了长周期大地电磁野外采集的数据质量，对长周期大地电磁勘探方法的发展有着较大的现实意义。

该方法主要包括：

S1：对获取的原始长周期大地电磁信号使用变分模态分解法进行分解，得到模态分量信号；

S2：基于步骤S1得到的模态分量信号，使用相关系数法自适应筛选分类模态分量，用于去除其中的噪声分量信号，保留剩下的模态分量信号；

S3：预先设置区间长度，利用中位数绝对偏差法对步骤S2中保留的模态分量信号自适应识别离群点，并使用上下限阈值替换该离群点；

S4：对步骤S3中已完成离群点替换的模态分量信号进行信号重构，得到重构降噪信号，即得到降噪后长周期大地电磁信号。

进一步地，步骤S1具体为：

S11：构造变分问题

将原始长周期大地电磁信号分解为有限个的固有模态分量，其本质就是构造变分问题和求解变分问题，其约束条件为每个模态分量中心频率带宽之和最小，且各模态分量重构后就是原始信号，其约束性变分问题可表示为下式：

其中，{u_q}＝{u₁,u₂...,u_K}表示将原始长周期大地电磁信号分解成K个模态分量信号，k＝1,2,...,K，f(t)为原始长周期大地电磁信号，{ω_q}＝{ω₁,ω₂...,ω_K}表示K个模态分量各自的中心频率，δ(t)表示脉冲函数；表示偏导计算，对函数求时间t的导数，*表示卷积，j表示虚数单位；

S12：变成非约束变分问题

引入二次惩罚因子α和拉格朗日乘法算子λ(t)来求解变分问题，使其变成非约束变分问题：

其中，K表示模态分量信号的个数；f(t)为原始长周期大地电磁信号；α为二次惩罚因子；<>表示内积运算；

S13：交替更新，求解变成非约束变分问题

通过交替更新来寻求变分问题的最优解，使得找到非约束变分问题的极小值点；

其中，对更新的表达式为：

对中心频率更新的表达式为：

对λ(w)ⁿ⁺¹更新的表达式为：

其中，τ表示步长更新系数，λ(t)表示Lagrange乘法算子，ε表示给定的判别精度，n表示当前迭代次数，f(ω)、λ(ω)ⁿ⁺¹分别表示第n+1次迭代时f(t)、λⁿ⁺¹(t)的傅里叶变换，ω表示中心频率；

直到满足约束条件时终止迭代，条件式迭代停止条件如下式所定义：

其中，表示第n次迭代时的傅里叶变换。

进一步地，步骤S2中，使用相关系数法筛选分类模态分量的过程为：

基于步骤S1中分解得到的K个模态分量信号，将各模态分量信号和原始长周期大地电磁信号进行相关系数的计算，得到每个模态分量的相关系数，模态分量IMF＝(x₁,x₂,...,x_m)与原始长周期大地电磁信号f(t)＝(y₁,y₂,...,y_m)的相关系数计算方法如下所示：

其中，ρ_xy表示相关系数，m表示原始长周期大地电磁信号的长度，μ_x表示一个模态分量IMF＝(x₁,x₂,...,x_m)的平均值，μ_y表示原始长周期大地电磁信号f(t)＝(y₁,y₂,...,y_m)的平均值，ρ_xy的绝对值越大，表示模态分量信号中包含的信息与原始长周期大地电磁信号的大尺度漂移噪声干扰越相关，通过比较各个模态分量计算出的相关系数，选择其中与大尺度漂移噪声干扰最大的相关系数所对应的模态分量信号，将其视作大尺度噪声干扰分量直接去除，保留剩下的模态分量信号。

进一步地，步骤S3中，实现离群点识别的过程为：

将离群点定义为在一定长度的时间序列中与中位数相差超过3倍MAD₀的元素，MAD₀的计算方法由下式所定义：

MAD₀＝c*median(|A_i-median(A)|)

其中，MAD₀表示绝对中位数，c为一个常数，i为需要预先设置中位数计算的区间长度，A表示当区间长度为i时，在IMF＝(x₁,x₂,...,x_m)中所选取的时间序列，A_i表示选取时间序列中的第i点，median()表示对括号内的时间序列求解中位数，median(A)表示对时间序列A求取其中位数；

通过上述计算，有效识别出长度为i的时间序列中存在的离群点，通过依次平移该时间长度，就可以完整地遍历模态分量IMF＝(x₁,x₂,...,x_m)中的时间序列，识别出所有离群点。

进一步地，步骤S3中，实现离群点替换的过程为：

基于遍历IMF＝(x₁,x₂,...,x_m)所识别出来的离群点，使用长度为i的时间序列中上限阈值或下限阈值来替换所识别出的离群点；

其中，上限阈值为：

Upper_Thr＝median(A)+3*MAD

其中，下限阈值为：

Lower_Thr＝median(A)-3*MAD

当长度为i的时间序列中某个离群点的数值大于该时间序列的中位数时，该离群点的数值大于上限阈值，直接使用上限阈值替换掉该离群点；当某个离群点的数值小于该时间序列的中位数时，该离群点的数值小于下限限阈值，直接使用下限阈值替换掉这类离群点。

进一步地，步骤S4中，重构的降噪信号为：

其中，IMFs_k为完成离群点替换的剩余模态分量，Out(t)为重构的降噪信号。

一种存储设备，所述存储设备存储指令及数据用于实现一种基于VMD-MAD的长周期大地电磁信号降噪方法。

一种基于VMD-MAD的长周期大地电磁信号降噪设备，包括：处理器及所述存储设备；所述处理器加载并执行所述存储设备中的指令及数据用于实现一种基于VMD-MAD的长周期大地电磁信号降噪方法。

本发明提供的技术方案带来的有益效果是：本发明充分考虑了长周期大地电磁信号中存在的两种噪声干扰的特征，针对两种噪声的不同特征，采用分步处理的方法，使用相关系数法去除大尺度干扰噪声，使用MAD法消除尖峰干扰，有效减小了这两类噪声干扰对长周期大地电磁中有效信号的影响，减少了原始信号中的噪声，提高了采集信号的数据质量，也提高了长周期大地电磁法这一物探方法勘探结果的准确性。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种基于VMD-MAD的长周期大地电磁信号降噪方法的流程图。

图2是本发明实施例中VMD分解方法的示意图。

图3是本发明实施例中原始仿真信号和含噪信号时域与频域的示意图。

图4是本发明实施例中VMD分解模态分量时域的示意图。

图5是本发明实施例中剩余模态分量MAD去除离群点时域的示意图。

图6是本发明实施例中重构信号时域与频域对比示意图。

图7是本发明实施例中各降噪方法处理结果时域波形对比示意图。

图8是本发明实施例中硬件设备工作的示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

请参考图1，图1是本发明实施例中一种基于VMD-MAD的长周期大地电磁信号降噪方法的流程图，首先对输入的大地电磁信号使用变分模态分解法(VMD)进行分解，接着使用相关系数法自适应筛选分类模态分量，去除其中的噪声分量信号，保留剩下的模态分量，预先设置区间长度，利用中位数绝对偏差法(MAD)对保留的模态分量信号自适应识别离群点，并使用上下限阈值替换该离群点，对完成离群点替换的各模态分量信号进行信号重构，得到重构降噪信号，并输出降噪后长周期大地电磁信号，最终完成对长周期大地电磁信号的降噪。该长周期大地电磁信号降噪方法的具体实现步骤如下：

S1：对获取的原始长周期大地电磁信号使用变分模态分解法(VMD)进行分解，得到模态分量信号；具体为：

S11：构造变分问题

将输入信号分解为有限个的固有模态分量(IMF)，其本质就是构造变分问题和求解变分问题，其约束条件为每个模态分量中心频率带宽之和最小，且各模态分量重构后就是原始信号，其约束性变分问题可表示为下式：

其中，{u_q}＝{u₁,u₂…,u_K}表示将原始长周期大地电磁信号分解成K个模态分量，k＝1,2,...,K，f(t)为输入信号，即原始长周期大地电磁信号，{ω_q}＝{ω₁,ω₂…,ω_K}表示K个模态分量各自的中心频率，δ(t)表示脉冲函数；表示偏导计算，对函数求时间t的导数，*表示卷积，j表示虚数单位。

S12：变成非约束变分问题

引入二次惩罚因子α和拉格朗日乘法算子λ(t)来求解变分问题，使其变成非约束变分问题，因此，上述的描述公式可以重新描述为下式：

其中：K表示模态分量个数；f(t)为输入信号；α为二次惩罚因子；<>表示内积运算。

S13：交替更新，求解变成非约束变分问题

通过交替更新λ(ω)ⁿ⁺¹来寻求变分问题的最优解，使得找到非约束变分问题的极小值点；

其中，对更新的表达式为：

对中心频率更新的表达式为：

对λ(w)ⁿ⁺¹更新的表达式为：

其中，τ表示步长更新系数，λ(t)表示Lagrange乘法算子，ε表示给定的判别精度，n表示当前迭代次数，f(ω)、λ(ω)ⁿ⁺¹分别表示第n+1次迭代时f(t)、λⁿ⁺¹(t)的傅里叶变换，f(t)为原始信号，ω表示中心频率。

直到满足约束条件时终止迭代，条件式迭代停止条件如下式所定义：

其中，表示第n次迭代时的傅里叶变换。

在实施例中，上述变分模态分解算法的具体实现流程如图2所示。首先，需要预设模态数K，并根据经验值设置二次惩罚因子α；然后初始化更新计算u_k、ω_k，对变量n和k进行自增，判断k的值是否到达预设值，然后得到λ后，判断是否满足迭代停止条件，当满足条件后最终得到K个模态分量。

具体地，在实施例中，将含有多个不同频率不同振幅的正弦波信号叠加形成合成信号，用该信号模拟长周期大地电磁信号，向合成信号中随机加入大幅度方波干扰和尖峰干扰来模拟采集信号中存在的大尺度干扰和尖峰噪声，将信号转换到频域中，能更准确看出信号的被干扰程度，具体参见图3中原始仿真信号和含噪信号时域与频域示意图。

进一步的，经过步骤S1中对上述的模拟信号使用变分模态分解法(VMD)进行分解，在本次分解中，将模态数K和惩罚因子α分别预设为5和1000，经过图2中所描述的分解算法流程将图3中含噪信号分得到五个模态分量信号，在时域上的波形如图4所示。

S2：基于步骤S1得到的模态分量信号，使用相关系数法自适应筛选分类模态分量，用于去除其中的噪声分量信号，保留剩下的模态分量信号；具体如下：

基于步骤S1中分解的K个模态分量信号，将各模态分量信号和原始输入信号进行相关系数的计算，得到每个模态分量的相关系数，模态分量IMF＝(x₁,x₂,...,x_m)与原始信号f(t)＝(y₁,y₂,...,y_m)的相关系数计算方法如下式所示：

其中，ρ_xy代表计算出的相关系数，m表示输入信号的长度，μ_x表示一个模态分量IMF＝(x₁,x₂,...,x_m)的平均值，μ_y表示输入信号f(t)＝(y₁,y₂,...,y_m)的平均值，ρ_xy的绝对值越大，表示模态分量信号中包含的信息与输入信号的大尺度漂移噪声干扰越相关，通过比较各个模态分量计算出的相关系数，选择其中与大尺度漂移噪声干扰最大的相关系数所对应的分量信号，将其视作大尺度噪声干扰分量直接去除，保留剩下的模态分量信号。

本实施例中，按照相关系数定义式的计算方法，分别计算图4中5个模态分量与图3中含噪信号的相关系数，其相关系数计算值如表1所示。

表1各模态分量相关系数计算值

从表中可以看出IMF5模态分量与原始信号的相关系数是最大的，这在图4分解出的模态分量中也可以看出，IMF5分量的时域波形与加入的大尺度噪声干扰非常相似，几乎不包含有效信号。依据步骤(S2)，本实施例中，直接将IMF5模态分量看成噪声模态分量，去除IMF5模态分量，保留剩余的四个模态分量作为有效分量信号，其中分解出的其他四个模态分量中出现了很多尖峰噪声的干扰，这些噪声的处理会在步骤(S3)进行。

S3：基于步骤S2中相关系数法筛选保留的模态分量信号，预先设置区间长度，利用中位数绝对偏差法(MAD)对保留的模态分量信号自适应识别离群点，并使用上下限阈值替换该离群点，具体如下：

S31：离群点识别

经过步骤(S2)的自适应筛选，得到保留后的模态分量信号。在实际采集的长周期大地电磁信号中存在很多的大尺度噪声和尖峰噪声，这些其中大尺度噪声经过步骤(S1)的VMD分解后，体现在模态分量中形式就是大尺度的平缓模态分量和大量分布在剩余模态的尖峰信号。通过步骤(S2)中的相关系数法可以筛选识别出前者的大尺度噪声，针对后者，可以使用离群点识别方法来去除。

其中，本方法中将离群点定义为在一定长度的时间序列中与中位数相差超过3倍MAD₀的元素。进一步的，MAD₀的计算方法由下式所定义：

MAD₀＝c*median(|A_i-median(A)|)

其中，MAD₀表示绝对中位数，c为一个常数，本实施例中其值为其中erfcinv()表示逆互补误差函数，i为需要预先设置中位数计算的区间长度，A表示当区间长度i为时，在IMF＝(x₁,x₂,...,x_n)中所选取的时间序列，A_i表示选取时间序列中的第i点，式中median()表示对括号内的时间序列求解中位数，median(A)表示对时间序列A求取其中位数。

通过上述的计算公式，能有效识别出长度为i的时间序列中存在的离群点，通过依次平移该时间长度，就可以完整地遍历模态分量IMF＝(x₁,x₂,...,x_m)中的时间序列，识别出所有存在的离群点。

(S32)离群点替换

基于步骤(S31)遍历IMF＝(x₁,x₂,...,x_m)所识别出来的离群点，为避免离群点对整体信号的影响，需要将这些离群点剔除，并填充上一个合适的数值以保证时间序列的长度不变，本方法中使用长度为i的时间序列中上限阈值或下限阈值来替换所识别出的离群点。

其中，上限阈值由下式所定义：

Upper_Thr＝median(A)+3*MAD

其中，下限阈值由下式所定义：

Lower_Thr＝median(A)-3*MAD

当长度为i的时间序列中某个离群点的数值大于该时间序列的中位数时，根据步骤S31中对离群点的定义，该离群点的数值必然也大于上限阈值，因此，直接使用上限阈值替换掉该离群点。同理，当某个离群点的数值小于该时间序列的中位数时，根据步骤S31中对离群点的定义，该离群点的数值必然小于下限限阈值，直接使用下限阈值替换掉这类离群点。

基于上述对两种离群点情况分类讨论的处理，就可以完成对全部离群点识别与替换。本实施例中，需要对步骤S2中筛选出的4个模态分量分别进行离群点识别与替换，即需要对图4中的IMF1～IMF4进行处理，处理结果如图5所示，从图5中的处理结果可以看出，图4中四个模态分量的尖峰噪声被有效消除。

S4：基于步骤S3中已完成离群点替换的各模态分量信号，对完成离群点替换的模态分量信号进行信号重构，得到重构降噪信号，并输出降噪后长周期大地电磁信号，完成对长周期大地电磁信号的降噪，具体如下：

基于步骤S11中描述的约束性变分问题的定义式，所有分解出的模态分量之和构成了原始的输入信号，可表示成下式：

其中，{u_q}＝{u₁,u₂...,u_K}表示将原信号分解成K个模态分量，f(t)为原始的输入信号。进一步地，在步骤S2中通过相关系数法识别出大尺度干扰噪声的模态分量，并去除了该分量信号，在步骤S31和步骤S32中利用MAD完成了对剩余模态分量的离群点识别与替换，步骤消除了大尺度噪声分解得到尖峰噪声和原本存在的尖峰噪声，这些完成离群点替换的剩余模态分量信号之和，即下式所定义：

其中，IMFs_k为完成离群点替换的剩余模态分量，由于步骤S2中取出来一个噪声模态分量，剩余的模态分量个数变成K-1，Out(t)即为完成降噪过程的输出信号，即重构的降噪信号。

在本实施例中，对图5中完成离群点去除的四个模态分量信号进行信号重构，得到图6中重构降噪信号，从图6中可以看出，重构降噪信号几乎与原始仿真信号相重合，在频域上也能明显地看出与原始仿真信号相同的频谱，有效地去除了加入的模拟干扰噪声。经过计算，重构降噪信号的信噪比从-19.92dB提升到了2.14dB，极大程度上提升了模拟长周期大地电磁信号的数据质量。

表2不同降噪方法对模拟长周期大地电磁信号的性能对比

表2为不同降噪方法在处理模拟长周期大地电磁信号所得结果的信噪比对比，各方法处理结果在时域上的波形对比如图7所示。其中，使用信噪比SNR和均方根误差RMSE作为算法性能的评判指标，EMD为经验模态分解法，EEMD总体经验模态分解法，然后将这两种方法与步骤S3中的MAD相结合对模拟长周期大地电磁信号进行处理。

请参见图8，图8是本发明实施例的硬件设备工作示意图，所述硬件设备具体包括：一种基于VMD-MAD的长周期大地电磁信号降噪设备401、处理器402及存储设备403。

一种基于VMD-MAD的长周期大地电磁信号降噪设备401：所述一种基于VMD-MAD的长周期大地电磁信号降噪设备401实现所述一种基于VMD-MAD的长周期大地电磁信号降噪方法。

处理器402：所述处理器402加载并执行所述存储设备403中的指令及数据用于实现所述一种基于VMD-MAD的长周期大地电磁信号降噪方法。

存储设备403：所述存储设备403存储指令及数据；所述存储设备403用于实现所述一种基于VMD-MAD的长周期大地电磁信号降噪方法。

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

Claims (5)

1.一种基于VMD-MAD的长周期大地电磁信号降噪方法，其特征在于：包括：

S1：对获取的原始长周期大地电磁信号使用变分模态分解法进行分解，得到模态分量信号；

S11：构造变分问题

将原始长周期大地电磁信号分解为有限个的固有模态分量，其本质就是构造变分问题和求解变分问题，其约束条件为每个模态分量中心频率带宽之和最小，且各模态分量重构后就是原始信号，其约束性变分问题可表示为下式：

其中，{u_q}＝{u₁,u₂…,u_K}表示将原始长周期大地电磁信号分解成K个模态分量信号，k＝1,2,...,K，f(t)为原始长周期大地电磁信号，{ω_q}＝{ω₁,ω₂…,ω_K}表示K个模态分量各自的中心频率，δ(t)表示脉冲函数；表示偏导计算，对函数求时间t的导数，*表示卷积，j表示虚数单位；

S12：变成非约束变分问题

引入二次惩罚因子α和拉格朗日乘法算子λ(t)来求解变分问题，使其变成非约束变分问题：

其中，K表示模态分量信号的个数；f(t)为原始长周期大地电磁信号；α为二次惩罚因子；<>表示内积运算；

S13：交替更新，求解变成非约束变分问题

通过交替更新λ(ω)ⁿ⁺¹来寻求变分问题的最优解，使得找到非约束变分问题的极小值点；

其中，对更新的表达式为：

对中心频率更新的表达式为：

对λ(w)ⁿ⁺¹更新的表达式为：

其中，τ表示步长更新系数，λ(t)表示Lagrange乘法算子，ε表示给定的判别精度，n表示当前迭代次数，f(ω)、λ(ω)ⁿ⁺¹分别表示第n+1次迭代时f(t)、λⁿ⁺¹(t)的傅里叶变换，ω表示中心频率；

直到满足约束条件时终止迭代，条件式迭代停止条件如下式所定义：

其中，表示第n次迭代时的傅里叶变换；

S2：基于步骤S1得到的模态分量信号，使用相关系数法自适应筛选分类模态分量，用于去除其中的噪声分量信号，保留剩下的模态分量信号；

使用相关系数法筛选分类模态分量的过程为：

基于步骤S1中分解得到的K个模态分量信号，将各模态分量信号和原始长周期大地电磁信号进行相关系数的计算，得到每个模态分量的相关系数，模态分量IMF＝(x₁,x₂,...,x_m)与原始长周期大地电磁信号f(t)＝(y₁,y₂,...,y_m)的相关系数计算方法如下所示：

其中，ρ_xy表示相关系数，m表示原始长周期大地电磁信号的长度，μ_x表示一个模态分量IMF＝(x₁,x₂,...,x_m)的平均值，μ_y表示原始长周期大地电磁信号f(t)＝(y₁,y₂,...,y_m)的平均值，ρ_xy的绝对值越大，表示模态分量信号中包含的信息与原始长周期大地电磁信号的大尺度漂移噪声干扰越相关，通过比较各个模态分量计算出的相关系数，选择其中与大尺度漂移噪声干扰最大的相关系数所对应的模态分量信号，将其视作大尺度噪声干扰分量直接去除，保留剩下的模态分量信号；

S3：预先设置区间长度，利用中位数绝对偏差法对步骤S2中保留的模态分量信号自适应识别离群点，并使用上下限阈值替换该离群点；

实现离群点识别的过程为：

将离群点定义为在一定长度的时间序列中与中位数相差超过3倍MAD₀的元素，MAD₀的计算方法由下式所定义：

MAD₀＝c*median(|A_i-median(A)|)

其中，MAD₀表示绝对中位数，c为一个常数，i为需要预先设置中位数计算的区间长度，A表示当区间长度为i时，在IMF＝(x₁,x₂,...,x_m)中所选取的时间序列，A_i表示选取时间序列中的第i点，median()表示对括号内的时间序列求解中位数，median(A)表示对时间序列A求取其中位数；

通过上述计算，有效识别出长度为i的时间序列中存在的离群点，通过依次平移该时间长度，就可以完整地遍历模态分量IMF＝(x₁,x₂,...,x_m)中的时间序列，识别出所有离群点；

S4：对步骤S3中已完成离群点替换的模态分量信号进行信号重构，得到重构降噪信号，即得到降噪后长周期大地电磁信号。

2.如权利要求1所述的一种基于VMD-MAD的长周期大地电磁信号降噪方法，其特征在于：步骤S3中，实现离群点替换的过程为：

基于遍历IMF＝(x₁,x₂,...,x_m)所识别出来的离群点，使用长度为i的时间序列中上限阈值或下限阈值来替换所识别出的离群点；

其中，上限阈值为：

Upper_Thr＝median(A)+3*MAD

其中，下限阈值为：

Lower_Thr＝median(A)-3*MAD

当长度为i的时间序列中某个离群点的数值大于该时间序列的中位数时，该离群点的数值大于上限阈值，直接使用上限阈值替换掉该离群点；当某个离群点的数值小于该时间序列的中位数时，该离群点的数值小于下限限阈值，直接使用下限阈值替换掉这类离群点。

3.如权利要求1所述的一种基于VMD-MAD的长周期大地电磁信号降噪方法，其特征在于：步骤S4中，重构的降噪信号为：

其中，IMFs_k为完成离群点替换的剩余模态分量，Out(t)为重构的降噪信号。

4.一种存储设备，其特征在于：所述存储设备存储指令及数据用于实现权利要求1～3任一项所述的基于VMD-MAD的长周期大地电磁信号降噪方法。

5.一种基于VMD-MAD的长周期大地电磁信号降噪设备，其特征在于：包括：处理器及存储设备；所述处理器加载并执行所述存储设备中的指令及数据用于实现权利要求1～3任一项所述的基于VMD-MAD的长周期大地电磁信号降噪方法。