1 引言

探地雷达(Ground Penetrating Radar, GPR)，又称透地雷达，通过发射天线向地下发射宽带电磁波，接收天线接收散射回波进行地下未知区域的无损探测，主要应用于工程地质探测和安全防卫等场合^[1,2]。宽带电磁波在地下传播的过程中，受到复杂多变的地下介质结构和外界随机噪声等多种因素的影响，会发生衰减、频散和多次散射干扰，这在很大程度上限制了探地雷达的探测分辨率和最终解释效果^[3]。为了提高GPR回波的数据质量，提高探测效果，需要对GPR回波信号进行去噪处理，抑制GPR回波信号中的噪声和干扰。

针对不同的探测场景和后续数据处理需求，研究人员提出了多种GPR信号去噪算法，包括快速傅里叶变化法、小波变化法、独立分量分析(Independent Component Analysis, ICA)、经验模态分析(Empirical Mode Decomposition, EMD)、完全总体经验模态分解(Complete Ensemble Empirical Mode Decomposition, CEEMD)等^[4–11]，这些算法各有特点。段炜^[8]提出对GPR信号使用小波变换特性对高频系数置零来去除噪声，实验验证了小波变换法比快速傅里叶变化法去噪效果好，但是也造成了高频细节丢失，同时该算法严重依赖于小波的形态。基于ICA理论, 张春成等^[9]提出了一种针对步进频率GPR的ICA与中值滤波处理相结合的去噪算法，但是没有解决ICA的相位不确定性问题；陈玲娜^[10]提出对GPR回波进行分频处理，通过选择固定频段的本征模态函数(Intrinsic Mode Function, IMF)进行重构从而抑制高频杂波和低频隐失波，再从原始信号中去掉高阶主成分分量进行去噪，但该算法需要人工剔除噪声后再重构信号，效率较低；Li等^[11]提出了一种基于CEEMD的GPR信号处理方法，通过对GPR回波信号进行CEEMD分解，人工方式进行噪声辨识和去除，进而实现信号去噪。运用Hilbert-Huang变换，证明了该算法比经验模态分解和总体经验模态分解的谱分辨率更高，但是也需要通过人工判读的方式对各IMF进行判别，再进行信号重构。

针对上述算法中相位不确定的问题和CEEMD分解后需要人工判别信号分量的问题，本文提出一种基于自动反相校正和峰度值比较的GPR回波信号去噪方法。GPR含噪回波信号与随机噪声进行随机拟合，ICA分解，得到两路峰度值不等的信号，根据峰度值的高低分别判定为信号和噪声。信号分量进行自动反相校正，再进行CEEMD分解，得到多个IMF分量。将噪声分量的峰度值作为阈值，对信号分量分解后得到的多个IMF分量进行自动分类。高于该阈值的IMF分量视为信号分量，累加实现信号重构，完成去噪处理。分别用仿真数据和实测数据对所提的去噪算法进行了实验验证，结果显示，该方法在自动相位校正和自动分量判别方面大大提高了运算效率，提高了去噪效果。

2 基本原理 2.1 独立分量分析ICA

ICA是基于信号的高阶统计量，研究信号间的独立关系，将数据变换到相互独立的方向上，使经过ICA算法变换所得到的各个分量之间不仅正交，而且相互独立^[12]。ICA处理的流程图如图1所示。

在图1中， ${{X}} = \left\{ {{x_1},{x_2}, ·\!·\!· ,{x_n}} \right\}$ 是 $n$ 个独立的原始信号 ${{S}} = \left\{ {{s_1},{s_2}, ·\!·\!· ,{s_n}} \right\}$ 通过混合矩阵 ${{A}}$ 的观测信号，则观测信号和原始信号之间的关系如式(1)所示。

其中， ${{A}} = \left\{ {{a_1},{a_2}, ·\!·\!· ,{a_n}} \right\}$ 为混合矩阵。

观测信号 ${{X}}$ 经过解混合矩阵 ${{W}}$ 得到信号为 ${{Y}} = \left\{ {{y_{\text{1}}},\ {y_{\text{2}}},\ ·\!·\!·,\ {y_n}} \right\}$ ，则

其中， ${{W}}$ 为解混合矩阵。

为了使信号 ${{Y}}$ 无限逼近原始信号 ${{S}}$ ，则由式(2)可得

所以ICA就是求解一个最优的矩阵 ${{W}}$ ，使得式(3)的误差最小。

2.2 完全总体经验模态分解(CEEMD)

CEEMD是一种基于EMD的通过噪声辅助的数据分析方法^[13–17]。将固定比例的高斯白噪声添加到原始信号中，形成新的待分解信号，对新的待分解信号进行EMD分解得到IMF分量，用不同的白噪声分别进行 $N$ 次分解，并将 $N$ 个1阶的IMF进行整体均值，即式(4)所示。

其中， ${I_1}$ 为原始信号 $x$ 的第1阶IMF分量， ${\omega _i}$ 是均值为0、方差为1的高斯白噪声， $\varepsilon $ 为加入噪声的比例， $N$ 为加入不同噪声实现EMD分解的次数， ${E_i}$ 表示产生第 $i$ 个IMF分量。

然后计算1阶残差

其中， ${I_1}$ 为原始信号 $x$ 的第1阶IMF分量。

接着继续分解 ${r_1} + \varepsilon {E_1}\left[ {{\omega _i}} \right]$ ，重复上述步骤，从而得到第2个IMF分量，再计算2阶残差，再分解。从而得到

继续进行筛选，直到残差的极值个数不超过两个(残差为常量或单调函数)时停止，得到

其中， $R$ 是最终的残差， $K$ 是原信号分解得到的IMF分量的数目。因此，原信号可以表示为：

所以CEEMD是在EMD的基础上加入高斯白噪声再将分解成从高频到低频的分量，它比EMD分解效果更好，并且解决了模态混叠的问题。

3 基于自动反相校正和峰度值比较的GPR去噪算法

本文提出一种基于自动反相校正和峰度值比较的GPR信号去噪算法，流程图如图2所示。

具体来说，包括以下步骤：

步骤1 含噪GPR信号与随机噪声的拟合　

设GPR获取的单道含噪数据为 $x\left( m \right)$ ，首先产生一道与GPR回波信号 $x\left( m \right)$ 长度相同的随机噪声信号 $n\left( m \right)$ ，设这两个信号的长度都为 $N$ ，将这两个信号组合成 $2 \times M$ 的矩阵，记为 ${{K}}$ ，该矩阵的第1行为 $x\left( m \right)$ ，第2行为 $n\left( m \right)$ 。为了让含噪信号和随机噪声混合成两个相差不大的混合矩阵，所以生成一个2×2的矩阵 $\left[ {\begin{array}{*{20}{l}}1&1\\1&2\end{array}} \right]$ ，记为 ${{L}}$ ，将 ${{L}}$ 左乘 ${{K}}$ ，得到矩阵 ${{U}}$ ，其第1行数据记为 ${U_1}$ ，第2行数据记为 ${U_2}$ 。本文采用的是基于负熵最大化的ICA算法，该算法是通过观察信号进行中心化、白化等预处理过程，然后计算最优的分离矩阵，将 ${U_1}$ 和 ${U_2}$ 运用ICA算法，因为混合信号是由两个不同信号混合而成，根据ICA算法的分离特性也会分离成两个不同的信号，输出峰度值不相等的两道长度为 $M$ 的信号，分别计算这两道信号的峰度值，峰度值计算公式如下：

其中， $x\left( m \right)$ 代表信号序列， $\bar x$ 代表 $x\left( m \right)$ 的均值。峰度值高的记为 $x'\left( m \right)$ ，峰度值低的记为 $n'\left( m \right)$ 。

步骤2 ICA分离和自动校正

判断 $x'\left( m \right)$ 与 $x\left( m \right)$ 是否反相。ICA分离混合信号得到噪声和信号，那么分离得到的信号跟原GPR含噪回波信号应该是同相，只是幅度不同，那么同相信号相减的绝对值会比原信号的绝对值要小，如果比原信号的绝对值要大，说明信号出现了反相，判断公式是

其中， $\left| \cdot \right|$ 表示对1维数组 $x'\left( m \right) - x\left( m \right)$ 中各元素取绝对值， $\max \left( \cdot \right)$ 表示取1维数组 $\left| {x\left( m \right)} \right|, $ $m = 1,2,·\!·\!·,M$ 中各元素的最大值。如果成立，置 $\alpha = - 1$ ；否则，置 $\alpha = 1$ ，输出信号为 $\alpha \cdot x'\left( m \right)$ 。

步骤3 CEEMD分解、自动阈值判别和重构恢复

对 $\alpha \cdot x'\left( m \right)$ 进行CEEMD分解，得到 $P \, $ 个IMF分量。计算每个IMF分量 ${y_p}\left( m \right),p = 1,·\!·\!·,P \, $ 的峰度值，记为 ${k_p},p = 1,·\!·\!·,P\, $ 。计算ICA算法后输出信号 $n'\left( m \right)$ 的峰度值，记为 $k$ 。因为ICA分离出的噪声 $n'\left( m \right)$ 与ICA分离出的信号 $x'\left( m \right)$ 中含有的噪声相似，所以峰度值相同 ，因为噪声是随机信号，所以信号的噪声比噪声的峰度值要大，因而将满足 ${k_p} > k$ 的所有信号分量累加，作为原始回波信号 $x\left( m \right)$ 去噪后的信号，记为

经过上述3个步骤，对原始含噪的GPR信号进行了去噪处理，压制了噪声分量，便于进行后续的成像等处理。

4 仿真与实测数据去噪实验

下面分别运用仿真和实测数据对所提的GPR去噪算法进行验证。仿真数据采用gprMAX2.0^[18]软件生成，实测数据为GPR对公路的探测数据。

4.1 仿真数据的处理

首先设置仿真场景，使用gprMAX2.0软件生成不含噪的GPR信号，人为地加入噪声，将该含噪信号作为原始数据运用本文所提的去噪算法进行去噪处理。去噪后的信号再与gprMAX仿真生成的不含噪信号进行对比分析，验证本文所提算法的有效性。

设定雷达探测场景如图3所示，地下介质有3层，第1、第2层的厚度均为15 cm，第1、第2、第3层的相对介电常数分别为20、10和15，一个半径为4 cm的金属管埋在第2层中，金属管的中心距第2层上表面6 cm，发射和接收天线位于金属管中心的正上方，距地表的高度为5 cm。发射信号为中心频率为900 MHz的Ricker子波，时窗取为40 ns，采样点数为6784，运行gprMAX得到GPR回波如图4所示。

该信号为无噪的回波信号，为模拟生成含噪的GPR信号，人为地加入了随机白噪声，幅度为4.3328，方差为1.0513, SNR=19 dB，加入噪声后的GPR回波信号如图5(a)所示。按照所提的去噪算法，生成与含噪GPR信号长度相等的随机噪声，如图5(b)所示。

将图5(a)含噪GPR信号和图5(b)随机噪声进行随机拟合，即与一个2×2的矩阵 $\left[ {\begin{array}{*{20}{c}} 1&1 \\ 1&2 \end{array}} \right]$ 相乘，得到两道混合信号，随机拟合结果如图6(a)和图6(b)所示。图6(a)和图6(b)输入ICA算法进行处理，得到两道输出信号，该两道信号的峰度值分别为81.6651和2.9763，峰度值高的记为 $x'\left( t \right)$ ，峰度值低的记为 $n'\left( t \right)$ 。 $x'\left( t \right)$ 和 $n'\left( t \right)$ 分别如图7(a)和图7(b)所示。

对图7(a)信号进行判断是否信号反相，该信号未反相，保持不变。再基于CEEMD算法的处理，得到图8所示的15个IMF分量。

对图8中的每一个IMF分量计算峰度值，峰度值分别为[2.1453, 2.1983, 2.7546, 2.8556, 75.0391, 41.8097, 25.3151, 16.4506, 9.3599, 7.3725, 3.2763, 3.7535, 3.9701, 1.6372, 1.7456]。计算ICA算法分离出来的图7(b)噪声信号的峰度值为2.9763。因此，将峰度值低于2.9763的IMF分量剔除，将剩余IMF分量进行累加得到的重构信号作为去噪后的信号，重构信号如图9所示。

为衡量本文所提算法的效果，采用均方误差评价因子，计算方法如下：

式中， $s\left( m \right)$ 代表原始的gprMAX模拟生成的不含噪信号， $z\left( m \right)$ 代表本算法去噪后的信号。

本仿真实验中，含噪的GPR回波的峰值为4.0418。运用本文提出的去噪算法得到去噪后的信号，用式(12)计算的去噪均方误差为0.0000283。

为定量对比分析本文所提的去噪算法与常规CEEMD去噪算法的性能，通过对原始gprMAX仿真得到的无噪回波信号加入不同幅度的噪声信号，生成不同SNR的含噪GPR回波信号，SNR的设置范围为[0, 20]，点数为41点。每一SNR情况下，分别运用本文的去噪算法和常规CEEMD算法对该含噪GPR回波进行去噪处理，得到去噪回波，再与图4所示的无噪GPR回波进行对比分析，计算相对误差值。遍历整个SNR设置区间，得到各SNR情况下本文所提去噪算法和常规CEEMD去噪算法的误差曲线，如图10所示。

从图10中可见，随着SNR的增大，所提的去噪算法和常规CEEMD去噪算法得到的结果的均方误差都随之降低。但在同一SNR情况下，本文去噪结果的相对误差更低，去噪性能高于常规CEEMD的去噪算法。

本文所提的算法，采用了自动相位校正和自动阈值判别，大大提高了计算效率。在CPU为Pentium^® Dual-Core 2.00 GHz，内存为2.00 GB的PC上运行，上述单道GPR回波去噪处理耗时358.35 s。虽然对单道GPR回波去噪处理耗时需358.35 s，但常规的CEEMD算法在耗时266.55 s的基础上还需要进行人工判别IMF分量，相对而言，本算法无需每道数据处理过程中的人工判别，效率更高。

4.2 实测数据的处理

采用中心频率为1 GHz的GPR对湖南省的邵怀高速进行沿线扫描探测，截取了其中一段数据，横向扫描长度为1.5 m，采样点数为187，纵向时窗为40 ns，采样点数为220。原始GPR B-Scan如图11所示。

运用本文所提的算法对该B-Scan剖面的187道数据进行逐道去噪处理，得到未去噪的B-Scan1和去噪后的原始B-Scan2剖面如图12所示。

通过图12对比，可以看出本文算法在去噪效果方面的有效性。为了进一步体现其去噪效果，从去噪后中选取出第50和第100道数据与未进行去噪的原始数据进行对比，如图13所示。

通过图13中第50和第100道数据与未进行去噪的原始数据进行对比，从对比结果来看，本算法的去噪效果不错，再一次验证了本文算法去噪效果方面的有效性。

5 总结

GPR回波信号的去噪处理决定着后续的数据解译质量。针对现有去噪算法的相位不确定和人工判别IMF分量的问题，本文提出基于自动反相校正和峰度值比较的去噪算法，设计了相位判别因子，实现了ICA分解后的信号相位自动判别和校正，设计了基于峰度值比较的IMF分量自动筛选，阈值选择为ICA算法分离出来的噪声信号的峰度值。该算法避免了ICA分解后的相位不定性，且在CEEMD分解后无需传统的人工方式进行各IMF分量的筛选。通过对仿真和实测数据的处理，验证了本文所提算法的有效性。后续研究工作将集中在CEEMD算法运算效率的优化上，进一步提高去噪效率。