1.   引言

近年来，由于国际和国内公共安全形势的严峻性，火车站、机场、地铁等公共安全场所的安检显得尤为重要^[1]。金属探测器、X光安检机等传统安检手段因其各自的局限性，不适合或者不能高效地用于人体安检^[2,3]。毫米波具有穿透能力强、对人体辐射小且分辨率高等先天优势，尤其适用于隐藏物体检测^[4,5]，因此毫米波成像为人体安检提供了新的选择，并且具有广阔的应用前景。主动式毫米波雷达3维成像基于合成孔径雷达(Synthetic Aperture Radar, SAR)成像技术^[6,7]，实现对人体快速3维成像并进行异物检测，毫米波3维成像安检系统受到了国内外广泛重视。

毫米波3维成像系统有多种实现方式，国内外研究比较深入且具备实际工程应用能力的主要是平面扫描成像和圆柱扫描成像^[7]，平面扫描成像虽然具备3维分辨能力，但受限于波束扫描范围，因此不能通过一次扫描实现全方位的观测角度成像^[1]。圆周扫描成像系统采用阵列SAR与圆迹SAR相结合的方法，从理论上解决了毫米波成像系统的3维分辨问题，并弥补了平面扫描方式的不足，在较小的系统代价条件下实现高性能3维成像^[8]。目前国内研究重点多集中于圆周扫描成像系统，文献[6]给出了一种主动式毫米波圆周扫描3维成像的系统模型以及重构算法，并且搭建出了基于该模型的实际成像系统。文献[9]基于圆柱面扫描给出了一种基于巴克码稀疏采样和频域压缩感知的毫米波人体安检3维成像方法，并设计了一种基于巴克码和收发分置模式的稀疏平面阵列，得到了较好的3维成像实时处理效果。本文针对圆周扫描3维成像系统的特点，结合上述研究背景，分析了3维ω-K算法和3维后向投影(Back Projection, BP)算法在工程应用中的限制条件，提出了一种将ω-K算法和BP算法相融合的3维图像重构算法。

2.   现有成像算法分析

如图1所示，在主动式圆柱扫描毫米波3维成像系统中，两个天线阵列以Z轴为轴心，以R为半径做圆周旋转，形成一个圆柱扫描面，其中一列线阵作为发射天线，另一列线阵作为接收天线，天线阵列位置与X轴的夹角定义为 $\theta$ , $\theta \in \left[ {0,{{π}} } \right]$ 。该系统在Y方向和X方向的分辨率依靠雷达的距离向分辨率和沿角度 $\theta$ 方向的合成孔径获得。在天线阵列旋转的同时通过收发开关来控制天线的收发时序，从而在竖直方向形成合成孔径获得Z方向的高分辨，Z方向合成孔径的长度与波束的俯仰角宽正相关，假设波束俯仰角宽为S，则竖直方向的合成孔径长度为 ${{\rm{L}}_z} = 2R\tan \left( {{\rm{S}}/2} \right)$ 。

图  1  圆柱扫描3维成像系统模型

Figure  1.  Three-dimensional imaging system model of cylindrical scanning

下载: 全尺寸图片 幻灯片

对上述模式成像系统重构算法的选取，主要应该考虑算法对雷达运动模式的适用性和成像算法的聚焦性能。BP算法在成像过程中没有近似，成像精度较高，同时对雷达的运动方式没有特殊要求，并且不受转动角度的限制，因此非常适用于本系统的工作模式^[10]。3维BP成像处理对于网格点数较多的3维场景成像运算量巨大，针对高实时要求的人体安检应用存在困难。文献[6,11]提供了一类基于3维频域(波数域)成像算法(3维ω-K算法)，可以快速实现3维成像，该类算法在计算3维逆傅里叶变换时要求在角度维均匀采样，非线性相位补偿比较复杂^[12]。在柱面扫描条件下，要求天线阵列扫描轨迹为理想圆柱，一定程度限制了系统工作方式，对天线阵列运动轨迹要求严格。本文结合两种算法的优点，提出了ω-K算法和BP算法相结合的方式来实现快速柱面扫描方式的3维成像。在阵列(高度)方向与距离方向采用ω-K算法进行解耦合与聚焦，减少后向投影的运算量；对于每一个阵元，在距离方向和角度方向采用BP算法进行合成孔径处理，避免角度维插值处理，同时适用于非理想圆柱扫描方式。

3.   融合ω-K和BP 3维重构算法

设某一脉冲时刻，雷达天线位置为 $\left( {R\cos \theta ,} \right.$ $\left. { R\sin \theta ,z'} \right)$ ，成像区域第 $i$ 个目标点的位置为 ${P_i}\left( {{x_i},{y_i},{z_i}} \right)$ ，散射系数为 ${\sigma _i}$ ，则成像区域第 $i$ 个目标点到雷达天线中心的斜距为：

该时刻雷达天线发射的线性调频信号为：

式中， $\widehat t$ 为快时间， ${T_{\rm{p}}}$ 为脉冲宽度， ${{\rm{f}}_{\rm{c}}}$ 为信号载频， $\gamma$ 为线性调频信号的调频率，则在 $\left( {t,\theta ,{z'}} \right)$ 域中表示的 ${P_i}$ 点的回波信号为：

式中， $\Delta t = 2{R_i}/{\rm{c}}$ 为雷达天线发射信号与成像区域目标点的双程时延。对式(3)进行解线频调(dechirp)处理，选取的参考信号为：

式中， ${R_{\rm{ref}}}$ 为参考距离。参考信号与回波共轭相乘，即作差频处理，回波变成单频信号，且其频率与回波和参考信号的距离差成正比。若令 ${R_{\Delta {{i}}}} = {R_i} -$ ${R_{\rm{ref}}}$ ，则式(3)与式(4)共轭相乘后的差频信号为：

实际上回波数据是由成像区域内的多个散射点的回波相叠加而成的，于是式(5)改写为：

式(6)即为实际采集的回波信号，式中相位项第1项为差频之后的单频项，对应频率为 ${f_i} =$ $- \gamma \left( {2{R_{\Delta {{i}}}}/{\rm{c}}} \right)$ , ${R_{\Delta {{i}}}}$ 的变化会使对应的距离项中的频率即 ${f_i}$ 发生改变，对其沿着快时间做傅里叶变换，即可将不同距离处的目标分离出来，即脉冲压缩得到1维距离像，其数学表达为：

式(6)第2项的相位变化使回波产生多普勒，对应同一个目标不同脉冲时刻的相位变化。第3项为解线频调带来的剩余视频相位(Residual Video Phase, RVP)，可以通过相位补偿的方式将其补偿掉。而在本系统中由于在量级上 ${\rm c} \gg {R_{\Delta {{i}}}}$ ，所以此项也可以忽略不计。式(6)为回波数据的时域表达形式，如果用波数域的表达，同时补偿掉RVP项，其数学表达式为：

式中， $\Delta {K_r} = 4{{π}} \gamma \widehat t/{\rm{c}} \in \left[ { - 2{{π}} \gamma {T_{\rm{p}}}/{\rm c},2{{π}} \gamma {T_{\rm{p}}}/{\rm{c}}} \right]$ , ${{{K}}_{{\rm {rc}}}} =$ $4{{π}} {{\rm{f}}_{\rm{c}}}/{\rm{c}}$ , ${K_r} = {{{K}}_{{\rm{rc}}}} + \Delta {K_r}$ 为频率在空间的表达，即波数表达方式，为了得到回波信号关于距离-高度两个维度的2维频谱，对式(8)沿着Z轴做傅里叶变换到竖直方向的波数域可得：

为了下述表示方便，用 ${R_{{\rm{xyz}}}} = {R_{\Delta {{i}}}} + {R_{{\rm{ref}}}}$ 表示在3维坐标系下的目标点到雷达天线的斜距，用 ${R_{{\rm{xy}}}} = \sqrt {{{\left( {R\cos \theta - x} \right)}^2} + {{\left( {R\sin \theta - y} \right)}^2}}$ 表示竖直方向某一天线阵元到该阵元所在XOY平面上目标点的斜距，则， ${R_{{\rm{xyz}}}} = \sqrt {R_{\rm xy}^2 + {{\left( {z' - {z_i}} \right)}^2}}$ 对式(9)乘以 ${{\rm{e}}^{ - {\rm{j}}{K_r}{R_{\rm{ref}}}}}$ 来补偿由于去斜参考时间延迟带来的相位偏移。采用驻定相位原理导出式(9)的频谱解析式：

式(10)的相位如下：

式(11)包含两项，第2项线性相位项代表目标高度位置的响应；第1项是距离和Z方向耦合产生的RCM(Range Cell Migration)，代表目标的倾斜距离对RCM空间的变化。使用ω-K算法中的变量代换对式(11)中的第1项耦合相位解耦合：

式(12)叫做Stolt映射(Stolt Mapping, SM)，通常采用插值来实现。由于SM移除了距离方向与高度方向的耦合，同时使得在波数域 ${K_{z'}}$ 与 ${K_{{\rm{xy}}}}$ 垂直^[13]。经过SM，式(10)变为式(13)形式：

式(13)通过距离IFT(Inverse Fourier Transform)处理可以将不同距离点目标在不同的距离单元区分开来，得到距离脉冲压缩的结果，再经过高度维IFT将信号变回慢时间域，实现高度维聚焦。对式(13)分别沿着距离和高度方向进行IFT，同时考虑到高度方向的窗函数，得到距离高度2维脉冲压缩结果具有式(14)形式：

式(14)中 ${A_i} = A{\sigma _i}$ , A表示信号处理过程中幅度系数。 $\tau = 2{R_{\rm {xy}}}/{\rm c}$ ，表征了目标散射点距离雷达相位中心的距离。通过上述ω-K算法，实现了距离与高度维的解耦合，此时可以将3维聚焦问题转化为在竖直方向的聚焦问题和在XOY平面内的投影问题。具体做法是：将3维成像区域沿竖直方向划分为N个XOY切片堆叠的形式，如图2所示。此时只需要对每个阵元对应的XOY切片都执行相同的聚焦操作，即可得到3维的聚焦图像，因此问题转化为对每个切片的聚焦问题。由于天线阵列的运动轨迹为圆周，因此选取对曲线运动适应性较好的BP算法来实现每个切片的2维SAR图像。

图  2  3维场景成像切片划分

Figure  2.  Imaging slice division of three-dimensional scene

下载: 全尺寸图片 幻灯片

BP算法的基本思想是先通过距离脉冲压缩，使不同距离的目标点分离开来，然后在方位向，通过对回波信号补偿相位进行相干积累，从而将处于同一距离的不同散射点区分开来，即实现方位向的聚焦。BP算法的具体实现方式如下，由于通过ω-K算法已经实现了竖直方向与距离方向的解耦合，因此对于每一个垂直于Z轴方向的2维XOY切片来说，式(14)可以简写为：

采用函数 $f\left( {x,y} \right)$ 来表示2维XOY切片上像素点位置 $\left( {x,y} \right)$ 的散射函数，由时域相关原理可知，位于 $\left( {{x_i},{y_j}} \right)$ 处点目标的散射函数可以由式(16)重建：

式(16)说明了位于 $\left( {{x_i},{y_j}} \right)$ 处点目标的散射函数的重建过程，对于解线频调信号由式(7)和式(14)可知，后向投影之前要补充相位因子才能重建散射点的sinc响应并完成相干积累。 ${R_{{{ij}}}}\left( \theta \right)$ 即是成像网格点 $\left( {{x_i},{y_j}} \right)$ 到雷达天线阵元的斜距， ${R_{{{ij}}}}\left( \theta \right) =$ $\sqrt {{{\left( {R\cos \theta - {x_i}} \right)}^2} + {{\left( {R\sin \theta - {y_j}} \right)}^2}}$ ，该值是天线阵元位置即 $\theta$ 的函数，对于每一个XOY切片，运用上述BP算法，即可重建该切片的2维SAR图像，然后将各切片沿Z轴方向堆叠，可以重建3维SAR图像。最后，为了便于对场景目标观察和检测，对重建的3维SAR图像按照不同的观测角度进行2维投影，即可得到全方位观测角度下的2维图像，完整的算法流程如图3所示。

图  3  3维场景重构算法流程

Figure  3.  Algorithm flow of three-dimensional scene reconstruction

下载: 全尺寸图片 幻灯片

4.   算法仿真与验证

4.1   多角度成像仿真

按照上述算法流程，使用Matlab对3维人体模型进行回波仿真并重构3维图像，以此验证算法的可行性。仿真时，天线阵列以固定的角速度沿圆周运动，设圆周半径为1 m，扫描角度范围为 $\theta \in \left[ {{{π}} /6,5{{π}} /6} \right]$ 。仿真所用3维人体模型及运动天线阵列3维轨迹如图4所示。

图  4  3维人体模型及阵列运动轨迹

Figure  4.  Three-dimensional human model and the trajectory of antenna array

下载: 全尺寸图片 幻灯片

成像时，对所有回波数据进行ω-K聚焦之后，BP投影采用了数据复用的方式实现多个视角下的成像，图5为使用仿真数据，采用文中所述重构算法得到的3维人体模型仿真结果。

图  5  多视角成像结果

Figure  5.  Multi-view imaging results

下载: 全尺寸图片 幻灯片

由图5可知，重构算法较好地完成了3维场景的重建，3维人体模型的轮廓清晰，直观上看，图像聚焦效果较好。并且不同角度下的3维人体模型重构图像质量没有损失，说明本文算法可以实现3维目标的多角度成像，完成了算法理论上可行性的验证。

4.2   非理想圆柱仿真

为了验证本文算法对天线阵列非理想圆柱扫描轨迹的适应性，仿真了一组阵列沿椭圆轨迹运动的3维人体回波数据，椭圆长轴设为1 m，短轴设为0.8 m。分别采用3维ω-K算法和本文算法进行对比，成像结果如图6所示。

图  6  成像结果对比

Figure  6.  Comparison of imaging results

下载: 全尺寸图片 幻灯片

由图6可知，在天线阵列圆周轨迹运动条件下，本文算法和3维ω-K算法均能较好地完成3维成像；但在天线阵列椭圆轨迹运动条件下，3维ω-K算法由于3维匹配滤波以及波数均匀化不准确等原因导致图像散焦和镜像，而本文算法依然能够完成精确聚焦，说明本文算法在非理想圆柱扫描模式下，依然能够较好地完成人体3维成像。

4.3   实测数据分析

为了验证本文算法在实际安检成像系统的工程可用性，对采集的人体回波数据进行了合成孔径成像处理，并分析了本文算法与3维频域算法成像结果的区别。所选实际安检成像系统工作在Ka波段，系统带宽为6 GHz, X, Y, Z 3个维度的理论分辨率分别为5 mm, 25 mm, 5 mm。圆周轨迹运动，扫描半径为0.8 m，脉冲重复频率为200 kHz，阵列圆迹运动采用先加速后减速方式，因此在等时采样条件下，角度维的采样间隔是非均匀的。成像结果对比如图7所示。

图  7  圆柱扫描实测数据成像结果对比

Figure  7.  Comparison of imaging results of measured data from cylindrical scanning

下载: 全尺寸图片 幻灯片

从两种算法的对比结果可知，两种算法均能够实现图像聚焦，人体口袋里隐藏的刀具和枪支模型等都可以清晰呈现。相比于3维频域算法，由于BP算法不需要在角度维进行插值操作，也不存在斜距上的近似，因此保留了更多的细节，成像之后场景更为丰富。同一观测角度下，图像中人体腿部的散焦问题也得到明显改善。

从运算量上分析，不同于3维ω-K算法的运算量主要集中于FFT，本文算法的运算量主要集中于BP投影。若假设X, Y, Z 3维成像区域像素点数分别为Nx, Ny, Nz，距离、竖直、角度3个维度的采样点数分别为Nr, Ne, Na，则3维BP算法的运算复杂度为 ${\rm{Nx}} \times {\rm{Ny}} \times {\rm{Nz}} \times {\rm{Ne}} \times {\rm{Na}}$ 量级。本文算法BP投影部分运算复杂度为 ${\rm{Nx}} \times {\rm{Ny}} \times {\rm{Nz}} \times {\rm{Na}}$ 量级。由于算法具备很好的并行结构设计，实测数据处理表明(3维成像网格大小为 $380 \times 200 \times 64$ )，在NVIDIA GTX1070平台上采用CUDA(Compute Unified Device Architecture)编程成像计算时间约为0.18 s，满足实时成像的需求。

5.   结语

主动式毫米波3维成像系统在安检中具有对人体安全，穿透性好，成像效果佳等优势，因此具备广泛的应用前景，也是目前3维成像系统的发展方向。本文结合主动式圆柱扫描毫米波3维成像系统的工作方式，提出了一种ω-K算法和BP算法相结合的3维场景重建算法，给出了其理论依据。对3维人体模型进行仿真验证，结果表明，本文算法能够完成3维场景的聚焦成像，具备理论上的可行性。同时对比了本文算法与3维ω-K算法在非理想圆柱条件下的成像效果，验证了本文算法对系统运动方式的普适性。然后进一步在实际的3维成像系统采集的回波数据进行了验证，确保了算法在工程应用中的可行性，为该体制下的3维成像系统提供了一种新的简单易行的场景重构算法。此外，将本文提到的算法进行了CUDA工程实现，实现了3维人体成像的实时处理。