基于变分模态分解算法的探地雷达信号去噪研究

冉利民， 李健伟， 杜娟， 程思远， 宋二乔， 刘四新

1.中石化经纬有限公司，郑州 450000；2.吉林大学 地球探测科学与技术学院，长春 130026

0 引言

探地雷达通过接收来自地下目标的反射回波探测地下地层，测量信号是一种非线性、非平稳信号，传统数据处理方法已经满足不了人们对探测精度越来越高的要求[1]。前人[2-9]对EMD方法不断进行改进，出现了总体经验模态分解(EEMD)、完全总体经验模态分解(CEEMD)等方法，但仍然存在不足之处，例如容易受到噪声干扰，存在模态混叠现象，依赖经验性等问题。

1998年，Huang et al.[2]提出了Hilbert-Huang Transform(HHT)，用来处理分析非平稳、非线性信号。经验模态分解(EMD)是HHT的关键步骤之一。本征模态函数(IMF)由EMD方法分解得到，对IMF进行希尔波特变换可以得到具有明确物理意义的瞬时信号[3-4]。2009年，Wu et al.[5]提出了EEMD方法，在一定程度上解决了EMD分解过程中出现的模态混叠问题。王宏等[6]采用了EEMD这种改进方法，并将之应用到穿墙雷达数据去噪中，相比于EMD方法，模态混叠现象有了明显改善。许军才等[7]把EEMD方法应用到探地雷达正演模拟数据去噪中，噪声得到有效压制，但仍存在加噪去除不干净的问题。2011年，Torres et al.[8]通过改进EEMD方法，提出了CEEMD方法，有效解决了EEMD加噪去除不干净的问题。雷文太等[9]为解决CEEMD方法中需要人为筛选有效IMF分量的问题，提出一种基于自动反相校正和峰度值比较的雷达数据去噪方法，该方法可以通过独立分量分析来筛选IMF分量。EMD方法虽然被不断研究改进，但仍存在一些缺点。例如EMD方法很容易受到噪声干扰，致使得到错误的IMF分量；EMD方法会出现模态混叠现象，这种情况下，IMF分量就失去了物理意义；EMD方法分解得到的IMF分量，没有严格的数学定义，依赖经验性。

2014年，Dragomiretskiy et al.[10]提出VMD算法。从原理上看，EMD方法是基于时域递归分解的，而VMD方法是基于频率完全非递归分解的，因此VMD方法能较好地解决EMD方法存在的问题[11-12]。VMD方法使用迭代方法来搜索变分模型的最优解，从而确定IMF分量的带宽和中心频率。该方法能自适应地剖分信号频域，从而使IMF被有效地分离[13]，这样得到的IMF具有很好的鲁棒性[14]。VMD现在已经是一种很成熟的方法，被应用到众多领域中，例如大坝信号监测、零件故障检测等[12-20]，在地球物理领域也有应用[21-27]。在地震数据处理中，Xue et al.[21]使用VMD方法将地震信号分解成若干IMF分量，这些分量又是非负平滑变化瞬时频率AM-FM信号，具有窄带特征。经过实际数据和模拟数据对比发现，相比于EMD方法，VMD方法具有更强的局部分解能力，且对噪声更加敏感。相比于小波变换等变换方法，VMD方法可以得到更高空间和时间分辨率的瞬时频谱，从而更有利于识别煤炭地震数据中的层位信息[21]。张杏莉等[22]把能量熵理论与VMD方法进行了结合，提出了一种自适应的微震噪声去除方法。相比于EMD方法，该方法去噪效果得到了显著的提高。

在探地雷达数据处理中，Zhang et al.[26]使用VMD方法把探地雷达信号分解成IMF分量，通过实验室模拟数据和实测数据验证，表明VMD方法可以有效地应用在雷达数据噪声处理中，从而可以更清楚地识别地下异常体。许军才等[27]在VMD方法的基础上，使用样本熵来选择高阶模量，从而成功地去除了探地雷达中的白噪声，降噪效果得到了明显的改善。Cheng et al.[28]把VMD方法应用到探冰雷达中，有效消除了噪音影响，获得了更加清晰的冰层结构。因此，笔者将VMD方法引入到机载探地雷达数据处理中。

本文首先介绍VMD方法的原理，然后通过对正余弦信号与高斯白噪声合成的含噪声信号分别用 EMD、EEMD、CEEMD和VMD算法进行分解，分析对比证明VMD算法的优越性。通过时域有限差分方法(FDTD)对一个水平层状冰盖模型进行模拟得到合成数据，并使用VMD方法对其进行去噪处理，最后使用VMD方法对一个机载雷达实测数据进行处理，经过处理的探地雷达信号峰值信噪比从15.4 dB增加到20.3 dB。

1 变分模态分解算法(VMD)

1.1 VMD 基本原理

VMD算法将一个复杂的解析信号f分解为预设尺度的k个具有特定稀疏特性的固有模态函数IMF[10-11]。IMF可以表示为：

uk(t)=Ak(t)cos(φk(t))

(1)

假设IMF在中心频率ωk(t)附近是有限带宽的，则通过估计每个IMF的带宽，建立使IMF带宽之和最小的变分约束模型：

(2)

式中：uk为分解得到的k个IMF分量，ωk为uk的中心频率，∂t为对函数求时间的偏导数，δ(t)为单位脉冲函数，j为虚数单位，*表示卷积。

引入二次罚函数项和拉格朗日乘法算子来求解该变分约束模型(式(2))，因而获得了扩展的Lagrange表达式：

L({uk},{ωk},λ)=

(3)

式中：α为带宽参数，通过设定足够大的正数，使信号的重构精度得到保证，λ为Lagrange乘子。式(3)可以通过交替算子法来求解，频率域上每个IMF的表达式为：

(4)

1.2 VMD主要步骤

(5)

(6)

(7)

式中：τ为噪声容限参数。当信号中含有强噪声时，为了获得良好的降噪效果，可以设定τ=0。VMD的详细完整算法可在[10-11]中找到。基于模式更新的维纳滤波使得VMD算法对噪声具有较强的鲁棒性。

1.3 含噪信号分析检验

首先合成一个含噪信号，由2个频率分别为5 Hz、20 Hz的正弦波信号和高斯白噪声合成的信号组成。然后应用 EMD、EEMD、CEEMD和VMD算法对含噪信号进行分解得到IMF分量。最后，为更好地比较信噪分离效果，把分解得到的IMF分量用频谱图(图1)和瞬时能谱图(图2)表示。由图1可以看到，5 Hz的频率成分被分解到IMF1分离中，20 Hz的频率成分被分解到IMF2分量中，从低频到高频信号被依次分解，模态混叠现象在VMD方法中被有效地避免。从图2可以看到，VMD方法比另外3种方法有更高的时间频率精确度和分辨率。

a.EMD；b.EEMD；c.CEEMD；d.VMD。图1 各方法的IMF分量频谱图Fig.1 IMF component spectrogram of each method

a.EMD；b.EEMD；c.CEEMD；d.VMD。图2 各方法的IMF分量对应的瞬时能谱图Fig.2 Instantaneous energy spectra corresponding to IMF component of each method

2 探地雷达合成数据处理

我们设计一个由冰下基岩、冰盖以及空气组成的模型。模型的大小设定为 45 m×40 m，空气和冰面分界线为0 m，-5～0 m为空气层，0～35 m为冰层(图3)。

正演方法选用时域有限差分(FDTD)。天线主频选择300 MHz，采样间隔0.101 ns，时间窗选择400 ns。网格尺寸为0.05 m×0.05 m，道间距为0.375 m，采集道数为121道。天线置于距地面 5 m 的空气中，收发距为0.1 m。

使用模拟数据(图4a)进行VMD分解的处理数据。设定k=5，可以得到5个IMF分量，有效信息在图4c和图4f中，噪声信息主要在图4b、图4d和图4e中。

a.介电常数模型；b.电导率模型。图3 合成数据模型(考虑实际情况增加了随机介质的干扰)Fig.3 Synthetic data model

图4 试验数据图和各IMF分量图Fig.4 Test data diagram and each IMF component diagram

表1为实验数据与各个IMF分量的相关系数，用于定量分析各分量的有效信息。通过在模拟数据中加入随机等效介质来产生噪声。当有效信息占主体，噪声信息比例较小时，定量分析以各分量与实验数据的相关系数为标准。如果IMF分量的有效信息占比较少，那么其相关系数也会越小，从而剔除该IMF分量。反之，如果IMF分析的有效信息占比较多，那么相关系数会比较大，该IMF分量需要保留下来。因此，笔者以相关系数>0.6为标准，对IMF分量进行重构。

表1 实验数据与IMF的相关系数

从表1可以看到，IMF2和IMF5的相关系数满足标准，对这两个分量进行重组，如图5所示。对比图5和图4a可以看到，剖面中的噪声明显减少，整体来看更加清晰。因此通过剔除噪声占比较多的IMF分量，重组有效信息占比较多的IMF分量，可以明显地提高信噪比。经过上面的对比，可以证明VMD方法在探地雷达信号去噪中的有效性。

图5 探地雷达模拟数据IMF2与IMF5分量组合图Fig.5 Combination of IMF2 and IMF5 components of ground penetrating radar simulation data

3 探地雷达实测数据处理

实测数据采用的是第32次中国南极科考期间，HiCARS机载雷达测量的F13R47a测线的一部分数据。该实测数据主频为60 MHz，每道采样点数为2 000个点，总共有1 000道。对其进行一些常规探地雷达数据处理，包括去直流、减平均等，其中带通滤波频率分别为1 MHz、20 MHz、150 MHz和300 MHz，常速度偏移中速度为0.167 m/ns。之后，使用VMD方法分解处理完毕后的数据。

图6 使用常规探地雷达数据处理方法处理后的数据(虚线表示的是测试道位置)Fig.6 Data processed by conventional ground penetrating radar data processing method

从图6可以看到，即使经过上述常规探地雷达数据处理，实测数据的信噪比还是比较低，存在较明显的干扰信号，这对等时层判断造成了严重的干扰。使用VMD方法对第405道实测数据进行处理，设定k=5，可以得到5个IMF分量(图7)。同样，使用VMD方法对整个剖面数据进行处理，k仍为5，得到5个IMF分量图(图8)。从图8可以看到，有效信息在图8c和图8f中，噪声信息主要在图8b、图8d和图8e中。

采取与合成数据中一样的评价方法，使用相关系数对各IMF分量中有效信息进行评价。基于该实测数据为有效数据的前提下，那么有效信号占比要明显大于噪声信号。因而把实测数据与IMF分量的相关系数作为评价标准是可行的。表2为实测数据与各个IMF分量的相关系数，计算获得的相关系数与根据IMF分量得到的结论是一致的。

图7 第405道数据和各IMF分量图Fig.7 The 405th channel data and each IMF component

图8 原始实测数据与各IMF分量Fig.8 Original measured data and each IMF component

表2 原始实测数据与IMF的相关系数

同样以相关系数>0.6作为标准，选择IMF1和IMF4分量进行重组，重组后的效果见图9。经过VMD分解，再挑选合适的IMF分量进行重组的结果剖面，要比原始实测数据的剖面更加清晰，可以较好地识别等时层与基岩面(图9)。我们使用峰值信噪比(PSNR)定量评价该方法的去噪能力。重组后数据的PSNR是20.3 dB，而原始数据PSNR只有15.4 dB，存在明显的差距。因此，可以进一步证实VMD方法在雷达数据去噪方面的有效性。

4 结论

(1)VMD方法可以有效避免模态混叠现象，相比EMD、EEMD和CEEMD等方法有更高的时间频率精确度和分辨率。

(2)VMD方法相比其他3种方法，可以更好地识别探地雷达数据中等时层与基岩面。原始数据只有15.4 dB的峰值信噪比，经过VMD方法处理，峰值信噪比可以提高到20.3 dB。

图9 原始数据剖面(a)与IMF1和IMF2重组后的剖面(b)Fig.9 Original data section (a) and reorganized section of IMF1 and IMF2(b)