稀疏轨迹毫米波雷达三维高分辨成像算法

马宇欣 海 宇 李中余 黄 鹏 王朝栋 武俊杰 杨建宇

(电子科技大学信息与通信工程学院 成都 611731)

1 引言

随着微波器件工艺和雷达系统技术的快速发展，毫米波雷达的应用得到了进一步的扩展。由于其优异的穿透性能[1]、较小的器件体积以及优异的成像分辨能力[2,3]，毫米波雷达被广泛应用在安检[4]、无损检测[5]、医学诊断[6]、自动驾驶[7]等领域。

毫米波雷达在应用于小目标成像识别时，为了达到多个维度的毫米级成像分辨，仅依赖发射高频段大带宽信号是难以实现的，因此许多学者提出了利用合成孔径原理，通过在高度和方位向合成二维实孔径的方式实现高的方位-高度分辨。根据成像分辨的计算式[8]，一方面，毫米波雷达通过发射中心频率较高的大带宽信号，利用增大带宽提升距离分辨；另一方面，通过采样操作，在方位向和高度向获得大的合成孔径同样能提高分辨。

实际应用中，为了达到高分辨成像的目的，往往利用平台搭载雷达，通过扫描合成大的方位-高度二维孔径，实现高度向和方位向的二维高分辨。当毫米波雷达平台通过多轨迹扫描采集回波数据时，增大轨迹间隔减少轨迹运动时间能降低操作的复杂度[9]，以上场景下采集的回波在高度维无法满足奈奎斯特定理，具有稀疏特性。这就导致当采用传统算法，如距离多普勒(Range Doppler,RD)、后向投影(Back Projection,BP)成像时，高度向的数据缺失会导致主瓣能量泄露、旁瓣抬升、目标混叠甚至无法成像。

为了解决该问题，常用的技术有压缩感知[10-12](Compressive Sensing,CS)和矩阵填充[13,14](Matrix Completion,MC)，其中压缩感知的概念于2006年被提出，它利用信号数据的冗余性，通过设计测量矩阵找到信号的稀疏表示，实现信号的稀疏恢复。在文献[15,16]中，多次采用压缩感知技术实现稀疏数据下的毫米波三维成像，但CS的稀疏恢复性能非常依赖测量矩阵的设计，当设计矩阵的稀疏基与场景中目标的分布不匹配时，恢复质量大打折扣。文献[16]提出了结合深度学习的恢复成像方法，对传统压缩感知的迭代软阈值算法(Iterative Soft Thresholding Algorithm,ISTA)展开设计unrolling网络，通过网络的训练，学习符合场景的测量矩阵，能够在数据稀疏度较大的场景高精度成像，不过网络的训练依赖大数量的样本并且耗时较长。

不同的是，MC是一种利用数据的低秩约束恢复缺失数据的技术，目前已经应用在医疗成像[17]、稀疏信道估计[18]、阵列信号处理[19]等领域。MC具有效率高、恢复精度高等优点，因而也被用于解决稀疏成像的问题。文献[20]提出了一种基于矩阵填充的二维稀疏毫米波三维成像方法，但矩阵填充技术不能从矩阵子集恢复原矩阵，即待恢复矩阵的每一行和每一列都要存在观测信息，因而不适用于某一维度的大量缺失。文献[21]提出了一种基于矩阵填充的稀疏线阵合成孔径雷达(Synthetic Aperture Radar,SAR)成像方法，通过对信号切面矩阵做旋转变换使得矩阵符合矩阵填充的工作前提，但由于仅考虑回波的低秩先验，矩阵填充恢复效果不佳，不适用于高度向稀疏度较高的场景。

考虑到上述方法在应用于稀疏轨迹毫米波雷达成像中仍存在不足，本文提出了一种基于Hankel矩阵构造和矩阵填充的成像算法，实现了高稀疏轨迹下的毫米波雷达回波重构以及高分辨成像。本文的主要工作和创新点主要体现在以下3个方面：

(1) 构建了稀疏轨迹毫米波雷达的回波模型，并分析了回波的稀疏-低秩先验特性，为成像算法的提出奠定了基础；

(2) 提出了一种基于稀疏-低秩先验的矩阵填充方法，利用回波高度-距离切面的稀疏-低秩特性，建立了基于截断Scatten-p范数与l0范数的双重约束优化模型，解决了传统矩阵填充由于无法对行列缺失矩阵恢复从而无法直接应用于回波重构的问题；

(3) 提出了一种基于Hankel矩阵变换和矩阵填充的高精度成像方法，通过引入Hankel矩阵构造解决了高度向缺失度较高时难以实现高精度数据恢复的问题。

2 稀疏轨迹毫米波雷达回波建模与信号特性分析

2.1 稀疏轨迹毫米波雷达回波建模

根据引言的介绍，毫米波雷达在实际的应用中，存在通过手持或机械平台搭载毫米波雷达扫描采集回波数据的场景，为了降低操作复杂度，通过增大轨迹间隔缩短轨迹时间，其扫描路径往往稀疏不确定性，其几何构型见图1。经过孔径校正后[22]，回波在方位向上均匀采样而在高度向上非均匀稀疏采样。其中X代表距离向，Y代表方位向，Z代表高度向，图中的圆圈代表毫米波雷达稀疏轨迹采样点，距离采样点数为K，方位采样点数为M，高度采样点数为N′。

理想情况下，为了满足奈奎斯特采样定理，毫米波雷达在方位向和高度向的采样点间隔需要小于等于半波长，三维回波采样布局图见图2，稀疏轨迹孔径校正后的高度向采样点采样间隔为半波长的整数倍，采样点数为N′，符合奈奎斯特采样定理的满采情况下，高度向采样间隔为半波长，所需的高度向采样点数目为N，因此轨迹稀疏度为I=N′/N，两种采样布局的高度孔径长度均为LN，因此理论可达的最高分辨一样，而数据量相对减少1-N′/N。

图2 采样布局图Fig.2 Diagram of sampling layout

毫米波雷达由于工作频段高，波长短，可以将目标视作多个散射点的集合。对于目标上的一个散射点Pp(xp,yp,zp)，采样点Pr(x0,ym,zn)发射线性调频信号，基于散射点模型，可以得到目标的回波信号表达式：

其中，t表示距离向时间，c是光速，fc是中心频率，P为目标散射点总数，σp表示第p个散射体的反射系数，wr(·)表示距离向窗函数，λ为信号波长，Kr为信号调频率，R(m,n,p)表示雷达采样位置Pr(x0,ym,zn)到目标上散射点Pp(xp,yp,zp)的距离历史，可以做以下近似：

其中，B是线性调频信号的带宽，由于不同方位和高度采样过程是独立的，通过对不同高度、方位的采样点采集回波可以组合得到最终的三维回波形式S3D∈CK×M×N：

其中，Z={z1,z2,...,zN′}为扫描轨迹的高度向采样点位置集合。

2.2 高度-距离切面信号的特性分析

同时，本节分析了采样点间隔为半波长时，三维回波高度-距离切面信号的先验特性。根据式(3)，对于第y0个方位时刻，经过距离压缩、距离徙动校正后回波的高度-距离切面写作：

其中，exp{-j2πfc(2R0/c)}常数可忽略，当在某一方位时刻时，exp[-j2π(/(λR0))]看作常数也可忽略，为了对式(5)中的高度-距离切面矩阵做进一步分解，对高度向的相位解耦合，即做相位补偿，补偿如下：

补偿后的矩阵可以分解为

其中，T ∈CK×P,Γ ∈CP×P,Z ∈CP×N，具体表达式如下：

2.3 Hankel矩阵构造及特性分析

为了充分利用切面矩阵的低秩先验特性以及高度向间相关性，引入Hankel构造变换，记作H(·)，变换示意图见图3。

图3 Hankel构造变换示意图Fig.3 Diagram of Hankel structural transformation

取高度-距离切面矩阵中的距离单元向量做Hankel构造变换，高度-距离切面的距离单元向量记作SH(t0,y0,z)=[SH(z1),SH(z2),...,SH(zN)]∈C1×N，经变换后为

其中，C=sinc[B(t0-2R0/c)]可以看作常数，d是Hankel矩阵的间隔参数，经过Hankel构造变换后的矩阵维度为(N -d)×(d+1)。利用范德蒙德分解[23]，式(9)矩阵分解如下：

由于式(11)中各项频率项相互独立不同，并且r ≪min{n1,n2}，Sfull是一个低秩矩阵，rank(Sfull)≤P。

除此之外，本文引入离散余弦变换(Discrete Cosine Transform,DCT)分析Sfull的稀疏性，DCT变换是在1974年由Ahmed等人[24]提出的，目前主要应用在图像压缩、音频压缩等领域，是一种类似于傅里叶变换的数学运算。当对Sfull做DCT变换时，变换后的矩阵幅度大小从左到右、从上到下递减，数据主要集中在左上角区域，具有明显的稀疏特性。综合以上分析，取一个距离单元向量做Hankel构造变换后的矩阵具有稀疏-低秩先验特性。除此之外，经过Hankel变换后的矩阵，具有更强的行列相关度。

3 一种基于Hankel-TSPN的三维成像方法

根据第2节的分析，满采样回波的高度-距离切面信号矩阵具有低秩性质，此外，对单个距离单元向量做Hankel构造变换可以增加待恢复矩阵的相关性，并且构造后的矩阵具有稀疏-低秩特性。

为了充分利用上述的低秩和稀疏先验信息，本文介绍一种融合低秩和稀疏先验的基于截断的Schatten-p范数(Truncated Schatten-p Norm,TSPN)的矩阵填充算法，假设S ∈CM×N为一个距离单元向量做Hankel构造变换后的待恢复矩阵，建立优化模型：

其中，M是高度向稀疏采样下观测矩阵，µ>0是稀疏正则化参数，PΩ(·)是投影算子：PΩ(M)m,n=,Ω 是观测矩阵的位置集合。由于低秩问题的求解是NP-hard的，这里采用截断Schatten-p范数[25]逼近秩，相比传统核范数可获得更好的近似效果，优化问题转化为

其中，tr(·)表示矩阵的迹，A,B是由S奇异值分解得到的，[A,Σ,B]=svd(S)，Ar和Br是通过对A,B取前r行得到，以达到截断的效果。将以上问题转化为无约束优化问题，对应的增广拉格朗日函数为

其中，β是惩罚参数，Y,Z为拉格朗日乘子，ρ>1确保惩罚参数递增。采用交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)[26]进行模型求解：

具体求解流程如下：

为了将变量X与变换算子DCT(·)分开，根据Parseval定理和DCT变换性质进行变换处理[27]，对式(17)应用逆变换可得

其中，IDCT(·) 是DCT(·)变换的逆变换，根据广义软阈值算法进行求解，通过已知参数，可得到阈值，如下所示：

对应的广义阈值收缩方程如下：

对式(22)求偏导，令求导结果等于0，可得到Dk+1更新方程：

根据以上分析，得到联合Hankel构造变换的TSPN矩阵填充算法流程，如算法1所示。该算法可以精确地实现回波高度-距离切面信号的重构。

为了对目标实现高分辨三维成像，通过循环方位向对高度-距离切面做Hankel-TSPN，实现对回波的三维重构，得到成像算法流程图见图4。

图4 成像方法流程图Fig.4 Flow diagram of imaging method

本文算法的核心在于利用回波的稀疏-低秩先验特性建立优化模型，通过求解优化问题实现对三维回波重构，重构的三维回波相当于在方位-高度向上以半波长间隔采样，从而提升成像质量。因此，从原理上讲，本文算法同样适用于高度向轨迹稀疏以外的其他缺失场景，如方位向稀疏采样或高度-方位随机缺失采样。

4 仿真与实验验证

本节将通过点目标仿真和实测数据验证上述稀疏轨迹毫米波雷达成像方法的性能，首先通过点目标仿真，初步验证了算法的可行性；进一步采用本文提出算法对多组实测数据进行处理，并与传统算法相比较，论证了算法的优越性。

算法 1 联合Hankel变换的TSPN矩阵填充算法(Hankel-TSPN)Alg.1 TSPN matrix completion algorithm combined with Hankel transformation (Hankel-TSPN)

4.1 仿真实验

为了验证本文方法的有效性，首先设置多点目标仿真，仿真参数如表1所示，点目标布局如图5(a)所示，共有5个点目标位于同一YOZ平面，每个点目标均具有相同的散射系数，采用BP算法处理成像结果见图5。

表1 仿真参数Tab.1 Simulation parameter

图5 多点目标布局及成像图Fig.5 Multi-point target layout and imaging diagram

首先，利用仿真验证Hankel-TSPN算法的优越性，TSPN算法由于结合了稀疏先验，解决了传统矩阵填充无法对行列缺失等结构性缺失矩阵恢复的问题，并且联合Hankel矩阵构造后可以进一步提升对矩阵的重构恢复效果。本文将本方法与一些经典的矩阵填充算法进行了对比。图6比较了奇异值阈值算法(Singular Value Thresholding,SVT),TSPN和Hankel-TSPN对一个高度-距离切面矩阵的恢复效果。从图中可以看出SVT算法由于仅考虑低秩先验无法对行缺失矩阵做恢复，低秩-稀疏先验的TSPN算法可以对行缺失矩阵一定恢复，但恢复效果不佳，而Hankel-TSPN恢复结果与原图基本一致。

图6 不同方法的高度-距离切面矩阵恢复图Fig.6 Different methods of height-distance section matrix restoration

进一步，分别采用30%,20%,10%稀疏轨迹验证本文方法的成像质量，图7给出了不同稀疏度的稀疏轨迹直接3D-BP成像结果和基于Hankel-TSPN的成像结果，图7(a)的轨迹稀疏度为30%；图7(b)轨迹稀疏度为20%；图7(c)轨迹稀疏度为10%。当直接利用稀疏轨迹回波BP成像时，在高度维会出现严重散焦，而采用本文所提出的基于Hankel-TSPN的算法成像时，目标清晰可分辨，且目标能量得到聚焦。

图7 点目标不同稀疏度轨迹下成像结果Fig.7 Imaging results of point targets with different sparsity trajectories

图8中，利用20%稀疏轨迹下单点的成像做高度向剖面图分析，当采用本文提出的方法处理时，目标旁瓣得到了有效的抑制。除此之外，图8中不同算法处理的信号剖面可以用峰值旁瓣比(Peak side lobe ratio,PSLR)和积分旁瓣比(Integral side lobe ratio,ISLR)两个指标进行分析，结果见表2，可以直观地说明采用本节方法处理后的回波成像结果具有更高的聚焦度。

表2 20%稀疏轨迹下高度剖面的峰值旁瓣比与积分旁瓣比(dB)Tab.2 Peak sidelobe ratio and integral sidelobe ratio of height profile at 20% sparse trajectory (dB)

图8 20%稀疏轨迹下中心点高度剖面对比图Fig.8 Height profile comparison of center points under 20%sparse trajectory

4.2 实测实验

为了进一步验证算法的优越性，对实测数据进行成像，目标是一个直径为8 cm的柠檬芯，回波信息通过THz通感一体(Integrated Sensing and Communication at THz band,ISAC-THz)样机获得[28]，图9给出了目标与样机示意图，实验参数见表3。

表3 柠檬芯实测目标参数Tab.3 Lemon core measured target parameter

图9 系统及目标示意图[28]Fig.9 System and target diagram[28]

考虑稀疏轨迹场景，稀疏度为20%的高度向轨迹位置见图10，同一稀疏轨迹采样方式不同方法的成像结果见图11，图中分别比较了理想成像、20%稀疏轨迹BP成像、压缩感知和本文算法的成像结果，其中压缩感知采用的是迭代最小稀疏贝叶斯重构

图10 稀疏度为20%的稀疏轨迹位置示意图Fig.10 Location diagram of sparse trajectory with a sparsity of 20%

图11 柠檬芯成像图Fig.11 Images of lemon core

(Sparsity Bayesian Recovery via Iterative Minimum,SBRIM)[29]，该算法在较多实测数据应用场景下均能有较好的成像质量。除此之外，传统压缩感知方法在对稀疏实测数据处理时，依赖稀疏基的设计，需要场景与目标的先验信息，对未知场景与目标的成像效果较差，而本文所提方法充分利用回波的先验特性，通过对回波先补全再成像的方式，可以保证三维成像精度，由此验证了本文算法的优势。

表4给出了柠檬芯方位-高度切面图的图像质量比较，表中分别从对比度、锐度和图像熵3个角度衡量图像质量，并采用结构相似度(Structural Similarity,SSIM)衡量不同方法在不同采样稀疏度下的成像结果相似度，其中采用本文方法获得的成像结果较于传统的压缩感知算法更好，且更接近满采的成像结果。由此可见，本文方法在对实测数据进行处理时，仍能有很好的成像效果。

表4 三维方位-高度切面图图像质量比较Tab.4 Comparison of image quality of 3D azimuth-height slice

除此之外，为了验证本文方法的鲁棒性，利用高分辨三维毫米波雷达数据集[30]，对几种复杂体目标做了进一步的验证。数据集获取模式为毫米波传感器沿二维平面扫描进行近场成像，扫描尺寸大小为0.4 m×0.4 m，雷达中心频率为78.8 GHz，发射信号带宽为3.6 GHz。不同稀疏度的采样方式见图12，图13给出了不同方法的成像图。

图12 不同稀疏度的稀疏轨迹位置示意图Fig.12 Location diagram of sparse trajectory with different sparsity

图13 不同目标不同稀疏轨迹下成像结果图Fig.13 Imaging results of different targets with different sparse trajectories

通过图13可知，本文方法在仅使用20%～30%的高度向数据时，便能实现接近满采的成像效果，当高度向为稀疏轨迹时，传统BP成像结果中成像目标的对应维度散焦。对不同场景的不同目标做成像，采用本文方法均能获得较高分辨的成像结果，由此体现了本文方法的优越性及广泛的适应性。

最后，本文对该方法的计算量进行了分析，设三维回波矩阵为RK×M×N，K为距离向采样点数，M为方位向采样点数，N为按半波长划分的高度向采样点数，投影矩阵的大小为Dx×Dy×Dz，不同算法的计算复杂度见表5。以柠檬芯实验为例，采用3D-BP成像需要332 s，采用基于SBRIM的压缩感知需要524 s，而采用本文算法需要618 s。因此本文算法和CS,BP在运算量上属于同一量级，但本文算法相比其他算法可以实现更高稀疏度下的高精度成像。

表5 不同算法的计算复杂度Tab.5 Computational complexity of different algorithms

5 结语

传统毫米波雷达为了实现高分辨三维成像，面临着操作复杂度高、数据量大、成本高的问题。本文主要针对，由平台扫描导致的稀疏轨迹构型。该构型通过增大轨迹高度维间距，缩短轨迹扫描时间，降低成像复杂度，然而高度维稀疏会导致成像散焦。本文提出了一种适用于该构型的三维高分辨成像方法。不同于传统的压缩感知算法，本文从回波特性入手，首先分析了回波的低秩-稀疏先验特性，因此构建优化模型并采用ADMM求解得到TSPN算法，进一步引入Hankel矩阵构造，增加待恢复矩阵的行列相关性，得到了一种基于Hankel-TSPN的成像算法，有效解决了传统算法对于单一维度高稀疏度下恢复成像效果不佳的问题。最后，本文通过一系列的实测数据处理验证了在轨迹稀疏度20%以上的场景，采用本文方法均能得到良好的成像结果，并通过与其他算法成像比较，充分体现了本文方法的优越性。