李家强, 陈德昌, 陈金立, 朱艳萍
(1. 南京信息工程大学电子与信息工程学院, 江苏南京 210044;2. 南京信息工程大学气象灾害预报预警与评估协同创新中心, 江苏南京 210044)
穿墙成像雷达(through the wall radar imaging)通过发射电磁波对墙体后静止或运动的目标进行探测与成像,广泛地应用于救援任务、领土安全、城市监控、军事侦察等方面[1-2]。在穿墙雷达成像过程中,很多学者对前墙的强杂波进行了深入的研究,并对其进行了有效抑制。Solimene等提出一种基于熵值的算法[3],能够有效地抑制前墙杂波,该算法是利用墙体与目标的差异来设置适宜的门限值,滤除前墙杂波。然而由于实际中侧墙的存在,接收的雷达信号不仅包括目标与侧墙间多径传播产生的回波信号,也包括侧墙与后墙之间反射产生的强杂波信号,这必将导致多径虚像与墙体杂波的产生,使得目标定位失准、真假目标难辨,降低穿墙雷达检测与成像性能,因此如何降低多径虚像与强杂波的干扰是很有意义的。
Setlur等[4]采用高斯加权函数的方法抑制多径虚像,该方法首先利用最小二乘法来定位虚像的位置,然后对于目标与虚像的位置分别设定不同的加权函数,使目标的像素值保持不变的同时将虚像的像素值置零。文献[5]提出用子孔径成像的方法来抑制多径虚像,分别找出不能接收左右两面侧墙体多径回波的子孔径,通过后向投影算法形成两幅子图像,最后再与全孔径成像图连乘,根据这三幅图像多径虚像的差异性来去掉左右墙的两个多径虚像。上述方法在虚像抑制之前都通过基于先验的背景相消[6]的方法消除了前后墙体杂波以及侧墙与后墙反射形成的强杂波,然而在实际工程应用中是无法直接测量墙体后的背景环境,所以对封闭墙体室内目标成像前预先滤除墙体杂波仍是穿墙雷达一个值得研究的方向。其中经典的去除墙体杂波的方法是奇异值分解法(Singular Value Decomposition,SVD)[7],其主要通过墙体杂波与目标回波信号能量之间的大小差异来确定墙体子空间,进而通过墙体子空间投影方法抑制墙体杂波,然而当侧墙存在时,侧墙与后墙反射产生的墙体杂波能量与目标接近且分布不均匀,因此难以滤除墙体杂波。
为了实现更为实用的方法对墙体后的目标进行检测,获得精度更高的成像,本文提出一种基于伽马校正的子孔径图像相乘的方法,该方法无需对墙体场景进行背景相消也能滤除虚像与墙体杂波。本文通过仿真实验模拟一个无目标的室内成像结果图分析得出两面侧墙与后墙两个强杂波产生与目标多径无关,然后基于子孔径成像相乘原理,抑制两面侧墙的两个虚像与前后墙体杂波。同时采用伽马校正的方法,根据每幅图像对比度因子的比例关系,得到每幅子图像指数权值,并对指数加权后的子图像进行连乘,以消除侧墙与后墙反射产生的强杂波,最后用阈值法滤除后墙产生的虚像残余。
本文基于FDTD的电磁仿真软件GprMax2D/3D[8],对穿墙雷达场景进行建模,探测模型如图1所示,设墙体由单层均匀介质组成,厚度为d,相对介电常数为εr,电导率为σ。目标为一理想的球状电导体,它的半径为r,圆心与前墙体的垂直距离为δ,全向天线与前墙体之间的垂直距离为h,沿着平行于墙体的测量线等间距扫描N次,发射天线发射超宽带窄脉冲信号,接收天线在相同的位置接收回波信号。发射信号为Ricker子波,即一阶高斯脉冲信号取负,其表达式为
(1)
天线接收到的回波信号模型表示为
e(t)=et(t)+ea(t)+ew(t)+em(t)
(2)
式中,et(t)为目标回波信号,ea(t)为天线失配噪声,ew(t)为前后墙体杂波信号与侧墙和后墙之间反射产生的强杂波,em(t)为由目标多径产生的多径回波信号。在本文中,假定全向天线为理想天线,因此不存在失配噪声,而式(2)中后两项是需要滤除和抑制掉的墙体间杂波和多径虚像信号,即为了得到干净的目标信号et(t),必须通过算法滤除墙体杂波信号ew(t)与目标多径回波信号em(t)。
本文将成像区域约束在封闭的墙体中,如图2所示,侧墙与前墙的长度分别为D1和D2,墙体的厚度与介电常数分别为d和εr,目标位于P(x,y),每根天线都是单发单收型,且沿着前墙平行的直线排列,每相邻两根天线的距离相等,共N根天线。第n根天线的位置是(xn,yn),n=1,2,…,N。
图2 多径传播模型
(3)
式中,c为电磁波的传播速度。设s(t)为发射信号,则得到目标回波信号与目标多径回波信号分别为
(4)
(5)
式中,Γpn为与传播路径path-p相关的目标散射参数。式(4)代表发射信号由第n根天线沿路径A到目标,再由目标沿此路径返回到天线,式(5)代表发射信号沿路径A到目标,再由路径B,C,D返回到天线。
本文提出的算法由3个步骤完成:第一步,找到两个合适的子孔径,即既要找到尽可能多的天线保证成像分辨率,又要确保只能接收单边侧墙的虚像;第二步,用两个子孔径图像与全孔径图像连乘来抑制侧墙的两个虚像与前后墙体杂波;第三步,根据伽马校正的方法对每幅图像设置一个指数权值并连乘以消除侧墙与后墙反射产生的强杂波。
如图2所示,以左墙为例:当作P点与左墙镜像对称点时,此对称点与前墙左端连线的延长线与天线有一个交点,这个交点距离天线左端有一段距离,根据电磁波反射的几何关系可知在这段距离中排放天线无法接收左墙的目标多径回波,然而在实际测量中无法预先探知目标的位置。选取子孔径的方法如下:天线从左端第1根开始排放天线,并进行BP(back projection)后向投影算法,然后依次累加一根天线,当第k根天线开始能接收到左墙虚像时,则第1根到第k-1根天线为子孔径1,同理可以得到子孔径2。
假设有N幅子图像,都是由BP后向投影成像算法产生,每个子图像的行列像素都相同,第i幅子图像表示为Pi,对其进行归一化处理:
Ii(x,y)=Pi(x,y)/max∀(x,y){Pi(x,y)}
(6)
式中,Ii(x,y)为图像Ii的第x行第y列的归一化像素值。将N幅图像进行相乘融合[9],即为归一化幅度图像的逐像素乘积,设融合后的图像像素值为If(x,y)。
If(x,y)=I1(x,y)*I2(x,y)…*IN(x,y)
(7)
鉴于本文多径回波模型取N=3,其中I1是全孔径成像图。I2无法接收左墙虚像,即子孔径1成像图,I3无法接收右墙虚像,即子孔径2成像图。在每个子图像中,杂波分布是不同的,左右两面侧墙目标虚像的位置也是相异的,但是目标的位置是相同不变的,所以用子孔径图像相乘可以去除侧墙的两个虚像与前后墙体杂波。由于两个侧墙与后墙产生的强杂波能量与目标接近,所以上述的方法处理后仍然有杂波残余,可进一步采用伽马校正[10]的方法进行消除。
在图像处理中,对图像进行伽马校正(又称幂律变换)的修订,能够使得图像信息的对比效果更加明显,伽马校正的基本公式为
s=c(r+ε)g
(8)
式中:s表示图像的输出数据;c表示常量系数,通常可以将其设为1;r表示输入图像的输入数据;ε表示为了有一个可适度的输出值而加入的一个偏移量。然而在一般情况下不在计算范围内,因此可以省略偏移量ε,于是公式可以简化为
s=crg
(9)
其中,g决定了像素变换的方式,g=1时代表线性变换,不会影响图像对比度;g≠1时代表非线性变换,会导致图像对比度发生改变。将常量系数c设为1,根据指数与对数的关系可知:
(10)
所以当指数权值g越大时,根据对数函数的性质,输入像素值r越大的区域有更宽的输出值范围,而像素值越小的区域则有更窄的输出值范围,所以在穿墙雷达图像中,一般来说目标区域的像素值最大,而其他杂波背景像素值较小,所以用伽马变换可以让穿墙雷达图像中的目标区域更加鲜明。
设有N幅子图像,每个子图像赋予值为gi(i=1,2,3,…,N)的指数权值,并用一个图像评估标准(NSSI)来约束它,把这个标准称为图像对比度[11],图像对比度与指数权值的关系为
(11)
gn=NSSIn/NSSImin
(12)
在所有子图像中对比度最低的gmin=1,其他子图像gi值通过比例关系算出,不同的子图像像素被赋予不同的指数权值然后相乘如下:
(13)
Iλ(x,y)=
(14)
式中,λth表示一个预定义的阈值。
为了方便比较虚像与杂波抑制的效果,现定义目标杂波比(Target-to-Clutter Ratio,TCR)如下[12]:
(15)
式中,I(·)为成像点所对应的归一化幅度值,At和Ac分别为目标区域和杂波区域(虚像也作杂波处理),Nt和Nc为各区域所对应的成像点数目。
根据图1所示的系统模型,建立仿真场景:全向天线阵列沿平行于墙体的测量线等间距进行26次扫描,且每根天线与前墙存在一段距离,天线的横向扫描范围为0.1~2.1 m。墙体是均匀介质,其厚度d=0.2 m,相对介电常数εr=6.4。目标为球状的理想电导体,其球心距离前墙体δ=1 m,半径r=0.1 m。通过GprMax2D/3D软件获取各个测量位置的回波信号,使用BP后向投影成像算法进行全孔径成像,结果如图3所示。
图3 全孔径成像图
图3中存在8个重点区域,如箭头所指。除了目标外,其余分别为虚像与墙体杂波。根据虚像定位的基本理论可知M1,M2,M3分别由左墙、右墙、后墙的目标多径回波产生的虚像。通过空场景成像图4所示,C1与C2是由侧墙与后墙的反射产生的杂波。
基于SVD去除四面墙体模型的墙体杂波时,如图5所示,该算法失效,而图6的基于背景相消法能够有效地滤除墙体杂波,但是该方法需要事先获取室内空场景回波信息,实际情况很难满足。
结合本文设定的模型,通过子孔径选取计算得到子孔径1为第1根到第8根天线,子孔径2为第19根到第26根天线,子孔径1,2的成像结果如图7、图8所示,图7与图8中分别抑制了虚像M1与M2,同时仿真实验表明目标图像始终存在,且位置不发生改变。
图4 空目标成像图
图5 SVD去墙体杂波
图6 背景相消去墙体杂波
图7 子孔径1成像
图8 子孔径2成像
将两个子孔径图像与全孔径图像直接相乘,得到结果如图9所示,虚像M1与M2及前后墙体杂波得到了有效抑制,但是仍然有大量的强杂波残余C1和C2。
图9 子孔径图像相乘
为了解决上述杂波残余的问题,采用了伽马校正的方法,即对每个子图像设置一个指数权值gi,gi的大小计算方法如上文介绍,根据得到的3个子图像的像素矩阵,通过比例关系得出g1=1,g2=4.11,g3=4.03,将三幅图像设置对应指数权值后连乘,并设置阈值λth=0.5得到如图10,从图中可以看出C1,C2及M3得到完全抑制。
图10 伽马校正后用阈值法处理
为了验证本文算法的性能,给出目标位于不同位置时,传统子孔径图像相乘方法与伽马校正之后的信杂比,如表1所示,当目标位置不同时,信杂比增量最小为24 dB。
表1 不同横坐标位置时输出目标图像TCR比较
由于在实际系统中噪声是不可避免的,因此需要在噪声的影响下评估提出算法的性能。图11给出在回波信号中加入不同信噪比时,原始成像、子孔径图像相乘、伽马校正之后目标区域TCR的变化,从图11可知通过伽马校正能提高TCR。图12给出伽马校正算法中3个指数权值随加入噪声的变化,可知g1值恒定不变为1,而g2与g3的值随着SNR增加而增大,这是由于SNR增加,图像对比度也随着增大,由伽马校正的公式可知每个指数权值也会相应变大,当SNR很小时,3个指数权值都接近1。
图11 加入噪声时信杂比的变化
图12 加入噪声时指数权值的变化
为了解决穿墙雷达中目标与侧墙、侧墙与后墙之间的多径虚像和杂波问题,首先对其产生原因进行了回波仿真实验分析。根据试验分析结果,提出基于伽马校正的子孔径图像相乘的方法,即通过子孔径图像相乘的方法消除由侧墙产生的虚像与前后墙体杂波,并以此为基础采用伽马校正的方法,计算出每个子图像的图像对比度。然后根据三幅图像对比度比例关系,将每个子图像设置相应的指数权值,并进行连乘,有效地去除了侧墙与后墙之间的强杂波。从仿真结果可以看出信杂比有了明显的提高,充分说明了本文所提方法对于实际工程中强杂波背景下穿墙雷达多径虚像抑制具有一定的理论指导意义。