1 引言

毫米波是指频率在30~300 GHz的电磁波,由于能够穿透衣服且无电离辐射,特别适用于人体安全检查。在毫米波频段,基于外差混频技术,毫米波收发电路能够实现对反射信号相位的直接测量,因而利用相位信息进行高分辨成像的全息成像技术在毫米波领域受到了广泛的关注^[1-2]。近年来,基于毫米波全息成像技术的人体安检仪发展迅速,在国内外的一些机场已逐步开始使用。

1995年,美国西北太平洋实验室的Collins等^[3]首次将毫米波全息成像技术用于人体安检,并发明了单频线阵扫描的毫米波成像装置。随后,该研究组^[4]将单频扩展到宽带,使成像装置具备了距离向的分辨能力,成像质量大幅度提升。为了实现快速成像,毫米波成像系统常采用线阵加机械扫描的扫描方式(按照扫描轨迹分为平面扫描^[5-6]和圆柱扫描^[7]),基于平面或圆柱面宽带全息数据对人体表面进行高分辨率三维成像。上述线阵加机械扫描的毫米波成像系统采用的是单站或准单站的成像模式,即仅有一个发射天线发射毫米波信号,且仅有一个与发射天线相邻的接收天线对反射信号进行接收(SISO)。除了单站成像模式外,也可以采用多发多收MIMO)模式进行毫米波全息成像^[8-9],即成像系统包含多个发射天线和多个接收天线,每个发射天线发射时可由多个接收天线同时对反射信号进行接收。基于MIMO天线阵列,2010年前后,Rohde & Schwarz公司研制成功了72~80 GHz基于稀疏阵列的毫米波全息成像系统^[10],由于采用全电子扫描,该系统的扫描时间可达157 ms且没有机械运动结构。

毫米波全息成像通过外差混频技术获得包含相位信息的复信号,采用图像重建算法进行数字聚焦。对于单站式成像,由于全息数据均匀排布在平面或圆柱面上,故可以采用快速傅里叶变换(FFT)重建算法直接获取高分辨的毫米波图像^[5]。而对于MIMO成像模式,由于收发数据排布不规则,快速傅里叶变换难以直接应用。近年来,针对毫米波MIMO成像的快速图像重建问题,部分研究者提出了采用一定近似或是针对特定收发阵列的快速傅里叶变换类MIMO图像重建算法。2012年,Zhuge等^[11]推导了适用于近场MIMO阵列的距离迁移重建算法,基于快速傅里叶变换和频率插值实现了图像快速重建,但是该算法要求MIMO阵列的发射和接收单元满足奈奎斯特采样条件,难以应用到一般的稀疏MIMO阵列成像中。2014年,Alvarez等^[12]提出了适用于接收阵列满阵排布的MIMO雷达成像系统的快速傅里叶变换重建算法,该算法通过平面波近似进行局部区域聚焦,采用并行多次重建的方法对整个区域进行成像,该方法需要对重建精度和重建速度进行折中处理。文献[ 13-14]提出针对一维线阵MIMO扫描成像的快速傅里叶类MIMO重建算法。除了快速傅里叶变换类重建算法,还有另外一种更加通用且精度更高的重建算法——反向传播(BP)算法^[15]。BP算法对于收发单元的排布没有任何要求,因而适用于单站式成像或任意排布的MIMO阵列成像。对于重建图像中的每个像素,该方法将所有全息数据补偿对应空间传播相位因子后叠加求和,因此可以避免傅里叶变换类算法中离散傅里叶变换以及插值所引起的边缘干扰和重建误差,且能够获得更高的图像信噪比。BP算法成像效果好、适用面广,但是计算量过大,限制了其在各类实际成像系统中的应用。

针对毫米波全息成像的图像重建问题,本文提出了一种基于降维策略的快速反向传播重建(DR-BP)算法,该方法能够在保证图像重建精度且不改变BP重建通用性的前提下大幅度减少计算量。本文首先对毫米波全息成像原理进行了介绍,然后给出了DR-BP算法的原理及具体实现方法,最后通过仿真和实验进行了算法验证。

2 毫米波全息成像原理

与焦平面成像方法不同,毫米波全息成像通过发散波束照射成像目标,成像目标反射的毫米波信号被接收天线接收。基于外差混频技术,接收到的信号包含毫米波空间传播以及成像目标反射的复信号。按照收发模式的不同,毫米波全息成像一般可以分为单站和多发多收两种。在单站式成像中,收发天线近似位于同一位置,可以通过线阵阵列和机械扫描的方式实现整个平面全息数据的采集。典型的平面扫描成像如图1(a)所示,假设成像目标的复反射率为f(x,y,z),收发天线的位置为(x₀,y₀,0),则接收到的全息数据可以表示为

U(x0,y0,ν)=∭f(x,y,z)exp(-j2kr)dxdydz,(1)

式中:ν为毫米波频率;k=2πν/c为毫米波波数;r为收发天线与成像目标上任一点之间的距离,表达式为

r=(x-x0)2+(y-y0)2+z2。(2)

对于MIMO阵列成像,发射天线和接收天线可以任意排布,示意图如图1(b)所示。假设发射天线和接收天线的位置分别为(x_T,y_T,0)和(x_R,y_R,0),则接收到的全息数据可以表示为

U(xT,yT,xR,yR,ν)=∭f(x,y,z)exp(-jkrT,R)dxdydz,(3)

其中T和R表示发射天线和接收天线的标号,r_T,R为发射天线和接收天线到成像目标上任一点的距离之和,表达式为

rT,R=(x-xT)2+(y-yT)2+z2+(x-xR)2+(y-yR)2+z2。(4)

毫米波图像重建算法就是基于(1)式和(3)式的成像模型由全息数据U反演获得成像目标的复反射率f(x,y,z)。

图 1. 毫米波全息成像示意图。(a)单站式成像;(b) MIMO成像

Fig. 1. Schematic diagram of millimeter-wave holographic imaging. (a) SISO imaging; (b) MIMO imaging

下载图片 查看所有图片

对于单站式成像,全息数据均匀排布在整个平面上,这一过程可以等效为成像目标虚拟场以半波长空间传播一定距离后得到的场分布^[16],进而可以通过角谱传播公式进行快速傅里叶变换重建,重建公式为

f(x,y,z)=FT3D-1FT2DUx0,y0,ν)]exp(-jkzZ1)},(5)

式中:FT_2D为针对(x₀,y₀)两个维度的数据进行二维快速傅里叶变换; k_z=4k2-kx02-ky02, kx0和 ky0为空间波数;Z₁为距离参数,用于选择聚焦平面;exp(-jk_zZ₁)表示进行相位补偿,补偿后可得三维频域数据;F T3D-1为三维快速傅里叶逆变换,即将频域数据变换为空间域三维图像。采用离散傅里叶变换进行聚焦成像时,若成像目标大于扫描孔径,则成像区域外的目标散射信号会对图像边缘造成严重的干扰。另一方面,在快速傅里叶变换重建中k_z分布不均匀,需要在频率空间进行插值,因而也会带来一定的重建误差,故图像的重建目标位置存在一定偏差^[17]。

对于MIMO阵列成像,发射天线和接收天线的排布不规则,快速傅里叶变换类的重建算法难以直接应用。不论对于单站式成像还是MIMO阵列成像,BP算法都能通过对全息数据补偿空间传播相位因子获得高质量图像。典型的反向传播重建公式为

f(x,y,z)=∑ν=1Nν∑R=1NR∑T=1NTU(xT,yT,xR,yR,ν)exp(jkrT,R)。(6)

由于(6)式是对所有全息数据相位补偿后求和,故反演得到的信号具有较高的信噪比。相比傅里叶变换类图像重建算法,BP算法重建得到的图像边缘干扰少、目标位置准确、信噪比高,同时适用于任意收发天线排布,因此在毫米波全息成像中具有重要应用。但是BP算法的计算量大,其时间复杂度为O(N_TN_RN_νN_xN_yN_z),与发射天线数目(N_T)、发射天线对应的接收天线数目(N_R)、频率点数(N_ν)以及待重建的三维图像像素数目(N_xN_yN_z)成正比,重建速度一直是亟待解决的难题。

3 基于降维策略的快速反向传播重建算法

为了加快BP算法的重建速度,提出了一种DR-BP算法。首先将(6)式写为

f(x,y,z)=∑R=1NR∑T=1NTPT,R,(7)

其中P_T,R表示发射天线T和接收天线R对像素f(x,y,z)的贡献,即

PT,R=∑ν=1NνU(xT,yT,xR,yR,ν)exp(jkrT,R)。(8)

对于一组收发天线形成的全息数据,其对不同像素的贡献值与r_T,R有关,即与发射天线、接收天线到像素点的距离之和有关。由(4)式可知:当r_T,R为定值时,符合该条件的像素分布在一个椭球面上,该椭球的两个焦点为(x_T,y_T,0)和(x_R,y_R,0);改变r_T,R,则椭球面长轴随之改变。在三维图像重建中,对于同一椭球面上的像素P_T,R值相同,因此可以将待重建图像的像素按照椭球面进行划分,如图2(a)所示。在二维图像重建中,平面与椭球面相交形成椭圆,二维图像的像素可以按照椭圆进行划分,如图2(b)所示。上述分析是基于MIMO阵列成像,而单站式成像是MIMO阵列成像的一种特例,在三维成像中椭球面变为球面,在二维成像中椭圆变为圆环,故图像像素仍然可以按照贡献因子P_T,R进行划分。

图 2. 三维以及二维重建降维策略示意图。(a)三维;(b)二维

Fig. 2. Illustration of the 3D and 2D dimensionality reduction strategies. (a) 3D; (b) 2D

下载图片 查看所有图片

基于上述像素划分方法,可以在毫米波全息成像反向传播重建中引入降维策略,算法步骤见图3,图中r_min为r_min组成的向量,f为由f(x,y,z)组成的向量。1)在DR-BP算法中,对于每组发射天线和接收天线,若重建区域覆盖收发天线x和y坐标,则可以快速计算r的最小值和最大值,最小值等于重建区域最小z值的2倍,最大值为收发天线与重建区域8个顶点之间距离和的最大值。2)在最小值和最大值范围内,基于dr进行离散化,并计算对应r_d的P_d值。3)将r按照最临近插值的方法找到对应r_d的坐标,从而获得该组收发天线对于像素(x,y,z)的贡献P_d。在上述算法中,主要的计算量在于对P_d的计算,相比BP算法,计算复杂度由N_xN_yN_z降维为N_d。

图 3. DR-BP快速重建算法步骤

Fig. 3. DR-BP fast reconstruction steps

下载图片 查看所有图片

由于最临近插值方法是将r近似为r_d,因此采用最临近插值方法的DR-BP重建精度取决于dr的大小。dr值可以采取经验选取的方法,一般而言,可以基于补偿的相位因子最大误差对dr有一定的约束,约束公式为

δφ=kdr/2≤π4,(9)

其中δφ为最大相位误差。在实际图像重建中,选取的dr越小,图像重建精度越高,但是N_d会增大,从而使计算量增加。为了进一步放松对dr的限制,减小误差,可以采用线性插值代替最临近插值,线性插值算法见图4。

图 4. DR-BP的线性插值计算方法

Fig. 4. Linear interpolation calculation method of DR-BP

下载图片 查看所有图片

4 实验验证

首先基于单站式成像实验对DR-BP算法进行验证。实验采用280~320 GHz高分辨亚毫米波全息成像系统^[18],其中喇叭天线为对角喇叭,天线波束半张角为6.1°,系统成像分辨率可达2.5 mm×2.5 mm×4 mm,该系统的原理图和实物图如图5所示。采用该系统对金属仿真手枪模型[如图6(a)]进行成像,扫描步长为1 mm×1 mm,扫描范围为15 cm×20 cm,成像距离约为30 cm,扫描频点为200个。获取全息数据后,基于MATLAB对不同重建算法进行对比和验证,数据处理所采用的计算机的配置为CPU Intel I7-6700HQ@2.6 GHz,内存为16 G。

在金属仿真手枪模型的成像实验中获取的全息数据包含151×201×200个复数据,对于单站式成像,收发天线排布均匀,可以采用快速傅里叶变换重建方法进行快速重建,重建结果的最大投影如图6(b)所示。作为对比,采用BP算法对该数据进行重建,重建像素个数为151×201×15,结果如图6(c)所示。对比两图可以发现,在图像边缘部分(方框区域)由BP算法重建的图像更加平滑和清晰,而在重建速度方面,快速傅里叶变换算法远优于BP重建算法(重建时间对比见表1)。采用本文提出的DR-BP算法,基于最临近插值在dr分别取λ/2、λ/4和λ/8的情况下,重建结果如图6(d)~(f)所示,基于线性插值在dr分别取λ/2、λ/4和λ/8的情况下,重建结果如图6(g)~(i)所示,对应的重建时间见表1。采用应用广泛的结构相似性评价指标(SSIM)对重建结果的图像质量进行定量分析^[19-20](参考图像为BP重建结果),结果见表2。采用最临近插值方法,当dr分别取λ/2和λ/4时,图像均存在模糊现象,SSIM分别为0.4027和0.7033,而当dr=λ/8时,重建的图像与BP算法没有明显差别,SSIM为0.9104,而此时的重建时间仅为765 s,相比BP算法,计算速度提升83倍。dr=λ/4时,采用线性插值方法的重建结果与BP算法相当,SSIM为0.9814,相应地计算量有所增加,重建时间为1030 s,相比BP算法计算速度提升63倍。

图 5. 高分辨亚毫米波全息成像系统原理图以及实物图。(a)原理图;(b)实物图

Fig. 5. Schematic diagram and photo of the high resolution submillimeter wave holographic system. (a) Schematic diagram; (b) photo

下载图片 查看所有图片

图 6. 仿真金属手枪的光学图像及重建结果。(a)光学图像;(b) FFT重建结果;(c) BP重建结果;(d)~(f) DR-BP最临近插值dr分别取λ/2、λ/4和λ/8(从左到右)时的重建结果;(g)~(i) DR-BP线性插值dr分别取λ/2、λ/4和λ/8(从左到右)时的重建结果

Fig. 6. Optical image and reconstruction results of the simulated metal pistol. (a) Optical image; (b) FFT based reconstruction result; (c) BP reconstruction result; (d)-(f) reconstruction results of DR-BP nearest interpolation (the values of dr equal λ/2,λ/4 and λ/8 from left to right); (g)-(i) reconstruction results of DR-BP linear interpolation (the values of dr equal λ/2, λ/4 and λ/8 from left to right)

下载图片 查看所有图片

进一步地,采用MIMO成像的电磁仿真数据对本文提出的DR-BP算法进行验证。仿真实验基于FEKO电磁仿真软件,收发天线排布采用正交阵列排布方式,水平排布的为发射天线,竖直排布的为接收天线,某一个发射天线发射毫米波信号时,所有接收天线同时对发射信号接收,仿真成像示意图如图7(a)所示。仿真实验中,电磁波频率设为70~80 GHz,扫描频点数为81个,收发天线采用偶极子全向天线各250个,相邻天线之间间隔为2 mm,成像距离设为1 m,电磁仿真算法采用大面元物理光学法。对于MIMO收发数据,无法采用基于快速傅里叶变换的快速重建算法直接进行重建,而采用BP或本文提出的DR-BP方法能够获得很好的重建结果。采用BP算法对仿真手枪的重建结果如图7(b)所示,重建区域像素数目为151×214×5,从图中可以看到清晰的仿真手枪图像,但是由于MIMO成像中采用了稀疏的天线收发单元,枪头及枪把等一些地方的细节有些模糊。采用DR-BP算法对该数据重建,重建结果如图8所示,图8(a)~(c)和图8(d)~(f)分别为dr取λ/2、λ/4和λ/8时采用最临近插值和线性插值的重建结果,对应的重建时间见表1、SSIM值见表2。从结果中可以看到,当dr取λ/4时两种方法得到的重建结果均与BP算法相当,SSIM分别为0.9259和0.9957(线性插值结果更加准确),此时相对BP算法计算速度提升60和36倍。

图 7. 正交阵列全息成像仿真示意图以及BP算法重建图像。(a)仿真示意图;(b)重建图像

Fig. 7. Simulation schematic diagram of the orthogonal array holographic imaging and reconstructed image of BP algorithm.(a) Simulated image; (b) reconstructed image

下载图片 查看所有图片

图 8. 正交阵列FEKO电磁仿真成像实验DR-BP重建结果。(a)~(c)最临近插值方法dr分别取λ/2、λ/4和λ/8(从左到右)时的重建结果;(d)~(f)线性插值方法dr分别取λ/2、λ/4和λ/8(从左到右)时的重建结果

Fig. 8. DR-BP reconstruction results of orthogonal array imaging based on FEKO simulation (the values of dr equal λ/2,λ/4 and λ/8 from left to right). (a)-(c) Nearest interpolation; (d)-(f) linear interpolation

下载图片 查看所有图片

综上,基于实际实验和仿真实验数据对DR-BP算法进行了验证。对于单站式成像DR-BP算法能够获得优于FFT算法的图像质量,同时相比BP算法重建速度有了大幅度提升;对于MIMO式成像,DR-BP算法能够实现对任意收发天线排布的成像数据进行重建,在获得同样图像质量的情况下,大幅度提升了BP算法的重建速度。

5 结论

提出了一种基于降维策略的毫米波全息成像反向传播快速重建(DR-BP)算法,该算法能够在不降低反向传播重建算法(BP)图像重建质量的情况下,提升了重建速度。通过实际实验以及仿真实验对DR-BP算法进行了验证,DR-BP算法对毫米波全息成像的成像模式以及收发天线排布没有要求,对于单站式成像,能够获得优于FFT重建算法的图像重建结果,对于MIMO式成像,能够大幅提升图像重建速度,因而更具有实用性。为了进一步提升算法速度,DR-BP同样可以结合图形处理器(GPU)加速技术进行并行计算。需要注意的是,DR-BP算法在每步循环中需要对待重建图像进行累加,因而在并行计算中需要使用并行技术中的原子累加操作,如CUDA里的atomicAdd函数。