基于变分模态分解和奇异谱分析的GPR信号去噪

戴前伟，丁 浩，张 华，张 豪

1.中南大学地球科学与信息物理学院，长沙 410083 2.有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083

0 引言

探地雷达(ground-penetrating radar, GPR)作为一种快速、高效的地球物理探测方法，通过接收并分析反射回波中携带的目标体信息来研究其形态大小、空间分布等物理特征[1-2]。因其具有操作简便、采集快速、适应性强和无损探测等优点，已在水文现象监测[3]、城市管线测量[4]、隧道检测[5-6]和地质灾害调查[7-8]等领域广泛应用。

在GPR的采集过程中，混入的噪声信号具有随机、非平稳的特点，对此，国内外学者已经采取了许多去噪技术来针对性地提高GPR信号的分辨率[9]。小波阈值去噪是一种常用的表征非平稳信号频率随时间变化的去噪分析方法。王勇等[10]对GPR信号的小波基进行了寻优，戴前伟等[11]提出一种基于小波分解频带划分与系数相关性相结合的阈值计算方法，邹根等[12]对阈值选取方式和阈值函数进行了优化，这些都取得了比常规小波阈值处理更好的效果。张梦殊等[13]在小波阈值去噪的基础上，针对低频小波系数进行维纳滤波，联合小波阈值和维纳滤波重构得到去噪信号，更接近原始信号。然而，使用小波函数降噪过程中，小波基、分解层数和小波阈值受人为选择影响容易使结果出现偏差。王宪楠等[14]发现在Shearlet域，GPR数据能够得到更加稀疏的表示，再通过阈值去噪的方法将随机噪声去除。Liu等[15]提出了一种主成分分析-奇异值分解(PCA-SVD)混合降噪方法，通过对信号子空间进行奇异值分解和主成分分析，将子空间的杂波噪音用低秩表示并去除，从而达到去噪的效果。刘雷等[16]在时域、频域分别使用主成分分析去噪算法对合成雷达剖面分别进行去噪分析，提高了GPR数据的信噪比。但基于主成分分析的方法在选择前K个特征值的问题上，对去噪的效果影响很大。K越小，目标信号越能够完整保留，抑制效果越差；K越大，抑制效果越好但会损失部分目标信号。K值的选择目前很大程度上仍依赖于经验计算或者根据数值模拟的结果进行人工选择。经验模态分解(empirical mode decomposition, EMD)[17]基于信号自身时间尺度特点，无需设定基函数或特征值参数，具有自适应时频处理信号的优势，适合于处理非线性非平稳信号。杨建军等[18]论证了EMD方法在GPR信号处理上更优于小波阈值去噪方法；冯德山等[19]使用EMD方法实现了在低信噪比GPR信号中突出雷达剖面中异常体特征，提高了信噪比和数据解释精度。然而，对于复杂的信号，采用EMD方法进行信号分解时往往会出现模态混叠及端点飞翼的现象。为了解决这一问题，许军才等[20]应用了集成经验模态分解(ensemble empirical mode decomposition, EEMD)对GPR信号进行去噪处理，得到了优于EMD的去噪效果。EEMD是人为添加白噪声，再进行平均得到各分量，解决了EMD的模态混叠问题，但在重构信号中仍存在残留白噪声且循环过程中增大了算法的计算量[21]。变分模态分解(variational mode decomposition, VMD)[22]是近年提出的一种时频分析方法，在故障诊断[23]、特征提取[24]、微震信号[25]及GPR[26]数据处理中都得到了成功应用。但该方法在模态数的选择和有效信号的提取上仍依赖于去噪效果对比和人工经验分析选择，且去噪后的信号中仍含有部分中低频噪声，使得信号整体存在振荡现象。

因此，本文提出一种自适应的VMD和奇异谱分析(singular spectrum analysis, SSA)[27]相联合的去噪方法，进行高效准确的GPR信号去噪。通过合成Ricker子波实验和合成雷达剖面模拟实验，对比EEMD、VMD和本文方法，检验该方法的正确性和有效性。最后，将该方法应用于实际工程的GPR信号去噪处理，以去除噪声，突出目标信号。

1 方法理论

1.1 变分模态分解

VMD是一种时频分析方法。该算法假定所有分量都是集中在各自中心频率附近的窄带信号，因此，根据分量窄带条件建立约束优化问题，从而估计信号分量的中心频率，能够将多分量信号一次性分解成多个单分量调幅调频信号，避免迭代过程中遇到的端点效应和模态混叠问题[22]。变分问题的构造步骤[28]如下。

1)采用Hilbert变换计算每个模态分量的解析信号，并获取其对应的单边频谱。

2)通过在各模态中加入指数项来调整各自估计的中心频率，将每个模态分量频谱调整至其对应的基带。

3)构建约束变分模型的表达式：

(1)

式中：uk为模态分量；ωk为uk的中心频率；K为模态数；δ(t)为脉冲函数；t为时间域变量；j为虚数单位；s为原信号；*为卷积运算。

4)在变分问题的求解中，引入二次惩罚因子α和Lagrange算子λ(t)，将式(1)转化为无约束方程：

(2)

式中，〈,〉为点积。

5)利用交替方向乘子迭代算法，结合傅里叶等距变换，优化得到相应的模态分量和中心频率，并搜索增广拉格朗日表达式的鞍点。式(2)的解为：

(3)

6)重复步骤5)直到满足迭代条件，并输出K个分量，判别式为

(4)

式中，ε为判断精度，且ε>0。

1.1.1 自适应模态数选择

进行变分模态分解需要预先设置分解的模态数，K值过大会导致过量分解，使得有效信号频率体现在多个分量中；过小会导致分解不完全，有效信号混叠在噪音信号中。所以选择最合适的模态数，是运用VMD方法的关键。

经变分模态分解后各分量是正交关系，因此各分量的能量和应与原信号的能量相等。但在实际操作过程中往往会因无法正交分解而产生能量的损失，损失的能量越少说明分解的正交性越好。因此，可以通过对信号损失能量进行定量分析选取分解的模态数K[29]。

设采集到有限时长L的GPR单道信号为si(i=1,2,3，…，L)，经变分模态分解后其各分量为uk(k=1,2,3,…,K)，则定义该信号总能量为

(5)

分量信号能量之和为

(6)

选取K=2,3,4,…,N，N为K的最大取值，综合统计规律和计算量，本文中N=15。此时分析其能量损失比η=(E-EIMF)/E，当η达到最小时，K为最合适的模态数。

1.1.2 皮尔逊相关系数法

相关系数法是通过计算各分量与原信号之间的皮尔逊相关系数，来判定有效信号与噪声信号。为了避免幅值较小的有效信号被去除，需先将所有的分量与原信号进行归一化处理。其归一化后的相关系数为[30]

(7)

以相关系数的标准差为判断阈值，其计算式为

(8)

若ri>ξ，则认为第i个分量为有效信号；否则为噪声信号。

1.2 奇异谱分析

经变分模态分解得到的有效信号中，往往会残留一部分噪声，导致信号振荡。为了进一步去除这部分噪声信号，改善振荡现象，引入SSA算法。

SSA是一种处理非线性时间序列数据的方法，通过对所要研究的时间序列的轨迹矩阵进行分解，计算时间序列矩阵的特征值，将较小特征值认定为噪声子空间，而保留较大特征值重构信号，从而实现对噪声的抑制[27]。

SSA可以概述为以下两个步骤[31]。

1)分解：对有效信号分量进行SSA去噪。首先需要选择合适的窗口长度m，通常情况下取m

(9)

式中，g为任一有效信号分量。计算其协方差矩阵S：

(10)

式中，0≤τ≤m-1。对S进行奇异值分解，得到m个特征值并排序：λ1≥λ2≥…≥λm≥0。其对应的特征向量为U1,U2,…,Um，且有：

(11)

(12)

式中：Uk,j为第k个特征值对应的特征向量中的第j个元素；0

2)重构：制定合适的阈值。选取前v个特征值，对应的{λv}构成新的子集，则可将m×n的矩阵g重构成长度为L的时间序列。具体重构式为：

(13)

式中，xk,i为第k个特征值对应的重构信号中的第i个元素。

1.3 算法步骤

运用VMD和SSA进行信号去噪的算法步骤如下：

1)预先设置K的取值范围(K=2, 3, …,N，N=15)，计算其对应的η，并选取最小η对应的K值。

2)对含噪信号进行变分模态分解得到k个本征模函数分量。

3)将含噪信号及各分量进行归一化，根据式(7)计算各分量的皮尔逊相关系数，并根据式(8)计算划分有效信号和噪声信号的阈值。

4)对有效信号分量进行奇异谱分析，进一步去除中低频噪声。

5)将经步骤4)处理后的分量进行重构，得到去噪后的信号。

综上所述，本文提出的基于VMD和SSA的GPR去噪方法的算法流程如图1所示。

附伪代码如下。

基于VMD和SSA的GPR去噪方法输入:GPR信号(M道),VMD参数,SSA参数输出:去噪信号for i=1: Mfor K=2: N 变分模态分解; 计算原信号能量及分量总能量;end计算损失能量比最小时的K值;变分模态分解;归一化:原信号及分量;根据式(7)和式(8)计算相关系数及阈值;识别并提取有效信号分量;SSA去噪;重构;end

2 数值实验

2.1 合成Ricker子波实验

为验证本文方法对低频含噪信号和高频含噪信号的去噪性能，分别用EEMD方法、VMD方法和本文方法(VMD-SSA)对合成Ricker子波信号进行去噪对比。合成信号由在20、50、80 ns处的Ricker子波信号(100和900 MHz)和高斯白噪声组成。合成信号表达式s(t)为：

(14)

式中：f0为Ricker子波中心频率；RN(t)为高斯白噪声。

为了更好地对比去噪效果，通过计算去噪前后的均方根误差(root mean squared error,ERMS)和信噪比(signal to noise ratio,RSN)进行定量分析：

(15)

式中，s′为降噪后的信号。RSN越大，且ERMS越小，表明信号中含有效成分越多，去噪效果越好。

图1 基于VMD-SSA的GPR去噪算法流程图

为了探究VMD模态数的自适应性选择，在100和900 MHz的Ricker子波合成信号中添加同一噪声信号，计算其能量损失比随模态数的变化曲线如图2所示。由图2可知，随着K增大，η先减小到最小值，其后总体呈现上升趋势。其中:对于100 MHz的Ricker子波信号，η值在K为5时达到最小；对于900 MHz的Ricker子波信号，η值在K为3时达到最小。表明对于含相同噪声的不同信号，依据能量损失比最小可以得到最适合该信号的K值。

向100 MHz的Ricker子波合成信号中添加不同信噪比的噪声信号，统计其规律如表1所示。随着信号含噪增多，K值呈现上升趋势，且其对应的η值也逐渐增大。

图3和图4分别为100和900 MHz Ricker子波信号去噪效果以及对应的时频谱图。对于加噪的Ricker子波信号，EEMD、VMD和本文方法都能有效压制噪声，恢复信号。从图3a、b、c、d和图4a、b、c、d中可知：针对不同频率的信号，EEMD和VMD两种方法去噪效果不稳定，表现为对于低频信号，前者相比于后者振荡现象更小，更贴合无噪信号，对于高频信号则相反，且此时EEMD重构信号中含噪更多；而本文方法对于不同主频Ricker子波信号都能最大程度压制噪声，使得信号整体更加光滑，更能凸显出有效信号。从对应的时频谱中图3e--h、图4e--h可以看出，经3种方法去噪后，杂波区域都有改善，且信号的中心频率都有恢复。在100 MHz Ricker子波去噪后的时频谱中，EEMD方法相比于VMD方法噪点更少，在900 MHz Ricker子波去噪后的时频谱中，结论恰好与之相反；而本文方法无论对于高频还是低频信号，杂波区域噪点都是最少的，中心频率所对应的能量更加突出，去噪效果最好。这与时域的结论是一致的。结合表2中信噪比分析可知，与EEMD和传统VMD方法相比，经本文方法处理后的信噪比最大提升13.587 8 dB。本文方法在高频和低频情况下都得到了最高的信噪比和最低的均方根误差，表明本文方法去噪效果更好，去噪后的信号更接近原始信号，具备有效性和稳定性。

a. 100 MHz Ricker子波；b. 900 MHz Ricker子波。

表1 100 MHz Ricker子波合成信号不同信噪比下的K值及其对应的η值

2.2 合成雷达剖面模拟实验

建立背景为均匀介质的空间，相对介电常数设置为6，电导率为0.005 S/m，相对磁导率为1。在介质层中设有3个异常体，分别为空气管线、金属管线及含水管线，其半径为0.1 m，管顶距空气层底0.2 m。每个异常体模型参数如表3所示，位置如图5所示。

设定计算区域大小为5 m×1 m，顶部空气层厚度0.1 m，计算网格大小为0.005 m，时间步长为0.01 ns，时间步为2 500步，时窗长度为25 ns。发射天线和接收天线均置于空气层，采取自激自收的收发方式。以布莱克曼-哈里斯脉冲作为激励源，主频为900 MHz，模拟采集200道，道间距为0.025 m，通过时域有限差分(finite-difference time-domain, FDTD)方法正演得到GPR剖面、并添加高斯白噪声后的合成GPR信号如图6a所示，该剖面背景噪声强烈，大部分反射回波都被淹没在噪声中，难以获取有效信息，分辨异常体特征。

图6中b、c、d分别为EEMD、VMD和本文方法去噪结果。经EEMD及VMD方法去噪后，可以从剖面中分辨出空气管线、金属管线、含水管线的上反射界面回波信息及含水管线的多次回波，但背景噪声依然十分强烈；本文方法压制了背景噪声，使得剖面整体更加平滑，突出异常体信息特征，可以清晰地获取到空气管线和金属管线的微弱回波信息，且含水管线的多次回波信息也更加丰富。

EEMD、VMD和本文方法对合成GPR剖面去噪后的信噪比如表4所示，本文方法处理后剖面的信噪比较EEMD和传统VMD方法分别提高了3.765 9 dB和2.655 7 dB，获得了最高信噪比，最大程度压制噪声，突出有效回波信息，这与合成Ricker子波实验结果是一致的。

取本文方法去噪后的GPR剖面第100道展示了更直观的去噪效果，计算η随K的变化曲线如图7所示。当K取3时，η为最小值。

图8中a、b、c、d分别表示第100道合成GPR信号及3种方法去噪结果，且将3～6 ns的波形进行放大显示。由图8b、c可知，EEMD方法和VMD方法在一定程度上都能恢复出该位置的反射回波，但波形振荡依然严重，重构信号中依然含噪较多。由8d可知，经本文方法去噪后的信号在反射回波处与无噪信号波形几乎重合，波形振荡得到了很好的改善，整体波形更加光滑，去噪效果最好。

a. 加噪信号；b. EEMD去噪信号；c. VMD去噪信号；d. 本文方法去噪信号；e. 加噪信号时频谱；f. EEMD去噪信号时频谱；g. VMD去噪信号时频谱；h. 本文方法去噪信号时频谱。

a. 加噪信号；b. EEMD去噪信号；c. VMD去噪信号；d. 本文方法去噪信号；e. 加噪信号时频谱；f. EEMD去噪信号时频谱；g. VMD去噪信号时频谱；h. 本文方法去噪信号时频谱。

表2 三种方法去噪效果对比结果

表3 异常体模型参数

图5 正演模型

3 实测数据处理

采用中心频率为400 MHz的GPR天线对广州市某地地下电力套管的目标区域进行扫描探测，得到的GPR剖面如图9a所示。选用点测模式进行测量，道间距为0.05 m，时窗为100 ns。图9b为使用本文方法去噪后的剖面，图10为单道去噪效果。由图9可知，50道左右，40～100 ns处有反射回波信息，可以判断此处存在套管。自60 ns起，存在较强噪声干扰。经本文方法去噪后，60～80 ns处的噪声基本去除，80～100 ns处噪声也明显受到压制，图像中噪声区域整体变得平滑，更加突出有效反射回波信息。为了更好地显示背景噪声的去除效果，截取第100道和第200道单道数据处理结果，VMDK值均为4。由图10可知：前50 ns无噪区域波形基本保持一致，直达波及部分反射回波的信号得到保留；从约60 ns起，杂乱的噪声波形得到了很好的压制，达到了去噪的效果。

a. 加噪剖面；b. EEMD去噪后剖面；c .VMD去噪后剖面；d. 本文方法去噪后剖面。

表4 剖面去噪后的信噪比结果

图7 第100道信号不同模态数对应的能量损失比

a. 加噪信号；b. EEMD去噪信号；c. VMD去噪信号；d. 本文方法去噪信号。

a. 实测剖面；b. 本文方法去噪剖面。

a. 第100道；b. 第200道。

4 结论

1)在变分模态分解过程中，分析了基于能量损失比的分解方法。对于不同中心频率的Ricker子波信号，随模态数增加，能量损失比呈现先减小后增大的规律；对于添加不同强度噪声的同一信号，模态数及最小能量损失比均随信噪比减小而增大。结果表明：针对不同的含噪信号，可以实现模态数的自适应选择，避免了人为选择的主观性。

2)本文利用SSA方法对变分模态分解后的信号进行二次滤波，结果表明，本文方法优于EEMD和传统VMD方法的去噪效果。其中，合成Ricker子波实验结果中，与EEMD和传统VMD方法相比，经本文方法处理后信噪比最大提高了13.587 8 dB，突出有效信号的同时使得信号整体更加平滑；合成雷达模拟实验结果中，本文方法处理后信噪比较EEMD和传统VMD方法分别提高3.765 9 dB和2.655 7 dB，剖面中能够分辨出更微弱的反射回波信息，且背景噪声得到更好的压制。

3)本文方法在实测资料的应用中，能够有效压制背景噪声和杂波干扰，突出有效信号，进一步证实了该方法用于处理GPR信号的适用性和鲁棒性，为GPR实际资料的数据处理提供了有价值的参考。