1.   引言

人体安检成像技术是近年来安检领域的研究热点^[1]，其中主动式人体安检成像系统因其受环境因素影响小、图像信噪比高等优点而被广泛应用^[2]。由于违禁物品往往具有隐蔽性，而毫米波能够穿透衣物和部分绝缘体遮挡、回波衰减率较小、对人体无害，因此利用毫米波作为载体的人体安检成像技术已成为目前最具潜力的安检成像技术之一^[3]。

2003年，美国L₃通信公司研制了ProVision毫米波人体扫描仪^[4]。该设备工作频率为 24～30 GHz，所采用的图像处理技术可获得5 mm横向分辨率，目前在美国、欧洲及亚洲等部分机场均有使用。2011年，德国埃尔兰根-纽伦堡大学微波光子研究所搭建了多基地毫米波快速成像系统^[5,6]。该原型系统工作频率为70～80 GHz，所获得的成像结果横向分辨率为2 mm；采用一个2 m×1 m大小的2维平面稀疏阵列，平板式数据获取时间可达到50 ms以下，可实现人体目标实时成像。

本文在上述研究背景下，首先基于柱面扫描成像模型，提出一种基于巴克码的稀疏采样方式和基于干涉处理^[7]频域压缩感知(Compressed Sensing, CS)的频域稀疏毫米波人体安检3维成像方法，以降低安检系统数据采集量，实际数据处理结果验证了算法的有效性。其次基于文献[5]的快速成像系统，设计一种基于巴克码和收发分置模式的稀疏平面阵列，进一步降低硬件成本。

2.   成像模型及信号处理

2.1   柱面扫描3维成像模型

毫米波人体安检柱面扫描3维成像模型如图1所示，图中以人体中心为坐标原点建立直角坐标系，其中X轴表示方位向(横向)，Y轴表示距离向，Z轴表示高度向(纵向)；1维天线阵列长度为L，且以坐标原点为圆心、测试距离R_a为半径进行扫描，形成圆心角为 $\varPhi$ 的柱面孔径。柱面扫描3维成像模型利用阵列天线结构获得人体高度向高分辨率图像，利用阵列天线沿圆弧上的机械扫描获得方位向高分辨率图像，利用各子阵发射并接收毫米波宽带信号实现距离向高分辨率成像。柱面孔径中任意采样点 $P_a$ 的坐标为 $({x_a},{y_a},{z_a}) \,=\,({R_a}\cos {\theta _a},{R_a}\sin {\theta _a},\,{z_a})$ $\left( {{\theta _a} \in \left[ {{\rm{ - }}\varPhi /2,\varPhi /2} \right],\,{z_a} \in \left[ {0,L} \right]} \right)$ ；人体成像区域任意点P_n的坐标为 $({x_n},{y_n},{x_n})$ ，散射强度为 ${ \sigma_n}$ 。与平面扫描3维成像相比，柱面扫描3维成像技术有多角度观察的能力。

图  1  毫米波人体安检柱面扫描3维成像模型

Figure  1.  The cylindrical scanning 3-D imaging model of a millimeter wave human security system

下载: 全尺寸图片 幻灯片

为减少数据获取时间、提高扫描速度，在方位向采用稀疏采样方式，但稀疏采样会使图像产生栅瓣和高副瓣。选择合适的采样准则可改善图像质量，本文采用巴克码伪随机序列^[8,9]作为稀疏采样准则。巴克码具有尖锐的自相关函数和良好的随机性，稀疏率约50%，其中长度为7的巴克码[1110010]对应的自相关函数图和方向图如图2所示。从图中可知以巴克码为稀疏采样准则可有效避免图像栅瓣和副瓣带来的能量泄露。

图  2  巴克码[1110010]的自相关函数图和方向图

Figure  2.  The Autocorrelation and Beam pattern diagrams for Barker code of [1110010]

下载: 全尺寸图片 幻灯片

2.2   回波信号处理

假设天线发射信号为 $p(t)$ ，则t时刻在采样点P_a处的成像区域回波信号为：

其中，D为人体成像区域范围。对式(1)在快时域进行傅里叶变换

其中， $P(\omega )$ 为发射信号的傅里叶变换， ${k_\omega }{\rm{ = }}\omega /{\rm c}$ , $\omega$ 为距离向波数，c为光速。定义k_x, k_y和k_z分别为X, Y和Z方向的波数分量， $k_x^2 + k_y^2 + k_z^2 = k_\omega^2$ ; k_r为X-Y平面的波数分量， $k_r^2{\rm{ = }}k_x^2 + k_y^2$ ，则式(2)中的指数项可写成平面波信号累加形式^[10]

假设 $P(\omega )$ 近似为狄拉克 $\delta$ 函数，则式(2)可写成关于距离向波数 $\omega$ 、高度向波数 ${k_z}$ 、方位向 ${\theta _a}$ 的函数

式中， ${F_\sigma }({k_r},\varphi ,{k_z}) \equiv {F_\sigma }(2{k_r}\cos \varphi ,2{k_r}\sin \varphi ,{k_z}) \equiv$ ${F_\sigma} ({k_x},{k_y},{k_z})$ 为空间域散射强度 ${\sigma _n}$ 的波数域谱。定义

则有

对 $g({\theta _a},{k_r})$ 与 ${F_\sigma }({k_r},{\theta _a},{k_z})$ 关于 ${\theta _a}$ 做1维傅里叶变换，则式(6)可写为：

由式(7)可知，构造形如 $1/G({\phi _a},{k_r})$ 的波数域匹配滤波器，对回波数据频谱 $S(\omega ,{\phi _a},{k_z})$ 进行匹配滤波，假设系统沿方位向均匀采样，对 ${\phi _a}$ 进行1维逆傅里叶变换，可得到散射强度 ${\sigma _n}$ 的3维波数域值 ${F_\sigma }({k_r},{\theta _a},{k_z})$ ，此方法称为波数域匹配滤波(Match Filtering, MF)成像算法。对 ${F_\sigma }({k_r},{\theta _a},{k_z})$ 进行Stolt插值，并进行3维逆傅里叶变换，即可得到空间域散射强度 ${\sigma _n}$ ：

(8)

3.   频域稀疏3维成像算法

当方位向采用稀疏采样方式时，采用上述基于波数域匹配滤波(Match Filtering, MF)的成像方法会引起图像混叠。对此本文提出一种基于干涉处理和频域CS的3维成像算法，在稀疏采样条件下重建图像。由式(6)可知，稀疏采样得到的回波信号可写成如下离散形式：

其中，M为方位向成像单元个数。在第i个距离向波数 ${\omega _i}$ ( $1 \le i \le {N_r}$ , ${N_r}$ 为距离向采样点数)、第j个高度向波数 ${k_{{z_j}}}$ ( $1 \le j \le {N_z}$ , ${N_z}$ 为高度向采样点数)处的回波信号为：

(11)

其中， ${k_r} = \sqrt {k_{{\omega _i}}^2 - k_{{z_j}}^2}$ 为距离向波数 ${\omega _i}$ 和高度向波数 ${k_{{z_j}}}$ 对应的X-Y平面波数分量。式(10)–式(12)构成了回波信号在频域方位向的CS线性测量模型，其中 ${{φ} _k}$ 为观测向量( $1 \le k \le {N_a}$ , ${N_a}$ 为方位向采样点数)， ${{α}}$ 为方位向目标散射强度的频域系数向量。CS线性测量模型写成矩阵形式为：

其中， ${{{S}}_{ij}}{\rm{ = }}{\left[ {S({\theta _{{a_1}}})\;S({\theta _{{a_2}}})\;\; ·\!·\!· \;S({\theta _{{a_k}}})\;\; ·\!·\!· \;S({\theta _{{a_K}}})} \right]^{\rm{T}}} \in$ ${\mathbb{C}^{{N_a} \times 1}}$ 为距离向波数 ${\omega _i}$ 、高度向波数 ${k_{{z_j}}}$ 处的方位向回波向量，H为基于巴克码采样准则的稀疏选择矩阵， ${{Φ}} {\rm{ = }}{\left[ {{\phi _1},{\phi _2},\; ·\!·\!· ,\;{\phi _k},\; ·\!·\!·,\; {\phi _K}} \right]^{\rm{T}}}\!\! \in \!\!{\mathbb{C}^{{N_a} \times M}}$ 为观测矩阵。

由于人体目标在空间域为连续分布，其图像在频域应具有稀疏性。但在毫米波人体复图像中，由于分辨单元间复散射系数不同、空间采样间隔通常远大于波长，其图像分辨单元初始相位是随机变化的^[11]，该随机初始相位主要是由于斜距方向分辨单元间复散射系数不同和斜距方向采样间隔较大而产生，由此使复图像频谱 ${{α}}$ 较宽难以压缩，因此无法利用CS理论求解。本文考虑利用稀疏采样中部分满采样回波数据构造参考复图像，并对全部采样数据对应的复图像进行干涉处理^[12,13]。假设待重建人体复图像3维频谱 ${{α}}$ 对应的空间域3维图像为 ${{σ}}$ ，采用基于波数域MF的成像方法对稀疏采样中部分满采样回波数据进行成像，得到参考复图像 ${{{σ}} _{\rm ref}}$ 。对两幅图像 ${{σ}}$ 和 ${{{σ}} _{\rm ref}}$ 进行干涉处理，得到干涉后的图像 ${{{σ}} _{\rm new}}$

由于两幅图像 ${{σ}}$ 和 ${{{σ}} _{\rm ref}}$ 对应的回波数据中天线视角相同，随机初始相位近似相等，故 ${{{σ}} _{\rm new}}$ 相位中已去除随机初始相位的大部分影响。假设图像 ${{{σ}} _{\rm new}}$ 对应的频谱为 ${{{α}} _{\rm new}}$ ，则 ${{{α}} _{\rm new}}$ 具备稀疏性，可通过频域CS进行恢复^[14,15]。

由式(14)可知，待重建复图像 ${{σ}}$ 可以表示为：

令 ${{σ}} {\rm{ = }}{{{F}}^{{\rm{ - 1}}}}{{α}}$ , F^–1为逆傅里叶基，则

即

综合式(13)–式(17)，可得

其中， ${{{α}} _{\rm new}}(i,j)$ 为3维频谱中在第i个距离向波数 ${\omega _i}$ 和第j个高度向波数 ${k_{{z_j}}}$ 处的方位向1维频谱， ${{P}}{\rm{ = diag\{ exp(}}{\rm j}\angle {{{σ}} _{\rm ref}}{\rm{)\} }}$ 为参考复图像相位构成的对角矩阵。经计算， ${{Φ}} {{FP}}$ 和F^–1的互相关系数低于0.0124，考虑稀疏采样矩阵H后二者的互相关系数可达0.0075，这表明上述模型中测量矩阵 ${{Φ}} {{FP}}$ 和F^–1是不相关的^[16–19]。由于3维频谱 ${{{α}} _{\rm new}}$ 具有稀疏性，基于CS理论，可通过求解如下所示的最优化问题，可求得干涉后的频谱，其中 $\varepsilon$ 为噪声上限。

在所有距离向波数域 ${\omega _i}$ 和高度向波数 ${k_{{z_j}}}$ 进行式(18)–式(19)求解，可得到人体安检复图像的3维频谱 ${{{α}} _{\rm new}}$ 。对频谱进行Stolt插值并逆傅里叶变换至空间域，可得到聚焦良好的人体安检3维复图像。

4.   柱面扫描实际数据处理

4.1   系统参数

本文所用的柱面扫描人体安检成像实际数据由华讯方舟科技有限公司提供，由该公司研制的柱面扫描型毫米波3维成像人体安检系统，如图3所示。系统中阵列天线长度约2 m，在方位向上的扫描角约66°，天线辐射单元在方位向和高度向波束宽度约60°，其余参数如表1所示。测试目标为一高约2 m、宽约0.60 m的人体模特，模特手持一块高度约0.30 m、宽度约0.50 m的金属分辨率板，胸前粘有金属细条。

表  1  系统参数

Table  1.  System parameters

参数	数值
测试距离Ra (m)	0.675
发射信号带宽Bs (GHz)	9.50
雷达中心频率f0 (GHz)	33.87
雷达工作波长 $\lambda$ (mm)	8.90
距离向采样率fs (GHz)	9.50
距离向采样间隔dr (cm)	1.58
方位向扫描角 ${{Φ}}$ (°)	66
方位向采样间隔da (°)	0.32
阵列长度L (m)	2
高度向采样间隔dz (mm)	4.80
距离向分辨率 (cm)	1.58
方位向分辨率 (mm)	4.50
高度向分辨率 (mm)	4.50

下载: 导出CSV 

| 显示表格

图  3  柱面扫描型成像系统和测试目标

Figure  3.  The cylindrical scanning imaging system and testing objectives

下载: 全尺寸图片 幻灯片

阵列天线以间隔d_a沿方位向均匀等间隔扫描，得到满采样数据；以扫描方向作为稀疏方向，采用长度分别为7, 11和13的巴克码伪随机序列作为稀疏采样准则，对满采样数据进行稀疏化(实际成像时可直接进行方位向稀疏扫描)。3种巴克码序列及性能如表2所示，基于7位巴克码的方位向稀疏采样示意图如图4所示，其中实心圆代表此处存在数据。

图  4  柱面回波数据稀疏方式示意图

Figure  4.  The sparse sampling method of cylindrical echo

下载: 全尺寸图片 幻灯片

表  2  3种巴克码序列

Table  2.  The 3 kinds of Barker sequences

巴克码长度	码型	稀疏率(%)
7	1 1 1 0 0 1 0	42.86
11	1 1 1 0 0 0 1 0 0 1 0	54.55
13	1 1 1 1 1 0 0 1 1 0 1 0 1	30.77

下载: 导出CSV 

| 显示表格

4.2   成像结果

满采样实际数据的直接成像结果如图5所示，其中图5(a)为人体3维成像结果的距离向投影，图5(b)为其中金属分辨率板的距离向投影，图5(c)为图5(b)红色框中5个金属条在距离向–0.12 m、高度向0.74 m处的1维切面。从图中可清晰看到人体和金属条形状，红色方框内的5个金属条方位向尺寸约为5 mm，间隔2 cm；红框左侧的6个金属条方位向尺寸约为3 mm，间隔1 cm。

图  5  满采样条件下的人体安检直接成像结果

Figure  5.  The direct imaging result under complete sampling

下载: 全尺寸图片 幻灯片

在扫描方向实施基于巴克码采样准则的稀疏采样，为了比较所提稀疏采样方法的有效性，本文同时采用等间隔稀疏采样方式。等间隔稀疏采样方式和基于7位巴克码准则的稀疏采样对应的直接成像结果分别如图6、图7所示。从图中可知，在等间隔稀疏采样方式下，图像方位向受栅瓣影响严重混叠，不能复原观测场景；在基于巴克码准则的稀疏采样下，图像无栅瓣影响，虽然方位向部分能量泄露至副瓣，但通过进一步的图像重建方法可恢复原始图像。

图  6  等间隔稀疏采样条件下的人体安检直接成像结果

Figure  6.  The direct imaging result under uniformly spaced sparse sampling

下载: 全尺寸图片 幻灯片

图  7  基于巴克码稀疏采样条件下的人体安检直接成像结果

Figure  7.  The direct imaging result under Barker code sparse sampling

下载: 全尺寸图片 幻灯片

在基于巴克码准则的稀疏采样下，采用本文所提频域稀疏3维成像算法对稀疏采样信号的重建结果如图7–图8所示，其中参考图像由7位巴克码中前3位对应的部分满采样数据获得。从图中可知，经过原始图像信号和参考图像信号的干涉处理，图像信号3维频谱已压缩变窄(图8)；对回波信号按照式(18)–式(19)进行频域CS稀疏重建，并将重建后的频谱逆变换至空间域，得到的成像结果没有高副瓣的现象，且分辨率接近满采样条件下的3维图像(图9)。

图  8  干涉处理前后图像3维频谱

Figure  8.  The 3-D frequency spectrum before and after interferometry

下载: 全尺寸图片 幻灯片

图  9  基于干涉处理和频域CS 3维成像算法的图像重建结果

Figure  9.  The reconstruction of 3-D imaging algorithm based on interferometry and CS in frequency domain

下载: 全尺寸图片 幻灯片

4.3   误差分析

上述成像结果表明，基于干涉处理和频域CS的频域稀疏3维成像算法能够在稀疏采样条件下恢复目标场景，且图像具有与满采样相当的分辨率水平。为了定量分析所提方法在稀疏采样下的图像重建性能，本文将满采样对应的图像近似作为目标真值，采用均方根误差(Root Mean Square Error, RMSE)、相关系数和结构相似度^[20](Structural SIMilarity, SSIM)3个指标对重建图像进行量化评价。均方根误差用来测量观测值和真值的偏差，定义为：

式中，A, A $'$ 分别代表满采样和稀疏采样下3维复成像的幅度矩阵。相关系数这一参数在雷达干涉领域经常被使用，定义为：

其中， ${\rm{cov(}}{{A}},{{A}}'{\rm{) \!=\! E}}\!\left[ {({{A}} \!-\! {\rm{E}}\!\left[ {{A}} \right])({{A}}' \!-\! {\rm{E}}\left[ {{{A}}'} \right])} \right]$ , ${\rm{E}}\left[ {\Large\cdot} \right]$ 为图像期望， ${\rm{D}}\left[ {\Large\cdot} \right]$ 为图像方差。两幅图像的相关系数越接近1，其线性相关性越强。结构相似度综合两幅图像的亮度(均值)、对比度(标准差)以及图像之间的相似性(协方差)来评价重建图像质量，其定义为：

(22)

其中，C₁, C₂为用来维持稳定的常数。结构相似度取值范围为[–1, 1]，两幅相同图像的结构相似度为1；两幅完全不相关图像的结构相似度为0；两幅完全相反的图像的结构相似度为–1。表3所示为基于7位巴克码准则稀疏采样条件下不同成像算法的图像重建性能评价结果。本文同时采用了长度为11和13的巴克码作为随机稀疏采样准则，对实际数据进行方位向稀疏化，并采用所提方法对稀疏回波数据进行图像重建；限于文章篇幅，仅列出不同长度巴克码稀疏采样下的图像质量评价结果(表4、表5)。从表中可知，在巴克码稀疏采样条件下，本文所提方法的重建图像接近满采样对应的成像效果，稀疏前后的图像相关系数大于0.9。

表  3  7位巴克码稀疏采样方式对应的图像质量评价

Table  3.  The image quality assessment corresponding to Barker code with length of 7

成像方法	相关系数	SSIM	MSE
等间隔稀疏采样直接成像	0.8698	0.8169	0.0431
7 位Barker码稀疏采样直接成像	0.8991	0.8973	0.0388
7 位Barker码稀疏采样干涉处理频域CS成像	0.9118	0.9036	0.0359

下载: 导出CSV 

| 显示表格

表  4  11位巴克码稀疏采样方式对应的图像质量评价

Table  4.  The image quality assessment corresponding to Barker code with length of 11

成像方法	相关系数	SSIM	MSE
11位Barker码稀疏采样直接成像	0.8768	0.8552	0.0491
11位Barker码稀疏采样干涉处理频域CS成像	0.9028	0.8676	0.0418

下载: 导出CSV 

| 显示表格

表  5  13位巴克码稀疏采样方式对应的图像质量评价

Table  5.  The image quality assessment corresponding to Barker code with length of 13

成像方法	相关系数	SSIM	MSE
13位Barker码稀疏采样直接成像	0.9159	0.8892	0.0343
13位Barker码稀疏采样干涉处理频域CS成像	0.9321	0.9052	0.0324

下载: 导出CSV 

| 显示表格

5.   用于快速安检成像的平面稀疏阵列设计

上述实际数据对应的成像系统基于柱面扫描3维成像模型，为减少数据采集量，增加安检成像通过速率，可采用基于实孔径平面阵列的毫米波快速成像系统。传统实孔径平面阵列为满足均匀半波长( $\lambda /2$ )间隔分布，系统需要大量天线辐射单元(假设雷达工作频率 ${f_0} = 30\;{\rm{GHz}}$ ，则1 m×1 m大小的阵列平面需要约40000个辐射单元)。为减少硬件系统规模，可采用收发分置的平面稀疏阵列布局。根据接收等效相位中心原理，收发分置的两个阵元中间位置会产生等效相位中心，从而可考虑对稀疏阵列天线布局进行优化。

5.1   2维稀疏阵列布局

基于频域稀疏成像方法、巴克码稀疏采样方式和收发分置工作模式，设计用于快速安检成像的稀疏阵列布局，其中阵列最小单元为6个L型天线阵列组成的正方形结构^[21]，如图10所示。图10(a)为一个L型收发分置天线阵列布局示意图，接收、发射天线阵元个数均为N，阵元间隔为d，在多发多收条件下共产生N×N个间隔为d/2的等效相位中心；图10(b)为阵列最小单元结构示意图，包含6个收发分置的L型天线阵列，在系统多发多收工作模式下产生满采样分布的等效相位中心。

图  10  2维稀疏阵列中最小单元结构示意图

Figure  10.  The structure of minimum unit in 2-D sparse array

下载: 全尺寸图片 幻灯片

基于图10所示最小结构，假设雷达工作波长 $\lambda {\rm{ = }}1\;{\rm{cm}}$ ，频率范围为25～35 GHz；最小单元大小为0.16 m×0.16 m，其中每条边的辐射单元个数为16个(L型天线阵列的发射单元和接收单元个数N=8)，间隔为d=1 cm，等效相位中心间隔为5 mm。为满足人体成像要求，需7×13个最小结构形成大小为1.12 m×2.08 m的2维阵列平面。为进一步减少辐射单元数量和系统复杂度，可在方位向采用基于巴克码准则的稀疏采样布局。以7位巴克码[1110010]作为稀疏准则的平面稀疏阵列布局如图11所示，其中阴影部分表示此处无天线收发阵元。

图  11  毫米波快速成像2维稀疏阵列布局示意图

Figure  11.  The 2-D sparse array structure for millimeter wave fast imaging

下载: 全尺寸图片 幻灯片

5.2   指标分析

本节对上述用于毫米波快速成像的平面稀疏阵列的主要性能指标进行分析。

(1) 成像分辨率

假设场景中心到平面阵列的测试距离为R₀，阵列平面大小为 ${L_x} \times {L_z}$ 。在人体3维成像模型中，天线发射/接收微波信号，利用合成的宽带信号获得距离向分辨率，通过天线孔径的平面阵列结构获得目标的方位向和高度向分辨率^[22]。对于工作波长为 $\lambda$ ，总带宽为B的天线发射信号，3维图像分辨率表达式如下：

其中， ${\rho _x}$ , ${\rho _z}$ 和 ${\rho _y}$ 分别代表方位向、高度向和距离向分辨率。假设工作波长为1 cm，测试距离R₀=0.675 m，总带宽B=10 GHz，则采用图10所示的2维稀疏阵列平面对应的图像方位向分辨率约为3.92 mm、高度向分辨率约为2.98 mm、距离向分辨率约为1.50 cm。

(2) 稀疏度

设稀疏阵列的物理辐射单元个数为N，稀疏阵列天线形成的等效相位中心均匀分布时所对应的空间位置个数为M，则稀疏阵列的稀疏率为：

当采用如图10所示的稀疏阵列和上述成像参数时，稀疏阵列的稀疏率约为94.6%。

6.   结论

本文提出了一种基于巴克码准则的稀疏采样方式以降低柱面扫描安检系统的信号采集时间。在稀疏采样条件下，提出一种基于干涉处理和频域CS的频域稀疏3维成像算法。利用干涉处理使复图像频谱具备稀疏性，在频域建立柱面3维成像CS模型，并对稀疏后的回波进行图像重建。实际数据处理结果表明：本文所提方法在数据采集量减少约50%的条件下，仍可获得接近满采样的图像分辨率和成像效果，稀疏采样前后的图像相关系数优于0.9，在大人流量场所安全检测中具有重要应用价值和广阔市场前景。

为减少扫描采样时间，增加安检成像通过速率，减少辐射单元数量和系统复杂度，本文设计了基于频域稀疏成像方法和巴克码稀疏采样方式的毫米波快速成像稀疏阵列布局，在收发分置工作模式下实现了人体实时安检成像。在保证成像质量的前提下，该2维稀疏阵列的稀疏率高达94.6%，对今后的安检成像系统具有重要意义。

致谢 本文研究工作得到国家自然科学基金(项目批准号：61271422)项目资助；本文所用的柱面扫描人体安检成像实际数据由华讯方舟科技有限公司提供。由华讯方舟科技有限公司研制的柱面扫描型毫米波3维成像人体安检系统的研制工作得到深圳市科技创新委员会(Science, Technology and Innovation Commission of Shenzhen Municipality, No.JSGG20160229114231872)和深圳市发展和改革委员会(Development and Reform Commission of Shenzhen Municipality Government, No. [2016]1588)项目资助。在此表示感谢。