基于低秩稀疏分解的GPR杂波抑制方法

陈 诚, 宋晓骥, 何志华, 刘 涛, 曹来保, 粟 毅

(国防科技大学电子科学学院, 湖南 长沙 410073)

0 引 言

探地雷达(ground penetrating radar,GPR)利用电磁波对介质的穿透性,可以实现对地下目标的成像探测,如地雷、管线等[1-3]。相比基于声学、光学、化学和核物理等的技术,基于 GPR 的探测技术得益于电磁波对介质的穿透能力、非接触的测量方式和有效获取目标散射特性的能力,能够更好地满足地下探测特别是地雷探测对准确性和安全性的要求,因而受到了研究者的密切关注[4-5]。

在GPR探测时,接收到的雷达回波中通常不仅包含目标响应,还包含了地表反射回波和天线间直接耦合等。尤其是地表反射回波,其往往远强于目标信号,因此在雷达图像中目标响应往往被掩盖于强杂波之中,导致目标检测十分困难[6]。杂波抑制方法研究是GPR探测领域中的一个重要方向,传统杂波抑制方法如平均值对消(mean subtraction, MS)法[7]、时间门技术[8]、奇异值分解(singular value decomposition, SVD)法[9-11]等,在地面平整均匀且参数已知的情况下杂波抑制效果较好,但当实际地表粗糙起伏时,往往由于环境复杂,难以准确估计地表回波时延或是杂波子空间,导致抑制效果变差甚至失效。形态学成分分析(morphological component analysis,MCA)[12-14]通过结合图像的稀疏表示理论,利用信号组成成分的形态差异将图像分解为目标和杂波,其性能依赖于稀疏表示字典的构建。

低秩稀疏分解(low rank and sparse decomposition,LRSD)是指将探测数据分解为一个低秩矩阵与一个稀疏矩阵的叠加,在模式识别、图像处理、雷达探测领域均有应用,主要用于提取低秩背景或稀疏目标[15-17]。近年来,稳健主成分分析(robust principle component analysis,RPCA)将LRSD理论应用于雷达穿透成像领域以实现杂波抑制[16-17],并展示了其应用潜力。由于GPR在探测时,目标响应往往被地表强反射淹没,因此杂波以地表背景为主,且由于目标(主要包括地雷等)数量有限,通常呈稀疏分布,因此在GPR领域,目标图像可以被视为稀疏成分,而杂波可以表征为低秩矩阵[17]。由于矩阵的秩是非凸非连续的,直接计算其解,将原始图像分解为目标和杂波信号是非确定性多项式(non-deterministic polynomial,NP)困难问题。RPCA通过求解核范数最小化,将其松弛为凸优化问题,但在核范数最小化时引入的SVD迭代导致计算量增大[11],降低了求解效率。此外,将核范数作为正则化项时,大奇异值在优化过程中起决定性作用,容易造成矩阵信息损失[18-19]。文献[20]提出加权核范数最小化(weighted nuclear norm minimization,WNNM)的思想,即用小权值约束大奇异值,用大权值约束小奇异值,有利于更好地恢复出低秩矩阵,但该方法仍存在计算效率低的问题[21]。文献[22]提出随机SVD算法,将高维矩阵投影到低维后再进行分解运算,能够在提高计算效率的同时保持精度。

针对上述问题,本文提出一种基于低秩稀疏分解的杂波抑制方法。该方法将加权核范数正则化项引入RPCA方法,通过控制核范数的权重,有利于精确估计出低秩矩阵,并利用随机SVD代替核范数计算时的奇异值分解以减少计算量,结合交替方向乘子(alternating direction method of multipliers,ADMM)法,建立带有增广拉格朗日项的优化方程求解目标函数,以进一步提高求解效率[23-26]。本文实验应用仿真与实测数据,通过与多种方法进行对比,定量和定性地分析了本文方法的性能优势,表明所提方法能够在GPR应用中有效实现杂波抑制。

1 LRSD及RPCA理论

根据低秩稀疏分解理论,对于GPR系统,可以将获取到的m×n维(假设m≥n)雷达数据矩阵D进行分解[27],即

D=L+S+N

(1)

式中:L表示含少数非零奇异值的低秩矩阵;S为含少数非零项的稀疏矩阵;N表示随机噪声,其Frobenius范数满足

(2)

RPCA方法将上述分解转换为如下易解的凸优化问题[15]:

(3)

2 基于LRSD的杂波抑制方法

鉴于核范数在低秩稀疏分解中的局限性,将加权核范数引入到RPCA方法中,通过控制核范数的权重来约束各奇异值,提高求解的精确性。矩阵L的加权核范数定义为

(4)

式中:w=[w1,w2,…,wn]T;σi(L)为矩阵L的第i个奇异值;wi表示该奇异值对应的权重。

因此,式(3)中的问题可转换为

(5)

在本算法中,定义加权核范数的权重值为

wi=μ/[σi(D)+ε]

(6)

式中,μ为正则化参数;ε为大于0的极小值。可以看出,当ε→0时,有如下关系:

表明定义的加权核范数能有效逼近数据矩阵的秩。

为了提高算法效率,在进行矩阵分解时,引入基于随机SVD算法[28]过程的随机SVD算法,即对m×n维雷达数据矩阵D,生成n×l维高斯随机矩阵X(l远小于n)。根据图1所示流程生成样本矩阵并进行正交三角分解,通过正交变换Q构建远小于原数据矩阵D的l×n维矩阵E。由于正交变换的稳定性,在变换过程中能够保持原矩阵的各类性质,因此对E进行奇异值分解能够在保持运算精度的条件下提高运算效率。

图1 随机SVD流程图Fig.1 Flowchart of randomized SVD

根据ADMM算法思想求解式(5),引入等式约束,得到增广拉格朗日函数为

(7)

式中:U表示拉格朗日乘子矩阵;(·)H表示对矩阵进行共轭转置;β为惩罚系数。通过引入拉格朗日乘子项为优化函数增加等式约束,能够将目标函数的解约束于可行域内,提升算法的求解效率。

在求解式(7)中的3个变量时,可通过分别固定其他两个变量、更新其中一个变量的方式,具体迭代过程如下:

(8)

(9)

Uk+1=Uk+β(D-Lk+1-Sk+1)

(10)

式中:k表示雷达数据第k次的迭代值。

综上,可提出本文基于LRSD的杂波抑制方法流程,如算法1所示。

算法 1 基于LRSD的杂波抑制方法输入 原始雷达矩阵D,正则化系数λ,β,最大迭代次数K,精度阈值η。输出 低秩矩阵L和稀疏矩阵S步骤 1 初始化L0=S0=U0=0。步骤 2 更新低秩杂波矩阵Lk+1=argminLLw,*+(β/2)D-L-Sk+(Uk/β)2F计算核范数最小化时,通过随机SVD算法对待分解矩阵进行降维处理。步骤 3 更新稀疏目标矩阵Sk+1=argminSλS1+(β/2)D-Lk+1-S+(Uk/β)2F步骤 4 更新拉格朗日乘子矩阵Uk+1=Uk+β(D-Lk+1-Sk+1)步骤 5 若不满足停止准则,继续步骤2～步骤4,若满足停止准则,跳出循环,得到输出低秩杂波矩阵L和稀疏目标图像矩阵S。

3 实验结果和分析

为了验证所提低秩稀疏分解方法的有效性,本节利用仿真与实测数据来验证其杂波抑制性能,并通过与MS、SVD、MCA和RPCA这4种方法处理后,将成像结果进行对比,证明算法在成像质量和速度上的优势。

3.1 仿真实验设计

仿真通过GPRMax[29]电磁仿真软件获取场景数据,实验场景为表面略有起伏(不超过10 mm)的土壤层,相对介电常数为5,土壤内部埋藏有两个相对介电常数为2.9的塑料圆柱,用于模拟埋藏在土壤中的地雷,如图2所示。左侧地雷直径为60 mm,高为40 mm,距离地面90 mm,右侧地雷直径为60 mm,高为40 mm,距离地面100 mm,两地雷间隔为25 mm。天线扫描面距离介质表面150 mm,扫描步进为10 mm,发射信号为雷克子波,频率为2.2 GHz。

图2 仿真实验场景示意图Fig.2 Schematic diagram of simulation test

3.2 实测实验设计

实验1通过如图3所示的雷达探测系统进行,实验目标为若干72式反步兵地雷,地雷直径为79 mm,高度为35 mm,等间隔埋藏于沙地中,埋藏深度约为100 mm,每两个地雷之间间隔为50 mm。天线距离地面约200 mm,扫描频率范围为0.75～4.75 GHz,空间采样间隔为20 mm。

图3 实验1场景图Fig.3 Schematic diagram of Experiment 1

实验2借鉴了乔治亚理工学院出版的雷达数据[30-31],实验频率范围为60 MHz～8.06 GHz,频段内共计401个频点,频点在频率范围内均匀分布,天线相位中心距离地面278 mm,空间采样间隔为20 mm。实验场景及目标布设情况如图4所示,这里选取x=45 cm处的数据(即图4中红线标注位置)进行数据处理及对比实验。对上述数据进行成像,在此基础上分别利用MS、SVD、MCA、RPCA与本文所提方法进行处理,处理结果如图5～图7所示。

图4 实验2场景图Fig.4 Schematic diagram of Experiment 2

图5 仿真实验结果图Fig.5 Simulation experiment results

图6 实验1结果图Fig.6 Experiment 1 imaging results

图7 实验2结果图Fig.7 Experiment 2 imaging results

3.3 实验结果与分析

图5所示为对仿真数据使用各算法处理前后的成像结果。其中,图5(a)为直接成像结果,可以看出,土壤回波对目标回波干扰严重。图5(c)和图5(d)表明,由于土壤表层的起伏,平均值对消法与SVD方法在抑制杂波时几近失效。从图5(e)可以看出,MCA方法尽管能够有效增强图像中的目标,但其构建字典的特性使得目标几何结构发生变化。图5(b)和图5(f)中目标图像明显,表明低秩稀疏分解方法能够在保留目标完整性的前提下有效抑制成像结果中的杂波,且图5(b)中残留的杂波明显弱于图5(f),表明本文方法抑制杂波的效果优于RPCA方法。

图6为对实验1数据处理前后的结果。其中,图6(a)为处理前的成像结果,从中可以看到明显的强地表回波,埋藏目标无法显现;从图6(b)～图6(f)可以看出,几种算法均能对成像结果中的杂波进行抑制。但从杂波抑制效果来看,本文方法效果最好,处理结果中残留的杂波相对于目标已经十分微弱,而平均值对消法、SVD方法与RPCA方法处理结果中仍存在着明显杂波;MCA方法处理后视觉效果较好,但在处理结果中引入了大量规则杂波。图7为实验2结果图,从中可以看出原始图像分别经过平均值对消法、MCA方法、SVD方法及两种低秩稀疏分解方法处理后,目标信息逐渐增强。对比图7与图4中的目标分布,RPCA方法处理结果中存在虚警目标,而本文方法能更好地反映出目标的真实分布情况。图8为对3组数据使用本文方法分别计算出的低秩杂波图像。

图8 本文方法处理各实验数据后的低秩杂波图像Fig.8 Low-rank clutter images with experimental data processed by the proposed method

为了量化评估各方法的细节表现,表1记录了各方法处理前后结果图的目标信杂比(signal-to-clutter ratio, SCR),即图像目标区域平均功率与杂波区域平均功率的比值,其定义为

表1 不同方法处理后成像结果SCR参数Table 1 SCR of imaging results processed by different methods dB

(11)

式中:Rt和Rc分别表示目标和杂波区域;Nt和Nc分别表示目标和杂波区域内的像素点个数;I(p)和I(q)分别表示目标和杂波区域各点的像素值大小。

从表1中的SCR评价指标来看,采用低秩稀疏分解方法处理后,图像的SCR明显高于其他方法,且本文方法SCR也始终高于RPCA方法,证实了所提方法抑制杂波的有效性。

为了进一步评估本文方法的运算效率,表2列出了各方法的处理时间及计算复杂度(设图像维度为m×n,迭代次数为r。由于MCA方法参数初始化过程复杂,未评估其计算复杂度)。由表2可以看出,对于各实验中的单幅图像,所提方法处理时间均不足1 s,且运算效率与RPCA方法相比提升了5倍以上。综上分析,本文方法能够兼顾精度与效率,在提升杂波抑制效果的同时,提高了运算效率。

表2 运算时间和计算复杂度对比Table 2 Comparison of running time and computation complexity

4 结 论

本文提出了一种改进的低秩稀疏分解杂波抑制算法。该方法在RPCA方法的基础上,引入加权核范数正则化项与随机SVD计算方法,并结合ADMM方法求解目标函数,实现低秩稀疏矩阵的准确、高效分解。多组数据的实验结果验证了本文方法抑制杂波的有效性与优越性,且证实本文算法运算效率优于RPCA方法。但当GPR探测仅关注特定目标时,该方法无法排除场景中由其他杂波物体造成的虚警,因此结合目标散射机理研究目标与杂波物体的特性差异将是后续研究工作的重点。