吴学礼 闫枫 甄然 武晓晶 孟凡华

摘 要：针对探地雷达原始图像中存在着大量以地表直达波为主的杂波噪声干扰问题，为了有效提取目标信号，提出了一种基于小波变换和K-Means奇异值分解的自适应双边滤波方法。将原始雷达数据进行小波分解，并应用K-SVD算法变换稀疏编码和更新原子，用更新后的稀疏系数和字典重构小块，将小块进行小波逆变换重构图像，然后对重构图像进行自适应双边滤波得到处理后的图像，以实现小波变换、K-SVD算法和自适应双边滤波方法的优势互补。与其他方法进行对比的结果表明，所提出的方法具有良好的杂波抑制效果，在峰值信噪比和目标图像熵上性能更好。研究结果在实际雷达勘测任务中具有一定的应用价值。

关键词：雷达工程;探地雷达;杂波抑制;小波变换;奇异值分解;自适应双边滤波

中图分类号：TP751;P631.8 文献标识码：A

doi：10.7535/hbkd.2021yx02003

Research on adaptive clutter suppression for ground penetrating radar based on wavelet transform and K-SVD

WU Xueli1，2，YAN Feng1，ZHEN Ran1，WU Xiaojing1，MENG Fanhua1

（1.School of Electrical Engineering，Hebei University of Science and Technology，Shijiazhuang，Hebei 050018，China; 2.Hebei Provincial Research Center for Technologies in Process Engineering Automation，Shijiazhuang，Hebei 050018，China）

Abstract：Aiming at the large amount of clutter and noise interference（mainly direct waves on the surface） in the original image of ground penetrating radar，an adaptive bilateral filtering method based on wavelet transform and K-Means singular value decomposition was proposed in order to effectively extract the target signal.The original radar data was decomposed by wavelet，and then the decomposed data was sparsely coded by K-SVD transform and updated atoms.The updated sparse coefficients and dictionary were used to reconstruct small blocks，and the small blocks were subjected to inverse wavelet transform to reconstruct the image.The reconstructed image was processed by adaptive bilateral filtering to obtain the processed image，so as to realize the complementary advantages of wavelet transform，K-SVD algorithm and adaptive bilateral filtering method.Compared with other methods，the experimental results show that the proposed method has a good clutter suppression effect and better performance in peak signal-to-noise ratio and target image entropy.The research results have certain application value in the actual radar survey tasks.

Keywords：

radar engineering; ground penetrating radar（GPR）; clutter suppression; wavelet transform; singular value decomposition （SVD）; adaptive bilateral filtering

探地雷達（ground penetrating radar，GPR）是近些年发展起来的一种探测地下目标的有效工具，通过从天线发射无线电波并获取从地下异常或掩埋物体反射的回波来探测地下[1]。其能够利用被测区域地下物质的介电性、导磁性不同，根据接收电磁波的动力学和运动学等特征，推断地下介质和目标的结构以及分布特性，从而实现地下检测目标的成像及定位。与其他常规地下探测工具相比，GPR以快速性、高效率、无损直观以及高分辨率等优点在探测工程隐患方面得到了广泛应用及认可[2]。例如在公路工程方面，利用探地雷达可以检测出潜藏于道路下方的空洞危险，以防止发生严重的道路坍塌事故[3]。

探地雷达系统一个主要的问题是杂波的存在，这会严重影响目标的成像和检测能力。造成杂波的原因有多种，例如地面的反射信号、发射天线和接收天线之间的耦合、地表直达波以及地下不同媒介物质引起的背景噪声等。其中直达波对回波数据影响最大，能量最强，是对目标信号影响最大的干扰信号[4]。由于目标信号会受到杂波噪声的干扰，而探地雷达的有效信号又往往会被杂波覆盖，因而有效提取探地雷达中的目标信号对实际勘测以及风险评估具有重要意义，抑制杂波对改善探地雷达回波信号质量和提高图像准确率非常重要[5]。

当前对于如何抑制探地雷达信号中的杂波问题，国内外学者提出了各种降噪方法从雷达回波信号中准确提取有效目标信号。WANG等[6]提出了一种通过信号分解的PCA-SVD（主成分分析-奇异值分解）混合降噪方法，获得了较好的降噪效果。基于子空间的杂波去除方法通过奇异值分解（SVD）、主成分分析（PCA）或独立成分分析（ICA）将原始数据矩阵投影到杂波和目标子空间上[7]。最主要的成分对应于杂波，其余成分用于表示无杂波的目标成分，基于子空间的杂波抑制方法是将目标信号进行数据矩阵的低秩表示，陈瑞鼎等[8]使用卡尔曼滤波器抑制雷达回波信号的杂波。原始信号被分解为目标信号、背景信号和噪声信号。然后，在空间方向上对目标信号进行快速傅里叶变换，获得了较好的目标信号特性。王超等[9]提出了一种基于希尔伯特-黄（HHT）变换的回声信号噪声抑制方法。该方法避免了传统傅里叶变换的约束，但是忽略了HHT变换中2个相邻雷达回波的数据影响以及周围噪声的连续性。2017年，SANTOS等[10]分析了时间反转（TR）技术与各种探地雷达预处理方法结合的性能，改善了地下目标的探测。杨建功等[11]指出，在对探地雷达信号进行去噪处理时，小波分析的去噪效果优于传统的去噪方法。

本文提出了一种基于小波变换和K-Means奇异值分解（K-SVD）的自适应双边滤波算法（Wavelet-K_SVD-BF），通过对探地雷达图像进行仿真实验，证明其有效性，并且通过和其他相关方法进行对比仿真，证明了其不仅可以有效抑制杂波干扰，还有效改善了探地雷达回波信号质量，提高了图像的准确率，为实际勘测工作提供有效帮助。

1 基于小波变换和K-SVD的自适应双边滤波的杂波抑制

1.1 小波变换的基本原理

小波变换是从时域到频域的信号转换，是一种信号的多尺度分析方法，也叫作多分辨率分析[12]。小波变换在信号的高频部分可以取得较好的时间分辨率，在信号的低频部分可以取得较好的频率分辨率，从而能有效从雷达初始信号中提取出有用信息。

连续小波变换的表达式：

WTf（a，b）=〈f（t），φa，b（t）〉=1a∫Rf（t）φ（t-ba）dt ，a≠0，（1）

式中：φa，b（t）是小波基函数，通过母函数φ（t）平移、伸缩得到;a，b分别为尺度因子和平移因子，可以用来调节时窗和频窗的大小，进而实现对时间和频率分辨率的控制，小波变换的时频窗结构可以基本满足对探地雷达信号去除杂波的要求。

连续小波变换逆变换表达式为

f（t）=1CφR1a2WTf（a，b）φa，b（t）dadb，（2）

式中：Cφ=∫Rψ（ω）2ωdω

小波变换处理数据时基函数的选择十分关键，小波波形越接近待处理的瞬态信号波形，处理效果越理想。综合时频域的分辨率来看，DB族小波是比较适合分析探地雷达信号的一种小波基函数[13]。

在大多数情况下，信号主要集中在低频部分，噪声部分基本分布在高频部分。所以，经过小波分解后，信号对应的小波变换系数要大于噪声的小波变换系数[14]。

1.2 K-SVD算法去噪原理

K-SVD算法用于解决图像去噪问题时，含噪雷达图像Y表示为

Y=X+B ，（3）

式中：纯净图像X被随机噪声B污染;X，B和Y代表的都是大小为N×N的方形图像，经过向量化后都是长度为N的列向量。K-SVD去噪算法是假设X是可稀疏的，而随机噪声B是不可稀疏的，利用学习到的字典对信号稀疏表示进行去噪。

图像越大对应的字典就越大。为了降低计算的复杂程度，字典学习通常在大小为N×N的数据块上进行。在去噪过程中首先将图像按照一定步长重叠分为若干小块。无噪图像的第i行第j个块可以表示为xij=RijX，xij∈Rn，Rij代表分块算子。在训练好的字典D∈Rn×p中，可以得到长度为p的稀疏系数。将图像块排列为列向量，每个干净的图像块可以由稀疏系数aij和字典D表示为

xij=Daij ，（4）

则去噪问题转化为求解字典D和稀疏系数的问题。K-SVD算法中字典D和稀疏系数通过求解如下最优化问题得到：

D^，α^ij，X^=argminD^，αij，Xλ‖X-Y‖22+∑ijμij‖αij‖0+∑ij‖Dαij-RijX‖22，（5）

式中：參数λ为正则化参数;实值μij表示第i行第j块对应的参数;0范数表明稀疏系数向量中非零值的个数，0范数值越小，稀疏性越好;α^ij代表无噪图像中第i行的第j个小图像块在字典D下的稀疏系数;式（5）第1项为保真项‖X-Y‖22，保证去噪后的图像要和含噪图像相似，不能丢失过多细节;第2项∑ij‖αij‖0要保证所得稀疏系数尽可能稀疏;最后一项∑ij‖Dαij-RijX‖22是训练好的字典和稀疏系数组合，和对应的图像块之间差越小越好，因为差越小去噪后残留噪声越少。

K-SVD去噪就是要从式（5）中求解字典、稀疏系数和纯净图像。对于一个式子中同时含有3个变量的最优化问题，采用交替方向乘子法（ADMM）。通过变量的交替更新和迭代，将式（5）求解分为3个子问题。

第1个问题是稀疏编码。给定字典和X，求解稀疏系数相当于求解如下代价函数的最小值：

α^ij=arg minα μij‖α‖0+‖Dα-xij‖22 。（6）

对每个小块都进行这步操作，相当于一个滑动窗走遍了含噪图像的所有像素，这个稀疏编码问题可以由正交匹配追踪算法求解。

第2个问题是更新字典。固定稀疏系数和X更新字典：

minD∑ij‖Dαij-RijX‖22 。（7）

第3个问题是计算无噪图像。交替进行以上2步计算直到分块算子遍布所有小块，并且在每个小块的计算过程中字典中原子全部参与了更新之后，便可以求出纯净图像X：

X^=argminx λ‖X-Y‖22+∑ij‖Dα^ij-RijX‖22 。（8）

K-SVD具有良好的图像去噪性能。该降噪算法可以得到使含噪图像尽量稀疏的字典，使得在该字典下图像是可稀疏的，随机噪声是不可稀疏的。

1.3 双边滤波

双边滤波是一种非线性滤波器，可以达到保持边缘、降噪平滑的效果[15]。结合空间信息和亮度相似性对雷达图像作进一步滤波处理，在平滑滤波的同时，能大量保留图像的边缘和细节特征，见式（9）。

f^（x）=∑y∈Ωf（y）ω（x，y）∑y∈Ωω（x，y），（9）

式中：f^为输出图像;f为输入图像;Ω是以像素点x为中心的领域窗口;ω为滤波核。双边滤波的滤波核由2部分乘积组成：空间核φ与值域核ψ。每个像素的权值用其与中心像素点的欧拉距离来计算[16]。2个滤波核通常采用高斯函数形式，如式（10）—式（12）所示。

ω（x，y）=φ（x，y）ψ（x，y） ，（10）

φ（x，y）=exp（-‖y-x‖22σd2） ，（11）

ψ（x，y）=exp（-f（x）-f（y）22σr2） ， （12）

式中：σd为空域高斯函数的标准差;σr为值域高斯函数的标准差。空域滤波系数由像素间的空间距离决定，距离越小，系数越大。值域滤波系数由像素间的相似度决定，像素值越接近系数越大。

在灰度变化平缓区域，值域滤波系数接近1，此时空域滤波起到主要作用，双边滤波器退化为传统的高斯低通滤波器，对图像进行平滑操作。在图像变化剧烈的部分（图像边缘），像素间差距较大，值域滤波起主要作用，因而能有效保持边缘信息。

1.4 基于小波变换和K-SVD方法的双边滤波

传统的SVD算法对探地雷达的A-scan数据时间窗口具有整体性，奇异值的大小反映了A-scan数据道之间在整个时间轴上的相关性，当探测环境较为复杂时，会导致剖面部分区域探测结果之间强相关性。采用传统的SVD算法处理探地雷达数据时，有用的回波信号可能会被滤除。

由于探地雷达回波信号是非平穩信号，因此采用小波变换来处理探地雷达回波数据。小波变换的方法主要是在时频域内滤除随机噪声，小波通过尺度伸缩可以改变时间和频率的分辨率，实现在不同尺度下的时频域处理[17]。二维小波变换将雷达图像分解为近似值、水平细节、垂直细节与对角线细节部分，对近似值进行多次分解，对图像的剖分就越细化。K-SVD算法是在奇异值分解算法上的改进，在字典中的K个代表中，求解时用到SVD分解。对二维小波变换得到各个部分进行K-SVD算法去噪去杂波，各部分被分解成一系列不同的特征成分，其中前1%的奇异值占了数据特征的90%左右，而前10%的奇异值占据了所有奇异值之和的99%以上，文中方法采用前10%的奇异值重构主成分，在迭代次数内稀疏编码和字典更新，最终选取稀疏的部分进行重构，对重构各个部分信号进行二维小波反变换，可以得到方法处理后的探地雷达目标的有效信号，实现二维小波域的K-SVD算法去噪。

基于小波变换的K-SVD算法对雷达原始图像信号去噪除杂波流程图，如图1所示。

基于改进的小波变换K-SVD算法的基本描述如下：

1）对原始探地雷达的回波信号进行小波变换，得到不同分解次数下的近似值、水平细节、垂直细节与对角线信号。根据雷达数据的实际情况，本文采用3次小波分解。

2）对小波分解后的各部分信号进行K-SVD算法处理，在迭代次数内稀疏编码和更新字典元素。

3）用最后一次迭代得到的稀疏系数和字典进行重构。

4）将重构处理后的不同信号进行小波逆变换，并进行进一步自适应双边滤波，得到去噪去杂波的图像。

2 仿真实验与分析

为了验证基于小波变换和K-SVD算法的自适应滤波算法的性能，进行了探地雷达模拟仿真实验。2005年，由爱丁堡大学Dr.Antonis Giannopoulos博士设计的一款探地雷达正演模拟软件GprMax2.0被国内外学者所认可[18]。实验中的雷达模型是用GPRMAX软件进行建模的，GPRMAX是一款基于FDTD算法和PML边界吸收条件的探地雷达正演模拟软件。探地雷达利用主频为数十兆赫兹到上千兆赫兹的高频电磁波，由发射天线Tx向地下发射高频带段脉冲，当遇到电磁性质不连续的地方时发生反射，反射信号由接收线圈Rx接收从而产生雷达图像。地下及地表介质的电磁性质主要参数有电磁率、介电常数、磁导率等，不同介质的波阻抗和对雷达波能量吸收的都有较大差异，即使是同种介质，在不同频率的电磁场作用下的特性也是不同的。

表1是常见介质的电磁常数。

本实验模型采用GPRMAX软件进行建模，选取一个半径为25 mm的完美导弹钢筋，位于厚度为600 mm、深度为75 mm的混凝土板中。设置混凝土的相对介电常数Er=6.0以及磁导率σ=0.01 S/m，用一个ricker源来模拟GPR天线，设置中心频率f=900 MHz，脉冲在比中心频率更高的频率处包含大量能量，所以最高频率应该是脉冲中心频率的3倍～4倍。脉冲频谱中的最高频率fm=3f=2 700 MHz。模型大小选择600 mm×300 mm，发射器和接收器间隔为50 mm，其中发射源为（0.075，0.2525），接收器为（0.125，0.2525）。发射源和接收器高度设置为距接口2.5 mm。具体仿真模型如图2所示。

对雷达图像的二维正演模拟图像进行灰度化处理，对原始雷达图像进行算法降噪去杂波处理。主要目的是除掉原始雷达图像中的背景噪声和地表直达波对于目标信号的干扰，有效提取出图像中目标信号的有用信息。实验分别采用均值滤波法、中值滤波法、小波阈值降噪法、SVD算法以及本文提出的基于小波变换的K-SVD算法进行降噪处理。不同方法进行的去杂波降噪的仿真对比见图3。

通过对比图3不同方法的杂波抑制效果可以发现：1）均值滤波法、中值滤波法、小波阈值变换法、SVD算法以及本文的Wavelet-K-SVD-BF算法都能夠有效去除杂波噪声干扰，其中均值滤波法和中值滤波法去噪去杂波效果相近，但是还存在较多的背景噪声和地表直达波;2）小波阈值变换法和SVD算法的杂波抑制效果相对较好，有效去除了大量背景噪声，但是还存在部分地表直达波;3）从图3 f）可以看出，本文提出的Wavelet-K-SVD-BF算法基本去除了原始雷达图像的背景噪声，地表直达波也基本得到了很好抑制，相比于其他方法，Wavelet-K-SVD-BF算法的杂波抑制效果更好。

本文通过峰值信噪比（PSNR）与目标图像的熵值比较不同算法间的杂波抑制效果。峰值信噪比是最常用的客观评价图像处理效果的量化指标[19]，计算式见式（13）。

PSNR（R0，R）=10×log25521M×N∑Mi=1∑Nj=1（R（i，j）-R0（i，j））2 ，（13）

式中：R0代表大小为M×N像素的原始雷达图像;R代表用算法处理后的雷达图像。在比较图像处理效果时，PSNR的值越大，表示目标图像中的信息损失越少，图像中的目标信号越明显，处理效果越好。

图像熵值Q表示整个目标图像的信息量，一般用图像熵值表示目标图像的清晰程度，图像熵的值越大，图像越混乱，杂波抑制效果越差[20]。熵值越小，探地雷达的图像效果越好。

图像熵值的计算如式（14）所示：

Q=∑Ii=1∑Jj=1K2（i，j）〗2∑Ii=1∑Jj=1K4（i，j） 。（14）

本文算法运行时，计算机处理器是主频为2.40 GHz的Inter Core i5-9300H CPU，运行内存为16 GB。表2从峰值信噪比、目标图像熵以及运行时间3个方面对这几种方法进行了对比。

表2中的数据表明，Wavelet-K-SVD-BF方法有效提高了图片的峰值信噪比，降低了目标图像的熵值，目标信号明显，自适应双边滤波在平滑滤波的同时，保持了图像的边缘和细节特征，图像更加清晰，与其他算法相比，Wavelet-K-SVD-BF方法在杂波抑制和图像目标清晰度处理方面具有良好的效果;在运行时间上，均值滤波和中值滤波法处理的时间最长，小波阈值变换的时间最短，Wavelet-K-SVD-BF方法处理时间较为适中。从图像直观效果以及处理后图像熵值、信噪比和运行时间进行综合对比，本文Wavelet-K-SVD-BF方法在探地雷达杂波抑制和噪声处理上都取得了较好的效果。

4 结 语

针对探地雷达在工程勘测时得到的原始图像存在大量噪声和地表直达波的干扰，进而造成无法有效提取目标信号的问题，提出了一种基于小波变换和K-SVD的自适应双边滤波方法。此方法将小波变换和K-SVD算法相结合，利用奇异值分解提高对雷达图像的整体分析能力，同时将小波变换对数据时频分析的能力相结合，有效提高了杂波抑制效果，并且加入自适应双边滤波，平滑滤波的同时大量保留了图像的边缘信息和细节特征，增加了目标图像的清晰度。实验结果表明，与其他方法相比，本文方法具有更好的峰值信噪比和图像熵，杂波抑制效果更理想，证明了该方法在原始雷达信号处理上具有更好的降噪去杂波的效果，也证明了方法的鲁棒性和可行性。虽然本文方法可以在提取探地雷达目标信号时有效去除了杂波和噪声的干扰，但是在保留图像的细节特征方面还存在不足。为了更好地研究地下目标信息，未来将进一步改进算法，实现在对目标信号去杂波、去噪声的同时保留更多的细节特征的目标。

参考文献/References：

[1] GIOVANNESCHI F，MISHRA K V，GONZALEZ-HUICI M A，et al.Dictionary learning for adaptive GPR landmine classification[J].IEEE Transactions on Geoscience and Remote Sensing，2019，57（12）：10036-10055.

[2] 纪丽静，施养杭.探地雷达在工程检测中的应用与发展[J].无损检测，2010，32（3）：201-204.

JI Lijing，SHI Yanghang.Review in engineering applications and development of ground-penetrating radar detection[J].Nondestructive Testing，2010，32（3）：201-204.

[3] 秦镇，张恩泽，吴海波.基于探地雷达的城市道路地下空洞探测研究[J].安徽理工大学学报（自然科学版），2018，38（5）：70-73.

QIN Zhen，ZHANG Enze，WU Haibo.Research on underground cavity detection of urban road based on ground penetrating radar[J].Journal of Anhui University of Science and Technology （Natural Science），2018，38（5）：70-73.

[4] 郑俊.一种用于探地雷达空洞数据预处理的Shearlet自适应阈值函数滤波方法[D].西安：长安大学，2019.

ZHENG Jun.An Adaptive Threshold Shearlet Filtering Function for GPR Voids Data Preprocessing[D].Xian：Changan University，2019.

[5] ZENG Z F，WU F S，HUANG L，et al.The adaptive chirplet transform and its application in GPR target detection[J].Applied Geophysics，2009，6（2）：192-200.

[6] WANG X，LIU H B，ZHENG M.A clutter suppression method of ground penetrating radar for detecting shallow surface target[C]//IET International Radar Conference.：，2016：1-4.

[7] VERMA P K，GAIKWAD A N，SINGH D，et al.Analysis of clutter reduction techniques for through wall imaging in uwb range[J].Progress in Electromagnetics Research B，2009，17：29-48.

[8] 陳瑞鼎，鹿琪，单子涵，等.基于卡尔曼滤波的超宽带穿墙雷达移动目标检测[J].地球物理学进展，2017，32（4）：1758-1763.

CHEN Ruiding，QI Lu，SHAN Zihan，et al.Moving target detection with the UWB through-wall radar based on Kalman filter[J].Progress in Geophysics，2017，32（4）：1758-1763.

[9] 王超，沈斐敏.一维HHT变换在探地雷达数据处理中的应用[J].工程地质学报，2015，23（2）：328-334.

WANG Chao，SHEN Feimin.Application of one-dimensional HHT transform to GPR data processing[J].Journal of Engineering Geology，2015，23（2）：328-334.

[10]SANTOS V，TEIXEIRA F.Application of time-reversal-based processing techniques to enhance detection of GPR targets[J].Journal of Applied Geophysics，2017，146：80-94.

[11]杨建功，江泉，刘文颖，等.探地雷达在堤防检测中的应用[J].水电能源科学，2012，30（1）：156-158.

YANG Jiangong，JIANG Quan，LIU Wenying，et al.Application of GPR in quality detection of dikes[J].International Journal Hydroelectric Energy，2012，30（1）：156-158.

[12]周克良，刘亚亚.新阈值小波变换的心音去噪[J].计算机工程与设计，2020，41（9）：2476-2481.

ZHOU Keliang，LIU Yaya.New threshold wavelet transform for heart sound denoising[J].Computer Engineering and Design，2020，41（9）：2476-2481.

[13]刘成明，王德利，王通，等.基于Shearlet变换的地震随机噪声压制[J].石油学报，2014，35（4）：692-699.

LIU Chengming，WANG Deli，WANG Tong，et al.Random seismic noise attention based on the Shearlet transform[J].Acta Petrolei Sinica，2014，35（4）：692-699.

[14]李永军，马立元，段永刚.一种改进的小波阈值去噪算法及其仿真分析[J].自动化仪表，2011，32（7）：21-24.

LI Yongjun，MA Liyuan，DUAN Yonggang.An improved wavelet threshold denoising algorithm and its simulation analysis[J].Process Automation Instrumentation，2011，32（7）：21-24.

[15]程卓越.一种高效的自适应双边滤波方法[J].数字技术与应用，2019，37（10）：121-123.

CHENG Zhuoyue.An efficient adaptive bilateral filtering method[J].Digital Technology and Application，2019，37（10）：121-123.

[16]郭贝贝，易三莉，贺建峰，等.改进的非局部均值滤波算法[J].计算机工程，2016，42（7）：227-231.

GUO Beibei，YI Sanli，HE Jianfeng，et al.Improved non-local means filtering algorithm[J].Computer Engineering，2016，42（7）：227-231.

[17]覃尧，黄春琳，陆珉，等.基于小波变换与主成分分析的探地雷达自适应杂波抑制方法研究[J].雷达学报，2015，4（4）：445-451.

QIN Yao，HUANG Chunlin，LU Min，et al.Adaptive clutter reduction based on wavelet transform and principal component analysis for ground penetrating radar[J].Journal of Radars，2015，4（4）：445-451.

[18]GIANNOPOULOS A.Modelling ground penetrating radar by GprMax[J].Construction and Building Materials，2005，19（10）：755-762.

[19]王康，杨瑞祺，杨依忠，等.基于二阶牛顿插值的图像自适应缩放设计及实现[J].计算机应用与软件，2020，37（9）：126-138.

WANG Kang，YANG Ruiqi，YANG Yizhong，et al.Design and implementation of Image adaptive scaling based on second order Newtonian interpolation[J].Computer Applications and Software，2020，37（9）：126-138.

[20]郑俊，刘婉萍，马念茹，等.探地雷达信号杂波抑制[J].电子测量技术，2018，41（21）：63-66.

ZHENG Jun，LIU Wanping，MA Nianru，et al.Clutter reduction for ground penetrating radar signal[J].Electronic Measurement Technology，2018，41（21）：63-66.