基于非负矩阵分解的GPR高频杂波抑制

苗 翠，原 达，王冬雨，李文生

山东省高校智能信息处理重点实验室（山东工商学院），山东 烟台 264005

探地雷达（GPR）是一种常用的探地无损检测技术[1]，对于地下埋藏物及地下结构的探索具有广泛的应用[2-4]。在采集到的探地雷达信号中，目标的反射频率往往是较强的，杂波可以定义为那些与目标散射特性无关但发生在同一采样时间窗中且具有与目标波长相似的光谱特性的信号[5]。当部分杂波反射更高强度能量并掩盖目标反射波时，称这类杂波为高频杂波，在图像中表现为较高的像素强度。其来源可能是由发射天线和接收天线之间的穿透以及天线和地面之间的多次反射引起的，地面特性阻抗的局部变化以及材料中包含的一组小反射源也会引起高频杂波。尤其对于浅层目标的检测，当目标埋藏物越浅时，使用的探测频率越大，这时所引起的高频干扰也更为复杂。如何从复杂的GPR反射波中去除高频杂波成分保留目标信号，是需要研究的重要课题。

对于探地雷达信号的杂波处理问题，基于主成分分析（principal component analysis，PCA）[6]及鲁棒主成分分析（robust principal component analysis，RPCA）[7]和奇异值分解（singular value decomposition，SVD）[8]等子空间技术往往表现出较好的性能，该类方法将GPR图像分解为对应于杂波和目标的分量，将主分量作为杂波，其余成分作为目标信号来进行处理，但该类方法在GPR图像含有多条杂波或多个目标信号的情况下杂波不能完全去除。近两年提出的形态学方法（morphological component analysis，MCA）[9]将GPR 图像的杂波和目标分量分别用曲波curvelet 和非采样离散小波变换（unsampled discrete wavelet transform，UDWT）字典稀疏表示来进行杂波抑制，虽然视觉效果表现不错，但算法复杂度较高，处理速度不够快。基于多分辨率的方法在GPR 去杂波的处理中也取得了较好的效果，如多尺度双边滤波（multiscale directional bilateral filter，MDBF）[10]借助于杂波和目标在反射波形中呈现出的几何差异进行多方向多尺度的分解。集成经验模式分解（ensemble empirical mode decomposition，EEMD）将探地雷达信号分解为一系列模函数（intrinsic mode functions，IMFs）并计算其置换熵（permutation entropy，PE），通过设置相关阈值来进行噪声和目标的区分，能有效提高目标的分辨率[11]，不过图像质量难以保证。基于径向基函数（radial basis function，RBF）神经网络的方法也同样被用于GPR 图像的杂波滤除，使用零偏移格林函数作为训练数据的期望输出，结果能提高探地雷达的垂直分辨率[12]，但对于目标成分的保留不够清晰。

非负矩阵分解（non-negative matrix factorization，NMF）由Lee和Seung[13]于1999年在自然杂志上提出，该算法与与PCA、ICA（independent component analysis）、SVD 等子空间分解技术不同的是，NMF 中对分解矩阵的非负约束能够使数据进行稀疏表示。分解的结果中不含有负值，使数据更具有可解释性和可使用性。并且由于其简便高效的特点，已经被广泛应用于数据降维、特征选择、图像融合、文本聚类、语音增强等方面[14-17]。对于NMF 在GPR 图像杂波处理的应用中[18]，将GPR 图像自身的像素强度值来进行NMF 分解，去掉第一主成分来进行杂波抑制。而后又提出了鲁棒非负矩阵分解（robust nonnegative matrix factorization，RNMF）[19]，在NMF 基础上加入误差矩阵进行分解，具有更好的实时性。在本文中，基于NMF 提出了概率非负矩阵分解（PNMF）应用于GPR杂波抑制，选取图像的水平梯度作为非负矩阵分解的输入，对于参数后验分布的求取，使用变分贝叶斯方法进行推理。变分贝叶斯是贝叶斯推理中的一种近似推理算法，对于概率模型中参数的求取具有速度快、稳定性好的特点，与传统的NMF 相比，具有更高的信噪比和更好的鲁棒性。

1 基于概率模型的非负矩阵算法

1.1 构建概率NMF模型

在NMF处理GPR图像时，像素强度矩阵F看作是目标成分Ftarget和杂波成分Fclutter的组合结构：

将矩阵F进行非负矩阵分解得到：

其中图像大小为M×N,F∈RM×N,W∈RM×K为基矩阵，H∈RK×N为系数矩阵，且W和H都具有非负约束条件，即s.t.W,H≥0。高频杂波的强度很大，一般选取前K个分量对应于杂波成分，K

其中D(•)表示距离度量函数，用来度量F分解为WH后与F的距离。为了衡量矩阵分解前后的相似性，使用KL 散度（Kullback-Leibler divergence）作为度量指标进行计算，此时的参数W、H的更新规则如下：

在基矩阵W和系数矩阵H的更新中分母1表示的是单位矩阵。上述使用拉格朗日乘子法进行非负矩阵求解时通常具有简单快速的特点，这种非概率的方法容易导致过拟合，不具有稳定性。而贝叶斯方法通过定义满足非负约束的先验分布，再结合实际观察的数据进行后验推理，能大大减少过拟合的现象，提高非负矩阵分解算法的性能。在概率模型下的NMF 中，将待分解的潜在基矩阵W和系数矩阵H看作是两个随机变量，数据F中的值来自于W和H的乘积，以及一些高斯噪声E，如下表达式：

1.2 基于变分贝叶斯分解非负矩阵

在进行变分贝叶斯推理之前，对变量的先验分布P(θ)进行假设，由于PNMF 的非负性约束，设W和H的先验服从指数分布：其中，ε(x;λ)=λexp(-λx)u(x)是指数分布的密度函数，λk>0，u(x)为单位阶跃函数。

噪声方差的先验设为逆伽马分布：

在PNMF 中，参数的后验分布使用变分贝叶斯（variational Bayes，VB）进行推理，通过求解一个近似的后验分布去逼近真实的概率分布。在复杂的多变量参数θ情况下，VB将其分解为一组相互独立的变量θi，表示如下：

其中，q(θ)对应于参数θ的分布。根据文献[20]，得到后验分布：

式（11）中，G(⋅)为伽马分布，密度函数表达式见式（9），TN(⋅)为截断高斯分布，概率密度函数表达式如下：

在截断高斯分布中，x<0 时为密度为0 的高斯分布，其中φ(⋅)是服从N(0,1)分布的累积分布函数。得到W和H以及相关参数的后验分布后，将前K项作为高频杂波对应的分量并进行如下分量表示：

其中K值代表矩阵的秩，通过该值可以得到杂波对应的低秩矩阵。对于高频杂波，前K项成分分量对杂波的贡献最大，K值的选择在实验部分进行了具体论证。

2 实验

2.1 实验数据背景

本文实验分别采用GPR模拟数据和应用环境采集的真实数据进行方法验证分析，并通过不同算法的信噪比及视觉效果来验证本文算法的有效性及鲁棒性。

2.2 模拟GPR数据

模拟数据通过gprMax软件[21]的时域有限差分方法（FDTD）进行模拟获得。在模拟过程中，将铝盘作为埋藏物，然后分别放置于3 种不同特性（干性土壤dry_sand、潮湿土壤damp_sand、湿润土壤wet_sand）的土壤环境下，埋藏物和3 种土壤的介质参数来自于文献[22]，其介电常数和导电率详见表1。

表1 材料的电磁特性Table 1 Electromagnetic properties of materials

获得的GPR 图像如图1 中的A1～A3 所示，对应于铝盘分别埋藏于干性土壤dry_sand、潮湿土壤damp_sand和湿润土壤wet_sand环境下的反射结果，可以从图中清楚的看到，它们的目标反射波形随着环境的改变而发生位置高度和向下开口大小的变化，但是都拥有相同特性的水平杂波。在使用本文算法PNMF进行实验时，首先使用图像自身的像素强度代入计算，得到图1 中B1～B3的处理结果，可以看到该水平杂波并没有完全去除，并且B1、B2 中的双曲波目标顶点部分也被削弱了。当通过图像像素值的水平梯度进行PNMF计算时，结果对应于图1 中的C1～C3，可以看出，处理结果不仅可以很好的保留双曲波信号，还对杂波进行了明显的抑制。

图1 不同PNMF对3幅模拟GPR图像的噪声抑制效果对比Fig.1 Comparison of noise suppression effects of different PNMF on three simulated GPR images

峰值信噪比（PSNR）是图像处理中最为普遍的一种客观评价指标，基于纯净参考图像I与处理后图像Î之间对应像素点(i,j)的误差来判断图像处理的质量好坏。通常通过均方差（MSE）进行定义：

其中M×N为图像I的大小，在模拟GPR数据中，无杂波的参考图像I通过含杂波的目标反射图像和只含杂波的图像作差来得到。在使用概率模型进行非负矩阵分解时，分别使用图像像素强度和水平梯度作为输入，通过PSNR 值的计算得到一组对比柱状图，如图2 所示。可以直观看到在使用梯度进行PNMF 探地雷达图像杂波处理时，得到的PSNR值更高。

图2 像素PNMF和梯度PNMF的PSNR对比Fig.2 PSNR comparison of pixel PNMF and gradient PNMF

由视觉质量和PSNR 对比得到，在PNMF 算法对GPR图像的去杂波过程中，选择其水平梯度作为输入矩阵进行处理效果更佳（后续的PNMF 指的都是梯度PNMF）。

在GPR 图像杂波处理中，杂波成分一般存在于矩阵分解后的前几个主分量中，由此分别选取K=1，2，3，4，5进行处理后PSNR值的比较，结果如表2所示。数据结果显示，干土壤中铝盘反射的GPR 图像在K取2 的时候去杂波后的PSNR值较高，而在潮湿和湿润的土壤环境反射下，K取1时PNMF的峰值信噪比较高。结合

表2 不同K 值下PNMF算法的PSNRTable 2 PSNR of PNMF algorithm with different K valuesdB

2.3 真实GPR数据

子空间技术中GPR图像杂波去除的经验[18]，选取第一主分量作为杂波成分，即K=1。

实验的真实数据是从含有地下管道的表面进行采集的，地质松软平坦，多含细沙碎石。使用GSSI公司提供的探地雷达设备，其中天线的中心频率为400 MHz，探地深度为1.5 m，数据采集模式使用距离采集，在地面做好测距标识之后拖动装置进行扫描，得到图3 的A1，并依次在地下浅层处埋藏空箱和地雷模型得到图3 的A2 和图3 的A3。周围有易产生电磁干扰的电线杆，会对采集的图像造成干扰杂波。

如图3的A1～A3所示。其杂波大体呈水平状、数量多且与目标成分部分重叠的特点。实验时选择了信噪比和视觉质量这两项评价指标来验证PNMF 算法的有效性和鲁棒性。

图像的信噪比是用于比较处理后图像与原图像质量的参数，信噪比的数值越大，图像的质量就越好。信噪比的计算往往需要借助于原始的纯净参考图像，但由于实际测量的GPR图像缺乏参考图像，本实验使用图像的均值与方差之比来计算信噪比进行度量。真实的GPR数据见图3中的A1～A3，分别使用NMF[18]、RNMF[19]和本文的PNMF 在GPR 图像上进行处理，其SNR 结果以折线图的形式呈现在图4。可以看出在GPR 图像的杂波处理中，RNMF 算法的信噪比高于NMF，而PNMF 算法得到的信噪比则比RNMF和NMF算法都要高。

图3 不同算法对3幅真实GPR图像的杂波抑制效果对比Fig.3 Comparison of clutter suppression effects of different algorithms on three real GPR images

图4 NMF、RNMF和PNMF算法处理后的SNRFig.4 SNR value by NMF，RNMF and PNMF algorithms

为了进一步验证本算法在GPR杂波处理中的鲁棒性，分别在真实图像中撒上密度为0.2、0.3和0.5的椒盐噪声，使用SVD、PCA、NMF、RNMF 和本文的PNMF 进行信噪比的计算，结果如表3所示。

表3 加入噪声后不同算法的信噪比值Table 3 SNR of different algorithms after adding noisedB

从表3中的实验数据可以看出，在实验图像中添加不同程度的噪声时，本文算法的去噪效果均要优于其他经典的去杂波算法，信噪比SNR 值相比其他算法都有一定的提高，说明本文算法对不同的噪声环境具有较强的鲁棒性。

同时，满足杂波抑制效果的前提下，算法的时间复杂度也是影响性能的关键因素。因此本文对比了基于NMF 和RNMF 算法运行所需的时间，运行时所使用的计算机操作系统为Windows 7，处理器为Inter Core i7-6700 3.40 GHz，8 GB 的安装内存和64 位的操作系统。运行结果的具体数值见表4。

从表4 中的数据可以看出，对于图3 的A1 和A2，PNMF的运行时间比NMF和RNMF都缩短了4倍左右，图3 的A3 的运行时间与RNMF 相差不大，但总体的运行速度还是有所提高，这是因为PNMF中的变分贝叶斯属于一种近似推理，能够在较低的时间复杂度下获得原问题的近似解。并且面对不同杂波类型时的抑制效果良好。

表4 算法的运行时间Table 4 Running time of algorithms

在图像处理领域，视觉质量是最直观也最重要的评价指标，效果图如图3 所示。针对3 种不同测量环境下的GPR 图像（图3（a）组）分别使用PCA（图3（b）组）、SVD（图3（c）组）传统的杂波去除方法和已经在GPR图像中应用过的NMF（图3（d）组）、RNMF（图3（e）组）方法进行实验，并与本文的PNMF（图3（f）组）进行比较。其中图3（a）组是真实的GPR数据，可以看出3张探地雷达图像的杂波形态不尽相同，但都属于像素强度高于目标反射的高频杂波。如图3的A1中有一条高频水平状杂波在经过SVD、PCA、NMF和RNMF算法处理后水平杂波都得到了一定程度的削弱，但只有本文的PNMF方法在去除杂波的同时保留了区域1处的双曲波，且区域2处的杂波处理也最干净。背景相对纯净的同时很好地保留了目标信号的双曲波。

图3的A2的GPR图像上方含有多条较细的高频水平状杂波，杂波之间相隔紧密，并覆盖双曲波顶部部分信息。在PCA（B2）中杂波的像素强度有所减弱但目标部分也变暗了；SVD（C2）中将第一条最强的杂波去除了，但后半段区域2处的杂波没有本文的PNMF算法效果纯净；在NMF 算法（D2）中，在PCA 的效果基础上对下面几条杂波也进行了部分减弱；该图像中RNMF（见图3 的E2）算法表现最差，没有起到杂波抑制的效果。PNMF对于该图中的几条杂波虽没有完全抑制，但相比之下的背景是较为干净的，且还对目标的边缘细节进行了增强，见图3的F2。

图3的A3中含有3条较粗的高频水平状杂波（方框区域），杂波之间相隔较远且与目标双曲波部分存在大量重叠。可以看到，PCA、SVD、NMF、RNMF对这3条杂波的去除效果都很好，见B3、C3、D3 和E3，但是覆盖在双曲波上的杂波干扰情况并未得到改善（如区域1和区域2）。PNMF（图F3）能准确高效地将杂波与目标分离，从而保留目标的双曲波信息，使处理结果清晰自然。

以上实验结果说明了在GPR图像的高频杂波抑制过程中，使用本文的PNMF 算法能获得较高的信噪比，且视觉上使GPR 图像背景清晰，突出了目标信息。从而验证了本文算法的有效性和鲁棒性。

3 结论

本文提出了一种改进非负矩阵分解的GPR图像杂波抑制方法PNMF，与传统的NMF不同的是，本文采用概率方法进行非负矩阵分解，得到的是分解后矩阵参数的后验分布，然后选取矩阵分解后的第一主成分作为杂波进行去除。对于真实的复杂环境，概率模型往往比非概率方法具有更强的稳定性。PNMF 使用变分贝叶斯作为非负矩阵分解的近似推理方法，追求速度的同时获得了较好的鲁棒性，并且在使用过程中，结合杂波特点提取了探地雷达图像的水平梯度作为输入而不是像素强度。最后的实验对比分析表明，在GPR 图像的杂波抑制问题上，本文算法具有良好的有效性和鲁棒性。后续工作将加入深度学习的知识，结合非负矩阵分解与深度卷积网络对图像数据进行深层次的特征提取，以此提升GPR图像背景去除的纯净度和目标识别的准确率。