0 引言
超宽带(Ultra Wideband,UWB)穿墙雷达成像(Through-the-Wall Radar Imaging,TWRI)是一种通过电磁波对隐藏目标进行检测、定位、成像和特征提取的技术,在人质救援、建筑物布局检测、隐藏目标探测等领域具有广泛的应用前景[1]。现有研究都是将均匀墙作为背景环境,然而当成像环境是周期性结构的空心墙时,电磁波穿过这种墙体会发生Floquet现象,使用现有的均匀墙环境下的成像算法就会产生多径虚假像[2]。
针对上述问题,文献[2]对周期性空心墙后的目标成像问题进行了研究,它强调如果在散射模型中不考虑多径传播现象,就会导致成像结果中目标像的扭曲和失真,因此采用基于玻恩近似的线性化散射模型来解决这个问题,但是从最终结果可以看出这种方法的去噪效果并不理想。文献[3]对周期性空心墙体结构导致的Floquet模式进行了仿真分析,结果证明周期性结构引起的高阶空间谐波是导致穿墙成像中杂波和“鬼影”出现的主要原因。在此基础上,比较了后向投影(Back Projection,BP)成像和传统的稀疏成像方法。BP成像通过将每条多径叠加到目标真实位置来消除Floquet现象,但这种方法数据量大、速度慢、主瓣宽且旁瓣多。而传统稀疏成像虽然精度上有所提高,数据量减小,但是仍然存在速度较慢的问题。不仅如此,上述方法都是针对点目标的成像,考虑到实际情况中更常见的是扩展目标,如果继续使用点目标模型对扩展目标进行成像,将会出现扩展目标边缘轮廓不清晰的问题。
全变分法(Total Variation,TV)是一种基于偏微分方程的变分去噪模型,最初应用在图像稀疏重构中[4]。文献[5]最先将TV法应用到扩展目标雷达成像中,通过增加一阶差分算子D在保留扩展目标边缘轮廓的同时减少旁瓣。文献[6]提出一种目标探测和杂波去除方法,它指出通过反复改变天线位置,直接路径的目标反射会有相同的回波,而间接路径的回波通常表示为随机噪声。因此,该方法可以从各类噪声中分开目标像点。文献[5-6]都能够对扩展目标进行成像,但都是在均匀墙环境下的成像,不能解决空心墙体导致的目标位置偏离问题,同时面临着需要人工设置参数的问题。
综上所述,针对穿过周期墙体的成像场景而言,要在保留点目标和扩展目标特性的同时,解决空心墙体导致的目标位置偏离和需要人工设置参数的问题。本文提出一种结合TV约束的自适应扩展目标稀疏成像方法,首先根据电磁波穿过空心墙体的电磁传播特性将空心墙体等效成三层均匀介质模型,基于贝叶斯推理准则,利用“Integrate-out”方法[7]得到点目标和扩展目标同时存在情况下的最大后验(Maximum a Posteriori,MAP)估计,然后采用交替循环迭代法对各个未知量进行循环优化直至达到最优值。最后通过仿真和实验验证了该算法的有效性。
1 信号模型
如图1所示,墙体相对介电常数为ε1,厚度为h,其中第一、二、三层的厚度分别为h1,h2,h3,墙体周期为d,空心结构之间的材料厚度为a。L个收发共置的天线构成的天线阵列进行接收和发射信号,天线阵列与墙体平行,发射由J个频率
组成的步进频信号。当天线发射的电磁波以β角穿过这种特殊的周期结构墙体时,将会在墙体另一侧的固定角度上产生与入射方向有关的Floquet模式,其第s阶Floquet模式输出角满足条件:
式中,k=2πf/c表示自由空间波数。图1中,天线在每一个频点发射的电磁波都有很多入射平面波,根据式(1)可以知道,每一个入射波都会产生几个Floquet透射波,这就会导致计算量的大大增加[3],然而只有直接路径产生的像是真实目标像,由间接路径产生的虚假像与真实目标相比旁瓣更大而能量却较弱,且形状很不规则,通常将这种间接路径导致的杂波表示为噪声[6]。在文献[8]中将这种复杂的周期性空心墙等效成三层均匀介质模型,如图2所示。
图1 电磁传播模型
图2 空心墙的等效三层模型
第二层的介电常数表示为ε2=1-g+g·ε1,其中g=a/d表示填充参数。这样墙体内传播时延就表示为三层均匀介质内时延的叠加,由此可以得到总的传播时延:
式中,xp和yp为像素点位的位置,xl和yl为第l个天线的位置,c为电磁波传播速度,θ≈arctan(|yp-yl|/|xp-xl|)。由此将接收天线接收到的回波信号表示为
式中,η为周围环境和墙体自身产生的噪声,σp为第p个像素点的反射系数,τp,l为第l个天线和第p个像素点之间的双向传播时延,NxNy为成像区域划分的像素点数。对式(3)进行矢量化表示,y∈CJL×1是包括所有JL次测量值y(j,l)的列向量:
式(7)中包括了JL的全部测量数据,可以使用较少的数据并从中恢复出σ。用一个MN×JL,(MN≪JL)的测量矩阵Φ对y进行采样[9],得到一个较稀疏的矩阵:
2 成像算法
对于不同天线位置的接收回波,直接路径产生的目标反射会有相同的回波,而间接路径产生的多径回波类似随机噪声,由于TV法通常只对噪声敏感而平滑信号的TV项通常较小,因此本文采用TV法从各类噪声中分出目标,在去除噪声的同时保留扩展目标的边缘信息[6]。然后采用交替循环迭代的方法[10-12],对各个变量进行循环优化,直至求出最优值。这种方法避免了人工设置参数和复杂的噪声功率求解,简化了求解过程。
2.1 MAP估计
采用MAP方法,假设噪声矢量η服从均值为零、协方差矩阵为βI的复高斯分布,则观测向量
的似然函数[13]表示为
由于散射体分布规律的稀疏性,σ通常满足Laplace分布,并用Jeffrey先验作为Laplace分布的先验参数[7],那么σ的先验信息表示为
式中,第一项表示点目标的散射特性,第二项表示扩展目标的散射特性,α1,α2为Laplace分布的先验参数,表示对散射体特征的增强程度,D为一阶差分算子,[D]i,i=1,[D]i,i+1=-1,1≤i≤NxNy-1,是为了保留扩展目标的边缘特性。Jeffrey先验由Fisher信息量I(·)定义[7]:
将式(11)、式(12)和式(13)代入式(10),采用“Integrate-out”方法[7],可得
式中,Γ(·)表示Gamma分布,基于贝叶斯准则,σ的后验概率表示为
在MAP估计中,通常将p(y)忽略,故σ的M AP估计表示为
为了简化求解,首先进行取对数操作,即
对上式进行整理,去掉常数项,式(16)可以表示为
2.2 交替循环法
采用传统的牛顿法[5]来求解J(σ,β)的最优化问题需要的计算量很大,通常目标函数的Hessian矩阵无法保持正定从而使这些方法失效,同时伴随着正则化参数的选择困难问题。因此,这里使用交替循环迭代的方法求解反射系数σ和参数β。
首先假定得到了第t次迭代的β(t),求出J(σ,β(t))关于σ的导数并令其为零,可求解出σ(t+1),即
式中,Λ1=diag{1/|σ|i},Λ2=diag{1/|Dσ|i},为了避免出现奇异值,在Λ1,Λ2中加入一个极小值κ,表示为
}。式(20)变为
式(25)中,逆运算的求解运算量较大,这里采用文献[8]的共轭梯度最小二乘法(Conjugate Gradient Least Squares,CGLS)对其进行求解,从而实现快速计算。
接下来,对J(σ(t),β)关于β求导,并令其为零,即
2.3 算法实现
算法流程如表1所示,首先采用BP算法得到初始值σ(0),初始值
,最大迭代次数设置为T。循环迭代式(25)和式(27),直到满足收敛条件。很明显,σ和β通过式(25)和式(27)的迭代更新,J(σ,β)会随着t的增加而逐渐减小。文献[11]中已经证明,交替循环迭代算法是收敛的,收敛条件为‖σ(t+1)-σ(t)‖2/‖σ(t)‖2＜δ,这里的δ为很小的正实数。
表1 算法流程
从算法的流程可以知道,每次迭代过程都需要对σ和β进行更新,其中对σ的求解过程需要进行逆运算,求解式(25)的逆运算计算复杂度为O((NxNy)3)。为了减小计算量,首先采用文献[12]的“Integrate-out”方法将σ无参数化,在这个过程中,充分利用了目标的先验信息,提高了成像结果的稳健性。然后采用文献[11]中的CGLS算法求式(25)的逆运算,CGLS算法在每次迭代只包括两个矩阵的乘积,因此每次迭代的计算复杂度仅为O(MN×NxNy),而且文献中提到CGLS算法的迭代次数只与接收信号的信噪比有关,这就在实现快速收敛的同时大大减小了计算量。
3 仿真及实验结果分析
3.1 仿真结果分析
使用XFDTD电磁仿真软件进行仿真。非均匀空心墙的介电常数ε1=5,电导率δ=0.01 S/m,厚度为0.2 m,第一、二、三层的厚度分别为0.04,0.12和0.04 m,墙体周期为0.15 m,空心结构之间的材料厚度为0.04 m。由21个收发共置天线组成的天线阵,与墙体平行且距离墙体0.5 m,阵元之间间隔0.05 m,发射频率2 GHz脉宽1 ns的调制高斯脉冲信号,提取频段1.5～2.5 GHz的DFT数据,共101个频点。图3给出了XFDTD仿真模型,图3(a)为长0.4 m、宽0.2 m的矩形单目标体,图3(b)为两个矩形目标体,各目标与墙体的距离均为1 m。
图3 XFDTD仿真模型
图4和图5分别给出了单目标体和两个目标体的成像结果,图像中的虚线方框表示目标位置。从图中可以看出,图4(a)和图5(a)的BP成像方法效果较差,存在大量的杂波,很难识别出目标位置,扩展目标轮廓不能很好保存。图4(b)和图5(b)采用传统的稀疏成像方法,杂波减少,但仍然存在多径虚假像。而图4(c)和图5(c)采用本文的方法成像相比其他两种效果明显改善,虚假像得到消除,目标准确且扩展目标轮廓十分清晰,同时本文的优化方法提高了计算速度。
图4 单目标成像结果
图5 两个目标成像结果
为了更直观地看出各算法性能上的差异,对存在一个目标的情况下,在程序运行时间(Program Running Time,PRT)、目标杂波比(Target-to-Clutter Ratio,TCR)、图像的熵(Entopy of the Full Image,ENT)和一阶差分均值(Mean of First Order Difference,MFOD)四个方面进行比较[14],结果如表2所示。PRT为算法的运算时间,PRT越小,运行时间越短。TCR为目标区域最大像素值和杂波区域平均像素值的比值,TCR值越大,表明目标相对背景杂波就越强,背景杂波抑制越充分,图像越清晰。ENT为整幅图的复杂程度,其值越小整幅图杂波越小、越干净。MFOD为分布目标的轮廓特性,MFOD值越小,说明检测区域像素值越接近,区域越平滑。从表中可以看出,BP方法的运算速度最快但背景杂波严重,稀疏方法的运算速度很慢,但是图像相对较清晰,本文方法运算速度介于BP方法和稀疏方法之间,但是其他各方面都有所改善。
表2 PRT、TCR、ENT和MFOD算法比较
3.2 实验结果分析
使用美国GSSI公司的探地雷达SIR2.0搭建穿墙实验场景,如图6(a)所示。使用的是长和宽分别为0.39 m,0.19 m的空心砖,第一、二、三层的厚度分别为0.035,0.12和0.035 m,墙体周期为0.18 m,相对介电常数为5。墙后矩形目标距离墙体1 m左右,目标的长和宽分别为0.7 m和0.2 m。选择2 GHz的喇叭天线,天线平行于墙体移动21次,合成孔径长度1 m,在每个位置测量2次,包括有目标和无目标场景。将采集到的数据波形进行平均处理。图6(b)给出了成像结果,虚线框为目标的真实位置。图中可以看出使用本文方法能保证在得到目标轮廓信息的同时抑制杂波,证明了本文方法的有效性。
图6 穿墙实验结果
4 结束语
本文提出的结合TV约束的扩展目标成像方法,相比较于常规算法解决了空心墙导致的Floquet杂波问题,能够在空心墙环境下进行高分辨成像,而且该算法能有效对扩展目标进行边缘成像。同时使用交替循环法进行迭代大大缩减了计算时间。在今后的工作中,将对该方法进行完善使其能应用到更复杂的场景中。
参考文献:
[1]YEKTAKHAH B,SARABANDI K.All-Directions Through-the-Wall Radar Imaging Using a Small Number of Moving Transceivers[J].IEEE Trans on Geoscience and Remote Sensing,2016,54(11):6415-6428.
[2]GENNARELLI G,SOLDOVIERI F.Radar Imaging Through Cinderblock Walls:Achievable Performance by a Model-Corrected Linear Inverse Scattering Approach[J].IEEE Trans on Geoscience and Remote Sensing,2014,52(10):6738-6749.
[3]CHANG P C.Physics-Based Inverse Processing and Multipath Exploitation for Through-Wall Radar Imaging[D].Columbus:The Ohio State University,2011.
[4]ZHAO Yao,LIU Jianguo,ZHANG Bingchen,et al.Adaptive Total Variation Regularization Based SAR Image Despeckling and Despeckling Evaluation Index[J].IEEE Trans on Geoscience and Remote Sensing,2015,53(5):2765-2774.
[5]SUN Bin,JIN Haibo.Feature Enhanced Synthetic Aperture Radar Image Formation[C]∥8th International Conference on Electronic Measurement and Instruments,Xi’an:IEEE,2007:2803-2809.
[6]HANDA H,MANSOUR H,LIU D,et al.Extended Target Localization with Total-Variation Denoising in Through-the-Wall-Imaging[C]∥IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing,Cancun:IEEE,2015:445-448.
[7]OGWORONJO H C,ANDERSON J M M,NGUYEN L.A Parameter-Free MAP Image Reconstruction Algorithm for Impulse-Based UWB Ground Penetrating Radar[C]∥IEEE Global Conference on Signal and Information Processing,Orlando,FL:IEEE,2015:742-746.
[8]BURKHOLDER R J,MARHEFKA R J,VOLAKIS J L.Radar Imaging Through Cinder Block Walls and Other Periodic Structures[C]∥IEEE Antennas and Propagation Society International Symposium,San Diego,CA:IEEE,2008:1-4.
[9]LEIGSNERING M,AHMAD F,AMIN M,et al.Multipath Exploitation in Through-the-Wall Radar Imaging Using Sparse Reconstruction[J].IEEE Trans on Aerospace and Electronic Systems,2014,50(2):920-939.
[10]LIM D,XU L,GIANELLI C,et al.Time-and Frequency-Domain MIMO FLGPR Imaging[C]∥IEEE Radar Conference,Arlington,VA:IEEE,2015:1305-1310.
[11]TAN X,ROBERTS W,LI J,et al.Sparse Learning via Iterative Minimization with Application to MIMO Radar Imaging[J].IEEE Trans on Signal Processing,2011,59(3):1088-1101.
[12]GLENTIS G O,ZHAO K,JAKOBSSON A,et al.Non-Parametric High-Resolution SAR Imaging[J].IEEE Trans on Signal Processing,2013,61(7):1614-1624.
[13]LIU Hongchao,JIU Bo,LIU Hongwei,et al.An Adaptive ISAR Imaging Method Based on Evidence Framework[J].IEEE Geoscience and Remote Sensing Letters,2014,11(6):1031-1035.
[14]晋良念,申文婷,钱玉彬,等.组合字典下超宽带穿墙雷达自适应稀疏成像方法[J].电子与信息学报,2016,38(5):1047-1054.JIN Liangnian,SHEN Wenting,QIAN Yubin,et al.Adaptive Sparse Imaging Approach for Ultra-Wideband Through-the-Wall Radar in Combined Dictionaries[J].Journal of Electronics&Information Technology,2016,38(5):1047-1054.(in Chinese)
