CN113376605A - 一种基于原子范数最小化的mmwave-MIMO雷达超分辨DOA估计方法 - Google Patents

Abstract

本发明公开了一种基于原子范数最小化的mmwave‑MIMO雷达DOA估计方法，首先根据mmwave‑MIMO雷达发射天线个数M，接收天线个数N，信源个数J，获取p×q×(M×N)的三维数据包D；其次，对三维数据包D先在对距离门做q点的傅里叶变换，给出目标所在距离门单元；依次从D中取出目标所在距离门单元下的二维数组B_i，将B_i按子脉冲进行脉冲累积得到单测量矢量b_i；将单测量矢量b_i等价于阵元数为M×N的均匀线阵下的DOA估计；给出基于时分mmwave‑MIMO雷达的超分辨DOA估计函数表达式，并转换成凸函数形式，最后求解原子范数最小化问题；得到目标角度参数

及其导向矢量与噪声子空间的正交性。

Description

一种基于原子范数最小化的mmwave-MIMO雷达超分辨DOA估计 方法

技术领域

本发明实施例涉及mmwave-MIMO雷达信号处理技术，尤其涉及一种基于原子范数最小化的超分辨mmwave-MIMO雷达超分辨DOA估计方法。

背景技术

随着毫米波雷达芯片集成电路和射频技术的逐渐成熟，毫米波多输入多输出(mmwave-MIMO)雷达已经成为智能交通场景下目标检测、目标探测和目标成像等领域的研究热点技术之一。与传统相控阵雷达相比，mmwave-MIMO雷达具有强抗干扰性、高分辨能力、集成度高、体积小、成本廉价等众多优点，因此使得mmwave-MIMO雷达技术的产学研一体化、民用化成为可能。

针对自动化智能交通民用技术领域，将mmwave-MIMO雷达安装至交通路口处的红绿灯杆上是行之有效的，即实现对车辆的测距、测速、波达角估计(DOA)和航迹等数据信息的智能感知。MIMO雷达按收发天线的配置可划分成单基地和双基地MIMO雷达，虽然单基地MIMO雷达和双基地MIMO雷达均能实现目标的高精度定位，但分布式MIMO雷达的设计布置复杂度和天线组件成本比共址MIMO雷达高。同时，MIMO 雷达按波形设计又可划分为时分、频分和码分MIMO雷达，值得注意的是时分方式生成的发射波形被接收天线接收后便于实现信号分离，且时分雷达对硬件要求低、有益于降低雷达造价成本，满足民用mmwave-MIMO雷达的要求。因此基于以上分析，本发明主要关注时分单基地mmwave-MIMO雷达的定位问题，即角度参数估计。

DOA估计是智能交通场景下mmwave-MIMO雷达实现车辆定位不可或缺的关键步骤。基于此，人们提出了大量有效的方法，例如西安电子科技大学在其申请的专利文献“基于雷达通道噪声功率的数字波束形成方法”(授权号CN 109932696 B)中公开了一种数字波束形成(DBF)方法。该方法可以用于解决由雷达通道接收噪声差异性导致DBF合成效率低的问题，从而提高目标检测性能，实现目标的角度定位。但是该方法对目标分辨性能受限于天线孔径，无法实现超分辨目标定位，尤其是在复杂交通堵塞场景下。此外，在文献Y.Zhang,G.Zhang and X.H.Wang,Computationally efficient DOA estimation for monostaticMIMO radar based on covariance matrix reconstruction,Electron.Lett.,49(3),(2017)111-113中提出了一种用于单基地MIMO基于协方差矩阵稀疏表示的DOA估计方法。该方法可以实现独立信源下的超分辨DOA估计，但是针对时分单基地MIMO雷达，由于接收回波信号是相干的，会导致协方差矩阵秩亏，进而错误估计目标角度。所以上述方法需要结合空间平滑算法以解相干，不可避免地，空间平滑是以牺牲孔径大小为代价，因此该方法在时分雷达上的应用必然导致目标分辨性能的降低。

发明内容

发明目的：本发明的目的是提出了一种基于原子范数最小化(ANM)的 mmwave-MIMO雷达DOA估计方法，该方法可以在不损失孔径大小的前提下，通过凸优化求解接收回波信号的原子范数最小化问题，高精度实现对多个相干信源的超分辨 DOA估计。

技术方案：本发明的一种基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，包括以下步骤：

S1、假设已知mmwave-MIMO雷达发射天线个数为M，接收天线个数为N，信源个数为J，其中J<M×N；对发射天线端采用时分方式发射正交波形，依次从第1根发射天线开始到第M根发射天线结束，每间隔T_c秒发射线性调频信号，T_c为一个脉冲周期，在一个脉冲周期内，发射p个子脉冲进行脉冲积累；并且对子脉冲做距离门上的q次离散采样，最终得到p×q×(M×N)的三维数据包D；

S2：对三维数据包D先取出距离门数据d＝D(p,:,M×N)，对d做q点的傅里叶变换，找出相对应的I个峰值点，即为目标所在距离门单元[q₁,q₂,…,q_I]，I≤J；不妨令循环变量i＝1,…,J，q_i为第i个目标所在距离门单元；若每个目标均对应不同的距离门单元，则 I＝J；反之，I<J；随后开始循环，初始化i＝1；

S3、判断循环变量i是否小于等于I，若是，则执行步骤S4；否则，则退出循环，执行步骤S9；

S4、从三维数据包D中取出目标所在距离门单元下的二维数组B_i，其中 B_i＝D(1:p,q_i,1:M×N)，将二维数组B_i按子脉冲进行脉冲累积得到单测量矢量b_i；

S5、将单测量矢量b_i等价于阵元数为M×N的均匀线阵下的DOA估计；

S6、根据原子范数理论，首先给出基于时分mmwave-MIMO雷达的超分辨DOA估计函数表达式，其次根据半正定规划理论将函数表达式转换成凸函数形式，最后通过 CVX凸优化工具包来求解原子范数最小化问题；

S7、对根据CVX凸优化工具包所求解得出的

进行特征值分解(EVD)，最终通过谱峰搜索得到目标角度参数

及其导向矢量与噪声子空间的正交性；

S8、令i＝i+1，返回执行步骤S3；

S9、结束循环，此时已经遍历目标所在距离门，整合得出Θ＝[Θ₁,...,Θ_I]及其每个距离门下的P_power。

进一步的，步骤S1中mmwave-MIMO雷达发射天线和接收天线别采用等间距6λ和1.5λ布局，λ为60GHz毫米波的波长。

进一步的，步骤S4中得到的单测量矢量b_i为：

其中，pp为取值从1到p的累加中间变量，

为复数集合。

进一步的，步骤S5具体为：

不妨假设第i个距离门单元上有x个目标，因此转换成以下带噪声的信号模型：

b_i＝A_is+n＝z_i+n (2)；

其中，

为由x个导向矢量张成的矩阵，

和

分别为第1列、第2列和第x列的导向矢量，

j为单位虚数，θ_x为目标所在角度，上标T为转置算子，

为复数集合，z_i＝A_is为阵元上的无噪接收信号，s为发射信号，n为阵元接收到的噪声，服从均值为0，协方差矩阵为σ²E，即：

n～N(0,σ²E) (3)；

其中，σ²为高斯白噪声的功率，E表示单位矩阵。

进一步的，步骤S6具体为：

基于时分mmwave-MIMO雷达的超分辨DOA估计函数表达式为：

其中，

为z_i的估计值，

为z_i的原子范数，数学符号s.t.指的是受限于，ξ为噪声功率的上限，因为

是凸的，根据半正定规划理论将公式(4)转换成以下凸函数形式：

其中，T(μ)、y、t分别为凸优化转换生成的中间矩阵、向量和变量，

分别对应为T(μ)、y、t的估计，||b_i-z_i||₂为残差向量b_i-z_i的l₂范数，ε为正则化参数，上标H为共轭转置算子，约束条件

表明矩阵满足半正定，理论上T(μ)是稀疏的并且rank(T(μ))＝x；T(μ)表示由向量

决定的Hermitian Toeplitz矩阵，其中μ₁为μ中第1个元素，μ_MN为μ中第MN个元素，

为复数集合，T(μ) 具体表现形式如下：

其中，

为μ中第2个元素的共轭转置、

为μ中第M个元素的共轭转置、

为μ中第M-1个元素的共轭转置、μ_M-1为μ中第M-1个元素；

最后通过CVX凸优化工具包来求解原子范数最小化问题。

进一步的，步骤S7具体为：

对根据CVX工具包所求解得出的

进行特征值分解：

其中

为正交矩阵，

为奇异值矩阵，

为正交矩阵的共轭转置；令

为噪声子空间，最终通过谱峰搜索得到目标角度参数

其中，目标角度参数

为真实角度θ的估计，P_power反映了导向矢量与噪声子空间的正交性，a(θ)＝[1,e^-j3πsinθ,...,e^{-j3π(MN-1)sinθ}]^T，a(θ)^H为a(θ)的共轭转置，

为

的共轭转置，由于已知该距离门上有x个目标，因此最终会得到x个不同的角度，不妨定义为

为第一个目标的角度，

第x个目标的角度。

有益效果：与现有技术相比，本发明通过引入原子范数最小化问题，实现了目标超分辨DOA估计。以3发4收60GHz毫米波雷达为例，在微波暗室环境下，传统DBF 的角度分辨率无法突破2.86°，而本发明所提方法仍可以对角度范围差小于2.86°的两个角反射器进行目标分辨。本方法通过利用接收数据的先验稀疏性，可以实现高精度DOA 估计。本发明方法用以在保证高精度估计性能的同时，提升相干信号下角度分辨率。实测结果表明本发明所提方法确实在在一定程度上提高了角度分辨率，达到了超分辨DOA 估计的效果。

附图说明

图1是本发明法方法流程图；

图2是本发明的3发4收mmwave-MIMO雷达收发天线及其等效虚拟天线示意图；

图3是本发明的两个角反射器实际放置场景，其中，(a)为角反射器和 mmwave-MIMO雷达系统探测场景，(b)为DBF角度分辨率场景；

图4是本发明针对图3所得实测数据DBF效果图，其中(a)为距离-方位二维图；(b)为谱峰搜索图；

图5是本发明针对图3所得实测数据ANM的二维图，其中(a)为距离-方位二维图，(b)为谱峰搜索图；

图6是针对单个角反射器实测数据所得的三维效果图，其中(a)为基于DBF的距离-方位-正交性三维效果图；(b)为基于ANM的距离-方位-正交性三维效果图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

如图1所示，本发明的一种基于原子范数最小化(ANM)的mmwave-MIMO雷达DOA 估计方法，包括以下步骤：

S1、假设已知mmwave-MIMO雷达发射天线个数为M，接收天线个数为N，信源个数为J，其中J<M×N。MIMO发射天线和接收天线别采用等间距6λ和1.5λ布局，λ为60GHz毫米波的波长。本发明对发射天线端采用时分方式发射正交波形，依次从第1 根发射天线开始到第M根发射天线结束，每间隔T_c秒(一个脉冲周期)发射线性调频信号。采用每次仅开启一个天线的目的是保证发射波形在时域正交。在一个脉冲周期内，发射p个子脉冲进行脉冲积累。同时，N个接收天线会依次接收M根发射天线发出的电磁波信号，最终获取到M×N个虚拟天线的回波信号。以M＝3，N＝4为例，图2为 mmwave-MIMO雷达收发天线及其等效虚拟天线示意图。若对子脉冲做距离门上的q次离散采样，最终得到p×q×(M×N)的三维数据包D。

S2、对三维数据包D先取出距离门数据d＝D(p,:,M×N)，对d做q点的傅里叶变换(FFT)，找出相对应的I个峰值点，即为目标所在距离门单元[q₁,q₂,…,q_I]，I≤J。不妨令循环变量i＝1,…,J，q_i为第i个目标所在距离门单元。若每个目标均对应不同的距离门单元，则I＝J；反之，I<J。随后开始循环，初始化i＝1。

S3、判断循环变量i是否小于等于I，若是，则执行步骤S4；否则，则退出循环，执行步骤S9；

S4、从三维数据包D中取出目标所在距离门单元下的二维数组B_i，其中 B_i＝D(1:p,q_i,1:M×N)。为了提高系统的抗噪性，将二维数组B_i按子脉冲进行脉冲累积得到单测量矢量b_i：

其中，pp为取值从1到p的累加中间变量，

为复数集合。

并且将数据的从二维转换成一维，降低了计算复杂度，便于硬件系统实时迅速处理信号，快速实现自动化目标定位。

S5、根据阵列信号处理相关理论知识，将单测量矢量b_i等价于阵元数为M×N的均匀线阵下的DOA估计。不妨假设第i个距离门单元上有x个目标，因此转换成以下带噪声的信号模型：

b_i＝A_is+n＝z_i+n (2)；

其中，

为由x个导向矢量张成的矩阵，

和

分别为第1列、第2列和第x列的导向矢量，

j为单位虚数，θ_x为目标所在角度，上标T为转置算子，

为复数集合，z_i＝A_is为阵元上的无噪接收信号，s为发射信号，n为阵元接收到的噪声，服从均值为0，协方差矩阵为σ²E，即：

n～N(0,σ²E) (3)；

其中，σ²为高斯白噪声的功率，E表示单位矩阵。

S6、借助文献(B.N.Bhaskar,G.Tang,B.Recht.Atomic norm denoising withapplications to line spectral estimation.IEEE Trans.Signal Process.,2013,61(23):5987-5999) 中所提的原子范数理论，首先给出基于时分mmwave-MIMO雷达的超分辨DOA估计函数表达式：

其中，

为z_i的估计值，

为z_i的原子范数，数学符号s.t.指的是受限于，ξ为噪声功率的上限，因为

是凸的，根据半正定规划理论将公式(4)转换成以下凸函数形式：

其中，T(μ)、y、t分别为凸优化转换生成的中间矩阵、向量和变量，

分别对应为T(μ)、y、t的估计，||b_i-z_i||₂为残差向量b_i-z_i的l₂范数，ε为正则化参数，上标H为共轭转置算子，约束条件

表明矩阵满足半正定，理论上T(μ)是稀疏的并且rank(T(μ))＝x；T(μ)表示由向量

决定的Hermitian Toeplitz矩阵，其中μ₁为μ中第1个元素，μ_MN为μ中第MN个元素，

为复数集合， T(μ)具体表现形式如下：

其中，

为μ中第2个元素的共轭转置、

为μ中第M个元素的共轭转置、

为μ中第M-1个元素的共轭转置、μ_M-1为μ中第M-1个元素。

本实施例通过CVX凸优化工具包(M.Grant and S.Boyd,“CVX:Matlab softwarefor disciplined convex programming,version 2.0beta,”Available online athttp://cvxr.com/cvx, Sep.2013.)来求解原子范数最小化问题。

S7、对根据CVX工具包所求解得出的

进行特征值分解(EVD)：

其中，

为正交矩阵，

为奇异值矩阵，

为正交矩阵的共轭转置；令

为噪声子空间，最终通过谱峰搜索得到目标角度参数

其中，目标角度参数

为真实角度θ(目标所在真实方位)的估计，P_power反映了导向矢量与噪声子空间的正交性，a(θ)＝[1,e^-j3πsinθ,...,e^{-j3π(MN-1)sinθ}]^T，a(θ)^H为a(θ)的共轭转置，

为

的共轭转置，一般的，从-90°到90°的方位范围内平均选取900个方位网格，即每隔0.2°进行谱峰搜索。由于已知该距离门上有x个目标，因此最终会得到x个不同的角度，为了方便书写，不妨定义为

为第一个目标的角度，

第x个目标的角度。图5(b)绘出了在目标所在距离门单元上的谱峰搜索图，可以明显看出该距离门上有两个亮点，即表示目标所在距离和角度信息。

S8、令i＝i+1，返回执行步骤S3。

S9、结束循环，此时已经遍历目标所在距离门，整合得出Θ＝[Θ₁,...,Θ_I]及其每个距离门下的P_power，绘出如图5(a)所示的基于ANM的实测数据二维图，横坐标表示角度网格数，纵坐标表示距离门单元，像素大小(0-1之间)表示导向矢量与噪声子空间的正交性。

本发明提出的一种基于原子范数最小化的mmwave-MIMO雷达超分辨DOA估计方法，针对图3所示的实际场景，本发明验证了所提方法的有效性，同时图4和图5分别画出了基于DBF和ANM的二维图和谱峰搜索图，其中图4中(a)和图5中(a)均为距离- 方位二维图，横坐标为方位网格，纵坐标为距离门单元，所在点的亮度表示是否存在目标，若亮度为0即黑色表示此处无目标，反之亦然；图4中(b)和图5中(b)均为谱峰搜索图，横坐标表示方位网格，纵坐标为正交性，曲线的峰值表示目标所在方位；对比基于DBF的二维图无法分辨距离较近(角度差小于2.86°)的两个角反射器，而本发明所提出的ANM可以分辨两个相距近的目标，实现了超分辨目标定位。图6绘出了针对单个角反射器的DOA估计对比图，其中图6中(a)为基于DBF的距离-方位-正交性三维效果图，x轴为方位网格，y轴为距离门单元，z轴为正交性，三维图的峰值所在点表示目标所在距离和方位信息；图6中(b)为基于ANM的距离-方位-正交性三维效果图，定义参照基于DBF的距离-方位-正交性三维效果图；虽然本发明所提方法和DBF所估计的目标角度一致，但是本发明利用了先验稀疏性，所以主瓣更窄，几乎没有副瓣，而DBF 主瓣很宽，受到了地杂波的干扰出现了多个副瓣。

Claims (6)

1.一种基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，其特征在于，包括以下步骤：

S1、假设已知mmwave-MIMO雷达发射天线个数为M，接收天线个数为N，信源个数为J，其中J<M×N；对发射天线端采用时分方式发射正交波形，依次从第1根发射天线开始到第M根发射天线结束，每间隔T_c秒发射线性调频信号，T_c为一个脉冲周期，在一个脉冲周期内，发射p个子脉冲进行脉冲积累；并且对子脉冲做距离门上的q次离散采样，最终得到p×q×(M×N)的三维数据包D；

S2：对三维数据包D先取出距离门数据d＝D(p,:,M×N)，对d做q点的傅里叶变换，找出相对应的I个峰值点，即为目标所在距离门单元[q₁,q₂,…,q_I]，I≤J；不妨令循环变量i＝1,…,J，q_i为第i个目标所在距离门单元；若每个目标均对应不同的距离门单元，则I＝J；反之，I<J；随后开始循环，初始化i＝1；

S3、判断循环变量i是否小于等于I，若是，则执行步骤S4；否则，则退出循环，执行步骤S9；

S4、从三维数据包D中取出目标所在距离门单元下的二维数组B_i，其中B_i＝D(1:p,q_i,1:M×N)，将二维数组B_i按子脉冲进行脉冲累积得到单测量矢量b_i；

S5、将单测量矢量b_i等价于阵元数为M×N的均匀线阵下的DOA估计；

S6、根据原子范数理论，首先给出基于时分mmwave-MIMO雷达的超分辨DOA估计函数表达式，其次根据半正定规划理论将函数表达式转换成凸函数形式，最后通过CVX凸优化工具包来求解原子范数最小化问题；

S7、对根据CVX凸优化工具包所求解得出的

进行特征值分解(EVD)，最终通过谱峰搜索得到目标角度参数

及其导向矢量与噪声子空间的正交性；

S8、令i＝i+1，返回执行步骤S3；

S9、结束循环，此时已经遍历目标所在距离门，整合得出Θ＝[Θ₁,...,Θ_I]及其每个距离门下的P_power。

2.根据权利要求1所述的基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，其特征在于，步骤S1中mmwave-MIMO雷达发射天线和接收天线别采用等间距6λ和1.5λ布局，λ为60GHz毫米波的波长。

3.根据权利要求1所述的基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，其特征在于，步骤S4中得到的单测量矢量b_i为：

其中，pp为取值从1到p的累加中间变量，

为复数集合。

4.根据权利要求1所述的基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，其特征在于，步骤S5具体为：

不妨假设第i个距离门单元上有x个目标，因此转换成以下带噪声的信号模型：

b_i＝A_is+n＝z_i+n (2)；

其中，

为由x个导向矢量张成的矩阵，

和

分别为第1列、第2列和第x列的导向矢量，

j为单位虚数，θ_x为目标所在角度，上标T为转置算子，

为复数集合，z_i＝A_is为阵元上的无噪接收信号，s为发射信号，n为阵元接收到的噪声，服从均值为0，协方差矩阵为σ²E，即：

n～N(0,σ²E) (3)；

其中，σ²为高斯白噪声的功率，E表示单位矩阵。

5.根据权利要求1所述的基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，其特征在于，步骤S6具体为：

基于时分mmwave-MIMO雷达的超分辨DOA估计函数表达式为：

其中，

为z_i的估计值，

为z_i的原子范数，数学符号s.t.指的是受限于，ξ为噪声功率的上限，因为

是凸的，根据半正定规划理论将公式(4)转换成以下凸函数形式：

其中，T(μ)、y、t分别为凸优化转换生成的中间矩阵、向量和变量，

分别对应为T(μ)、y、t的估计，||b_i-z_i||₂为残差向量b_i-z_i的l₂范数，ε为正则化参数，上标H为共轭转置算子，约束条件

表明矩阵满足半正定，理论上T(μ)是稀疏的并且rank(T(μ))＝x；T(μ)表示由向量

决定的Hermitian Toeplitz矩阵，其中μ₁为μ中第1个元素，μ_MN为μ中第MN个元素，

为复数集合，T(μ)具体表现形式如下：

其中，

为μ中第2个元素的共轭转置、

为μ中第M个元素的共轭转置、

为μ中第M-1个元素的共轭转置、μ_M-1为μ中第M-1个元素；

最后通过CVX凸优化工具包来求解原子范数最小化问题。

6.根据权利要求1所述的基于原子范数最小化的mmwave-MIMO雷达DOA估计方法，其特征在于，步骤S7具体为：

对根据CVX工具包所求解得出的

进行特征值分解：

其中

为正交矩阵，

为奇异值矩阵，

为正交矩阵的共轭转置；令

为噪声子空间，最终通过谱峰搜索得到目标角度参数

其中，目标角度参数

为真实角度θ的估计，P_power反映了导向矢量与噪声子空间的正交性，

a(θ)^H为a(θ)的共轭转置，

为

的共轭转置，由于已知该距离门上有x个目标，因此最终会得到x个不同的角度，不妨定义为

为第一个目标的角度，

第x个目标的角度。

Images