CN109765624A - 一种基于变分模态分解的频率域航空电磁数据去噪方法 - Google Patents

Abstract

本发明公开了一种基于变分模态分解的频率域航空电磁数据去噪方法，主要解决航空电磁数据噪声源多且特性复杂并与有效信号混叠的问题。该方法通过检测频率域航空电磁仪器系统采集的原始数据，得到测线电磁剖面数据；采用变分模态分解对航空电磁剖面数据滤波处理，抑制信号中的人文噪声和天电噪声。通过上述方法步骤，本发明可以同时在频率域和空间域上对天电噪声和人文噪声进行有效抑制。此外，在保证提高数据信噪比的情况下有较高的运算效率。因此，具有很高的使用价值和推广价值。

Description

一种基于变分模态分解的频率域航空电磁数据去噪方法

技术领域

本发明涉及航空电磁技术领域，具体地说，是涉及一种基于变分模态分解的频率域航空电磁数据去噪方法。

背景技术

航空电磁法通称航空电法。是用来快速普查良导电金属矿的航空物探方法。它是通过研究由人工或天然形成的电磁场对地质体感应激发产生的异常场特征和规律(即应用交变电磁场的感应原理)来寻找矿体或解决某些地质问题的。主要是用来快速普查良导性金属矿体(富铜、富铁，对大面积地质填图、圈定近地表的基岩起伏，研究地下水和冰冻层等方面也有一定效果。

航空电磁探测具有勘查速度快、探测范围广等优势，在我国具有广泛的应用前景。但是其机载的飞行探测方式，能够引起发射线圈、接收线圈晃动，同时飞行速度、飞行姿态等变化使得接收信号的信噪比较低，特别是雷电和人文设施所产生的干扰往往幅值较大且和有效信号的频带重叠，严重影响数据质量及成像精度，制约航空电磁探测系统对地下深部异常体及小异常体的反演解释，减小勘探深度。由于航空电磁数据的噪声源多且特性复杂，所以现有的许多滤波方法并不能很好的抑制在时间域及空间域上和瞬变电磁有效信号相重叠的天电噪声和人文噪声。

近年来，随着航空电磁法的推广，国内外学者围绕航空电磁数据开展了大量的处理和解释工作。黄鄂(1998)年提出了经验模态分解(Empirical Mode Decomposition，即EMD)方法，并引入了基于Hilbert变换的Hilbert谱的概念和Hilbert谱分析的方法，即HHT方法，该方法作为一种新的具有自适应性的信号处理方法，已经被用于非平稳、非线性电磁信号分析，包括探地雷达数据、大地电磁数据和航空电磁数据等。但EMD存在一个主要缺陷及模态混叠，当模态混叠出现时，不同频率的分量混叠出现在一个或几个IMF分量之中，这严重影响了基于EMD方法得到的时频分布。为了解决这一问题，针对EMD的改进一直在进行，期间也相继提出了一些其他的方法，如通过添加白噪声，使得经验模态分解比较稳定，以达到在一定程度上降低模态混叠效应、改进信号重构精度的目的，但是经验特性仍然是EMD及其改进版本的本质。基于此Dragomiretskiy于2014年提出了一种完全非递归的信号分解算法即变分模态分解(Variational Mode Decomposition，即VMD)，它可看作是经典Wiener滤波的推广，可自适应产生多尺度频带。该方法已被刘伟等学者成功运用于地震信号分析和探地雷达信号分析。变分模态分解方法具备完善的数学理论基础，就是可以有效缓解EMD信号分解后的模态混叠效应，且运算效率较高，但如何将该方法用于解决航空电磁数据噪声源多且特性复杂并与有效信号混叠的课题任未得到解决。

发明内容

本发明的目的在于提供一种基于变分模态分解的频率域航空电磁数据去噪方法，主要解决航空电磁数据噪声源多且特性复杂并与有效信号混叠的问题。

为实现上述目的，本发明采用的技术方案如下：

一种基于变分模态分解的频率域航空电磁数据去噪方法，包括如下步骤：

(S10)将航空电磁数据输入基于频率域的航空电磁系统；

(S20)航空电磁系统对输入数据进行读取与预处理，获得未去噪的电磁剖面数据，电磁剖面数据即航空电磁法对于地质勘测时在水平方向和竖直方向的磁场响应分量；

(S30)对所获得的电磁剖面数据采用变分模态分解算法进行处理，分离出预设个数的固有模态函数分量；

(S40)对分离出的各个模态分量波形进行Hilbert变换，并将获得的Hilbert变换结果进行有效模态重构，获取所需的去噪后的电磁剖面数据。

进一步地，所述步骤(S20)具体为读取基于频率域的航空电磁系统获得的电磁数据，对数据按若干测点形成的测线分布进行排列整理，并对其进行预处理，其中测点即为按比例尺选定的航空电磁系统观测点位，测线为沿一条直线布置的观测点组成的观测线。通过数据长度和归一化计算其采样频率和采样间隔，获得对应的时间序列和频率序列。

进一步地，所述步骤(S30)具体如下：

(S31)固有模态函数定义为：

u_k(t)＝A_k(t)cos[φ_k(t)] (1)，

式中A_k(t)为u_k(t)的瞬时幅值函数，t为采样时间间隔，φ_k(t)为瞬时相位函数：ω_k(t)＝φ'_k(t)，ω_k为u_k(t)的瞬时频率；

(S32)确定由航空电磁数据分解出的固有模态函数分量个数K，则由变分模态分解建立的约束变分数学模型表达式为：

式中，σ(t)为脉冲函数，j为虚数单位，s.t.为约束目标函数，f代表进行分解的航空电磁数据信号，引入扩展拉格朗日乘子函数，得到下式：

式中α为二次惩罚因子，λ为Lagrange算子，运算符号||||₂表示求该函数的欧几里德范数；

(S33)通过不断更新模态分量u_k、中心频率ω_k和Lagrange算子λ求得式(3)的最优解；更新公式如下：

上式中，ω为采样信号频率，带角标的分别为其对应的傅里叶变换；若ε＞0或n>N，则停止更新；其中，ε为判别约束条件，n为实际迭代次数，N为最大迭代次数。

进一步地，所述步骤(S40)中模态分量的Hilbert变换步骤具体如下：

(S41)根据步骤(S20)中得到的时间序列和频率序列，从而对由步骤(S30)得到的K个固有模态函数模态分析进行Hilbert变化，将每个固有模态函数设为x(t)，Hilbert变换按下式定义：

其中，y(t)即为x(t)的Hilbert变换，x(t)和y(t)构成Hilbert变换对；

(S42)根据式(6)、式(7)得到相应的解析信号z(t)和其幅值谱α(t)为：

相应的瞬时频率如下式所示：

进一步地，所述步骤(S40)中对获得的Hilbert变换结果进行有效模态重构，获取所需的去噪后的电磁剖面数据步骤具体如下：

(S43)根据得到的每一阶固有模态函数做Hilbert变换后的瞬时频率和幅值谱，得到原始信号的Hilbert时频谱：

(S44)对Hilbert时频谱的时间积分，得到Hilbert边际谱h(ω)如下：

式中，T为信号的总长度；边际谱表达了每个频率在全局上的幅度或能量，它代表了在统计意义上的全部累加幅度；

(S45)根据所获得的Hilbert边际谱，得到到原始信号的有效信号和噪音信号所对应频率，结合所获得的各模态分量的时域波形图，确认有效信号对应的模态，将有效模态进行叠加，得到去噪后的航空电磁数据。

与现有技术相比，本发明具有以下有益效果：

(1)本发明采用变分模态分解去除噪音，该方法具备完善的数学理论基础，通过变分模态分解将信号按不同主频进行分解，可以有效减缓信号分解后的模态混叠效应，且运算效率较高，对于航空电磁噪音数据的处理、分析和解释都有较大的优势，能够同时在频率域和空间域上对天电噪声和人文噪声进行有效抑制。

(2)本发明利用Hilbert变换得到信号的边际谱，表示每个频率在全局上的能量贡献，它代表了在统计意义上的整组数据在每个频率点的累积能量分布。该频率成分不必在信号的整个时间段都存在，其发生的时刻，则在Hilbert谱中给出了精确的定位。其不受Fourier变换所需的常幅值和常频率的限制，特别适用于航空电磁数据类非线性、非平稳的电磁信号。

附图说明

图1为本发明中基于VMD分解频率域航空电磁数据去噪方法流程图。

图2为本发明仿真信号及变VMD分解二阶IMF时域波形图。

图3为本发明根据VMD所获得的模态得到的Hilbert信号谱。

图4为本发明实测信号及变VMD分解六阶IMF时域波形图。

图5为本发明实测信号及变EMD分解六阶IMF时域波形图。

图6为本发明根据VMD分解后模态得到的Hilbert信号谱。

图7为本发明根据EMD分解后模态得到的Hilbert信号谱。

图8为本发明所示方法实测数据单道测线进行变分模态分解去噪前后对比图。

图9为本发明实测频率域航空电磁数据去噪前剖面图。

图10为本发明实测频率域航空电磁数据去噪后剖面图。

具体实施方式

下面结合附图说明和实施例对本发明作进一步说明，本发明的方式包括但不仅限于以下实施例。

实施例

如图1所示，本发明公开的一种基于变分模态分解的频率域航空电磁数据去噪方法，包括如下步骤：

(S10)将航空电磁数据输入基于频率域的航空电磁系统；

(S20)航空电磁系统对输入数据进行读取与预处理，获得未去噪的电磁剖面数据，电磁剖面数据即航空电磁法对于地质勘测时在水平方向和竖直方向的磁场响应分量；

(S30)对所获得的电磁剖面数据采用变分模态分解算法进行处理，分离出预设个数的固有模态函数分量；

(S40)对分离出的各个模态分量波形进行Hilbert变换，并将获得的Hilbert变换结果进行有效模态重构，获取所需的去噪后的电磁剖面数据。

所述步骤(S20)具体为读取基于频率域的航空电磁系统获得的电磁数据，对数据按若干测点形成的测线分布进行排列整理，并对其进行预处理，其中测点即为按比例尺选定的航空电磁系统观测点位，测线为沿一条直线布置的观测点组成的观测线，通过数据长度和归一化计算其采样频率和采样间隔，获得对应的时间序列和频率序列。

所述步骤(S30)具体如下：

(S31)固有模态函数定义为：

u_k(t)＝A_k(t)cos[φ_k(t)] (1)，

式中A_k(t)为u_k(t)的瞬时幅值函数，t为时间，φ_k(t)为瞬时相位函数：ω_k(t)＝φ'_k(t)，ω_k为u_k(t)的瞬时频率；

(S32)确定由航空电磁数据分解出的固有模态函数分量个数K，则由变分模态分解建立的约束变分数学模型表达式为：

式中，σ(t)为脉冲函数，j为虚数单位，s.t.为约束目标函数，f代表进行分解的航空电磁数据信号，引入扩展拉格朗日乘子函数，得到下式：

式中α为二次惩罚因子，λ为Lagrange算子，运算符号||||₂表示求该函数的欧几里德范数；

(S33)通过不断更新模态分量u_k、中心频率ω_k和Lagrange算子λ求得式(3)的最优解；更新公式如下：

上式中，ω为频率，带角标的分别为其对应的傅里叶变换；若ε＞0或n>N，则停止更新；其中，ε为判别约束条件，n为实际迭代次数，N为最大迭代次数。

所述步骤(S40)中模态分量的Hilbert变换步骤具体如下：

(S41)根据步骤(S20)中得到的时间序列和频率序列，从而对由步骤(S30)得到的K个固有模态函数模态分析进行Hilbert变化，将每个固有模态函数设为x(t)，Hilbert变换按下式定义：

其中，y(t)即为x(t)的Hilbert变换，x(t)和y(t)构成Hilbert变换对；

(S42)根据式(6)、式(7)得到相应的解析信号z(t)和其幅值谱α(t)为：

相应的瞬时频率如下式所示：

所述步骤(S40)中对获得的Hilbert变换结果进行有效模态重构，获取所需的去噪后的电磁剖面数据步骤具体如下：

(S43)根据得到的每一阶固有模态函数做Hilbert变换后的瞬时频率和幅值谱，得到原始信号的Hilbert时频谱：

(S44)对Hilbert时频谱的时间积分，得到Hilbert边际谱h(ω)如下：

式中，T为信号的总长度；边际谱表达了每个频率在全局上的幅度或能量，它代表了在统计意义上的全部累加幅度；

(S45)根据所获得的Hilbert边际谱，得到到原始信号的有效信号和噪音信号所对应频率，结合所获得的各模态分量的时域波形图，确认有效信号对应的模态，将有效模态进行叠加，得到去噪后的航空电磁数据。

实施实例一

为测试变分模态分解的实际应用效果，采用含有两种成分的测试信号进行分解，该信号的表达式为：

y＝[(1+0.3)sin(2π15t)]cos[(2π50t+0.5sin(2π15t))]+sin(2π150t)

信号由150Hz的正弦波和调制频率为15Hz、基频为50Hz的调频调幅成分叠加而成。选取IMF个数K为2个，据保真度均衡参数α为20000，收敛判断参数Tol为1e-7，DC分量数为0，初始中心频率Init为1。图2即为已知信号和该信号经变分模态分解出来的2阶IMF的时域波形图，共包括3条曲线，图3为根据IMF分量所得到的Hilbert边际谱。从图中可以看出，“c0”曲线代表已知信号，“c1”曲线代表一阶序列，对应Hilbert边际谱中的150Hz的正弦波的高频信号，“c2”曲线为二阶序列，代表Hilbert边际谱中基频为50Hz的调频调幅信号，这说明了VMD的正确性。

实施实例二

选取由频率870Hz的同轴装置实测得到的频率域航空电磁数据进行处理，该数据由九条测线组成，每条测线含4536个测点。选取IMF个数K为6个，据保真度均衡参数α为20000，收敛判断参数Tol为1e-7，DC分量数为0，初始中心频率Init为1。图4即为本发明中根据步骤(S30)，实测数据中单道测线在进行变分模态分解后所得到的时域波形图，图6则为本发明中根据步骤(S40)所得到的单测线的Hilbert信号谱，同时为对比VMD相对于EMD的优越性，图5为实测数据中单道测线在进行EMD分解后所得到的时域波形图，图7则为EMD分解后所得IMF分量得到的Hilbert信号谱。可以看出相较于EMD分解，变分模态分解可以有效地缓解模态混叠现象。在VMD分解所得IMF分量对应的Hilbert信号谱中可以清楚地看出该测线中存在5条能量很强的谱线，其中组成高频谱线的IMF即为要去除的包含噪音成分的模态分量，而EDM分解所得IMF分量对应的Hilbert信号谱中则无法分辨出数据的主要频率成分。图8即为本发明所示方法实测数据单道测线进行变分模态分解去噪前后对比，从图中可以明显看出去噪后的航空电磁数据没有了其他噪声信号的干扰，达到了一个较好的去噪效果。图9、图10为本发明所示方法实测频率域航空电磁数据去噪前后剖面图，其中图9为实测频率域航空电磁数据去噪前剖面图，图10为实测频率域航空电磁数据去噪后剖面图，从图中可以看出，频率域航空电磁数据在水平方向和竖直方向的电磁分量去噪后，噪声源与有效信号的混叠减少，表明了本方法能够有效去除噪声源。

通过上述方法，本发明能够在保留信号幅值、有效信息的同时，很好的抑制航空电磁数据的天电噪声和人文噪声。因此，本发明具有突出的实质性的特点和进步。

上述实施例仅为本发明的优选实施方式之一，不应当用于限制本发明的保护范围，但凡在本发明的主体设计思想和精神上作出的毫无实质意义的改动或润色，其所解决的技术问题仍然与本发明一致的，均应当包含在本发明的保护范围之内。

Claims (5)

1.一种基于变分模态分解的频率域航空电磁数据去噪方法，其特征在于，包括如下步骤：

(S10)将航空电磁数据输入基于频率域的航空电磁系统；

(S20)航空电磁系统对输入数据进行读取与预处理，获得未去噪的电磁剖面数据；

(S30)对所获得的电磁剖面数据采用变分模态分解算法进行处理，分离出预设个数的固有模态函数分量；

(S40)对分离出的各个模态分量波形进行Hilbert变换，并将获得的Hilbert变换结果进行有效模态重构，获取所需的去噪后的电磁剖面数据。

2.根据权利要求1所述的一种基于变分模态分解的频率域航空电磁数据去噪方法，其特征在于，所述步骤(S20)具体为读取基于频率域的航空电磁系统获得的电磁数据，对数据按若干测点形成的测线分布进行排列整理，并对其进行预处理；其中测点即为按比例尺选定的航空电磁系统观测点位，测线为沿一条直线布置的观测点组成的观测线；通过数据长度和归一化计算其采样频率和采样间隔，获得对应的时间序列和频率序列。

3.根据权利要求2所述的一种基于变分模态分解的频率域航空电磁数据去噪方法，其特征在于，所述步骤(S30)具体如下：

(S31)固有模态函数定义为：

u_k(t)＝A_k(t)cos[φ_k(t)] (1)，

式中A_k(t)为uk(t)的瞬时幅值函数，t为采样间隔时间，φ_k(t)为瞬时相位函数：ω_k(t)＝φ'_k(t)，ωk为u_k(t)的瞬时频率；

(S32)确定由航空电磁数据分解出的固有模态函数分量个数K，则由变分模态分解建立的约束变分数学模型表达式为：

式中，σ(t)为脉冲函数，j为虚数单位，s.t.为约束目标函数，f代表需进行分解的航空电磁数据信号，引入扩展拉格朗日乘子函数，得到下式：

式中α为二次惩罚因子，λ为Lagrange算子，运算符号||||₂表示求该函数的欧几里德范数；

(S33)通过不断更新模态分量u_k、中心频率ω_k和Lagrange算子λ求得式(3)的最优解，更新公式如下：

上式中，ω为采样信号频率，带角标的分别为其对应的傅里叶变换；若或n>N，则停止更新；其中，ε为判别约束条件，n为实际迭代次数，N为最大迭代次数。

4.根据权利要求3所述的一种基于变分模态分解的频率域航空电磁数据去噪方法，其特征在于，所述步骤(S40)中模态分量的Hilbert变换步骤具体如下：

(S41)根据步骤(S20)中得到的时间序列和频率序列，从而对由步骤(S30)得到的K个固有模态函数模态分析进行Hilbert变化，将每个固有模态函数设为x(t)，Hilbert变换按下式定义：

其中，y(t)即为x(t)的Hilbert变换，x(t)和y(t)构成Hilbert变换对，p为常数；

(S42)根据式(6)、式(7)得到相应的解析信号z(t)和其幅值谱α(t)为：

相应的瞬时频率如下式所示：

5.根据权利要求4所述的一种基于变分模态分解的频率域航空电磁数据去噪方法，其特征在于，所述步骤(S40)中对获得的Hilbert变换结果进行有效模态重构，获取所需的去噪后的电磁剖面数据步骤具体如下：

(S43)根据得到的每一阶固有模态函数做Hilbert变换后的瞬时频率和幅值谱，得到原始信号的Hilbert时频谱：

(S44)对Hilbert时频谱的时间积分，得到Hilbert边际谱h(ω)如下：

式中，T为信号的总长度；边际谱表达了每个频率在全局上的幅度或能量，它代表了在统计意义上的全部累加幅度；

(S45)根据所获得的Hilbert边际谱，得到到原始信号的有效信号和噪音信号所对应频率，结合所获得的各模态分量的时域波形图，确认有效信号对应的模态，将有效模态进行叠加，得到去噪后的航空电磁数据。