CN105093185A - 基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法 - Google Patents

Abstract

本发明提供一种基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法，首先建立单基地MIMO雷达系统的接收信号模型，构造降维转换矩阵进行降维处理；然后利用酉变换矩阵将降维后的接收数据矩阵变为实域的，设计实值扩展数据矩阵并获得其协方差矩阵；根据Khatri-Rao积，将实值协方差矩阵向量化以解决多测量矢量(MMV)问题，并得到稀疏表示框架下的相应模型；最后设计权值矩阵获得估计参数并构造实值l₁范数最小化框架，得到恢复矩阵，寻找恢复矩阵中的非零行，实现对MIMO雷达系统中目标DOA的估计。本发明计算复杂度明显降低，具有更高的角度分辨率和更好的角度估计性能，并且具有最低的SNR临界值。

Description

基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法

技术领域

本发明涉及MIMO雷达系统技术领域，特别涉及MIMO雷达系统DOA估计的应用，具体说是一种基于高运算效率稀疏表示的单基地MIMO雷达DOA估计方法。

背景技术

多输入多输出(MIMO)雷达在雷达领域引起了很大的关注并且成为了一个热门研究课题(IEEESignalProcessingMagazine，2007，24(5)：106-114)。MIMO雷达可以分为两类，一种是统计MIMO雷达，另外一种是相干MIMO雷达包括双基地MIMO雷达和单基地MIMO雷达。在双基地MIMO雷达中，发射阵列和接收阵列是彼此分开的，但是在单基地中相反，它们是相互紧密的。因此，双基地MIMO雷达波离角(DOD)和波达角(DOA)的方向是不同的，而在单基地MIMO中它们是相同的。本发明研究的是单基地MIMO雷达角度估计。

在由接收阵列中的噪声而被损坏的测量值中得到的多目标波达方向(DOAs)估计是阵列信号处理和MIMO雷达实际应用中的重要方面，已经提出的一些MIMO雷达角度估计算法，有Capon算法，多重信号分类(MUSIC)算法，旋转不变子空间(ESPRIT)算法(ElectronicsLetters，2008，44(12)：770-771)和基于张量分析的算法(IEEETransationsonSignalProcessing，2010，58(11)：5693-5705)等。根据均匀线阵单基地MIMO雷达虚拟阵列的特殊结构，人们提出了RD-ESPRIT(ElectronicsLetters，2011，47(4)：283-284)和RD-Capon的DOA估计算法，提高了角度估计性能，尤其是在低SNR情况下。利用非圆信号的特点而增大单基地MIMO雷达的孔径提出了一种共轭ESPRIT角度估计算法，比传统子空间算法提供了更好的角度估计性能。另一方面，最近稀疏表示领域的出现获得了极大的关注，给DOA估计提供了新视点。人们已经提出了基于稀疏表示的DOA估计方法，例如l₁-SVD(IEEETransationsonSignalProcessing，205，53(8)：3010-3022)，l₁-SRACV(IEEETransationsonSignalProcessing，2011，59(2)：629-638)和CMSR等算法。所有的仿真结果证明，与传统基于子空间的方法相比，这些算法能够更好地适应很多情况，例如连续目标，低SNR和有限快拍数。此外，一种实域l₁-SVD方法与l₁-SVD算法相比具有更低的计算复杂度和较好的角度估计性能。然而，在单基地MIMO雷达中基于稀疏表示的方法大多数有两个缺点：1)在信号重建过程中，由于高维完备字典和多测量矢量(MMV)问题造成了高计算量。2)MIMO雷达孔径的原因，使得角度分辨率是有限的。

发明内容

本发明的目的在于提出一种计算复杂度明显降低，具有更高的角度分辨率和更好的角度估计性能，并且具有最低的SNR临界值的基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法

本发明的目的是这样实现的：

1、建立单基地MIMO雷达的接收信号模型，并设计降维矩阵对接收数据进行降维处理。

考虑一个由M个发射天线和N个接收天线组成的窄带单基地MIMO雷达，其发射阵列和接收阵列均为阵元距为半个波长的均匀线阵，并且所有元素都是正交的。在发射阵列，M个发射天线发射M个具有相同带宽和中心频率的正交窄带信号。假设有P个远场不相关目标，θ_p表示第p个目标关于发射阵列和接收阵列的DOA，则接收阵列匹配滤波器的输出可以表示为

$x\left(t\right)=\[a_t\left({\mathrm{\θ}}_1\right)\⊗a_r\left({\mathrm{\θ}}_1\right),...,a_t\left({\mathrm{\θ}}_p\right)\⊗a_r\left({\mathrm{\θ}}_p\right)\]s\left(t\right)+n\left(t\right)---\left(1\right)$

其中是接收数据矢量，是零均值随机复高斯信号数据矢量，中β_p(t)和f_p(t)分别是反射系数和多普勒频率。其中表示Kronecker积操作，是一个具有零均值并且协方差矩阵为σ²I_MN的随机复高斯白噪声矢量。发射导向矢量a_t(θ_p)和接收导向矢量a_r(θ_p)分别表示为

$a_t\left({\mathrm{\θ}}_p\right)=\left\{\begin{array}{c}\[e^{-j\mathrm{\π}\frac{M-1}{2}{\mathrm{sin\θ}}_p},...,e^{-j\mathrm{\π}\frac{1}{2}{\mathrm{sin\θ}}_p},e^{-j\mathrm{\π}\frac{1}{2}{\mathrm{sin\θ}}_p},...,e^{-j\mathrm{\π}\frac{M-1}{2}{\mathrm{sin\θ}}_p}\],M=2k\\ \[e^{-j\mathrm{\π}\frac{M-1}{2}{\mathrm{sin\θ}}_p},...,1,...,e^{-j\mathrm{\π}\frac{M-1}{2}{\mathrm{sin\θ}}_p}\],M=2k+1\end{array}\right.$

和

$a_r\left({\mathrm{\θ}}_p\right)=\left\{\begin{array}{c}\[e^{-j\mathrm{\π}\frac{N-1}{2}{\mathrm{sin\θ}}_p},...,e^{-j\mathrm{\π}\frac{1}{2}{\mathrm{sin\θ}}_p},e^{-j\mathrm{\π}\frac{1}{2}{\mathrm{sin\θ}}_p},...,e^{-j\mathrm{\π}\frac{N-1}{2}{\mathrm{sin\θ}}_p}\],N=2k\\ \[e^{-j\mathrm{\π}\frac{N-1}{2}{\mathrm{sin\θ}}_p},...,1,...,e^{-j\mathrm{\π}\frac{N-1}{2}{\mathrm{sin\θ}}_p}\],N=2k+1\end{array}\right.$

通过收集J个快拍，定义 $A=\[a_t\left({\mathrm{\θ}}_1\right)\⊗a_r\left({\mathrm{\θ}}_1\right),...,a_t\left({\mathrm{\θ}}_p\right)\⊗a_r\left({\mathrm{\θ}}_p\right)\],$ 则

X＝AS+N(2)

其中是接收数据矩阵，S＝[s(t₁),s(t₂),...,s(t_J)]和N＝[n(t₁),n(t₂),...,n(t_J)]分别是信号数据矩阵和噪声矩阵。

根据单基地MIMO雷达的特殊结构，只有Q＝M+N-1个不同元素，因此导向矢量可以表示为

$a_t\left(\mathrm{\θ}\right)\⊗a_r\left(\mathrm{\θ}\right)=G_{MN\×Q}b\left(\mathrm{\θ}\right)---\left(3\right)$

其中G_MN×Q和b(θ)是降维转换矩阵和导向矢量，分别为

$G_{MN\×Q}=\left[\begin{array}{c}{J}_0\\ J_1\\ .\\ .\\ .\\ J_{M-1}\end{array}\right]---\left(4\right)$

$b\left(\mathrm{\θ}\right)=\left\{\begin{array}{c}\[e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}},...,e^{-j\mathrm{\π}\frac{1}{2}\sin \mathrm{\θ}},e^{-j\mathrm{\π}\frac{1}{2}\sin \mathrm{\θ}},...,e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}}\],Q=2k\\ \[e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}},...,1,...,e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}}\],Q=2k+1\end{array}\right.---\left(5\right)$

其中J_m＝[0_N×m,I_N,0_N×(M-m-1)]，m＝0,1,...,M-1。根据公式(4)，我们定义矩阵F＝G^HG，如下所示

根据噪声矩阵N的特征，降维矩阵定义为W＝F^-(1/2)G^H，满足WW^H＝I_M+N-1。将接收信号X与W相乘，有

$Y=F^{(1/2)}BS+WN=\hat{B}S+WN---\left(7\right)$

其中B＝[b(θ₁),b(θ₂),...,b(θ_p)]，

2、利用酉变换矩阵，将降维后的接收数据矩阵变成实域，构造实值扩展数据矩阵并获得其协方差矩阵。

将降维后的接收数据矩阵变成实域的，如下所示

由公式(7)可知，利用降维转换后，接收数据对应一个具有权值矩阵F^(1/2)的线性阵列。另一方面，导向矩阵变成了且满足

${\mathrm{\Π}}_{M+N-1}{\hat{B}}^{*}=\hat{B}---\left(8\right)$

其中Π_K表示具有反对角元素为1，其他元素为0的K×K交换矩阵，(·)^*表示共轭操作。根据公式(8)可以得出结论，在降维转换后线性阵列也是一个中心对称阵列。因此，酉变换能够使复值导向矩阵变成实值的。酉变换矩阵定义为

$U_{2K}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_K& {\mathrm{jI}}_K\\ {\mathrm{\Π}}_K& -{\mathrm{j\Π}}_K\end{array}\right]---\left(9\right)$

和

$U_{2K+1}=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_K& 0& {\mathrm{jI}}_K\\ 0^T& \sqrt{2}& 0^T\\ {\mathrm{\Π}}_K& 0& -{\mathrm{j\Π}}_K\end{array}\right]---\left(10\right)$

然后把酉变换矩阵U_Q利用到接收数据Y上，新的接收数据矩阵可以表示为

$\hat{Y}=U_Q^{H}Y=U_Q^H\hat{B}S+U_Q^{H}WN---\left(11\right)$

其中(·)^H表示共轭转置。根据公式(11)，经过酉变换后，复导向矩阵转换成为实值导向矩阵 $\stackrel{\‾}{B}=U_Q^H\hat{B}.$

构造实值扩展数据矩阵，如下所示

数据矩阵可以分为实数部分和虚数部分，有

$\mathrm{Re}\left(\hat{Y}\right)=\stackrel{\‾}{B}\mathrm{Re}\left(S\right)+\mathrm{Re}\left(U_Q^{H}WN\right)---\left(12\right)$

$\mathrm{Im}\left(\hat{Y}\right)=\stackrel{\‾}{B}\mathrm{Im}\left(S\right)+\mathrm{Im}\left(U_Q^{H}WN\right)---\left(13\right)$

其中Re(·)和Im(·)分别表示实数部分和虚数部分，利用公式(12)和(13)，可以构造实值扩展数据矩阵，表示为

$\stackrel{\‾}{Y}=\[\mathrm{Re}\left(\hat{Y}\right),\mathrm{Im}\left(\hat{Y}\right)\]=\stackrel{\‾}{B}\stackrel{\‾}{S}+\stackrel{\‾}{N}---\left(14\right)$

其中 $\stackrel{\‾}{S}=\[\mathrm{Re}\left(S\right),\mathrm{Im}\left(S\right)\]$ 和 $\stackrel{\‾}{N}=\[\mathrm{Re}\left(U_Q^{H}WN\right),\mathrm{Im}\left(U_Q^{H}WN\right)\].$

获得扩展数据矩阵的实值协方差矩阵，如下所示

$\stackrel{\‾}{Y}$ 的实值协方差矩阵表示为

$R_{\stackrel{\‾}{Y}}=E\[\stackrel{\‾}{Y}{\stackrel{\‾}{Y}}^H\]=\stackrel{\‾}{B}{R}_S{\stackrel{\‾}{B}}^H+{\mathrm{\σ}}^2{I}_Q---\left(15\right)$

其中E[·]表示期望，σ²是噪声功率。

3、将实值协方差矩阵向量化以解决MMV问题，并获得稀疏表示框架下的相应模型。

针对MMV问题，将实值协方差矩阵向量化，如下所示

基于Khatri-Rao积操作，对进行向量化操作，有

其中⊙表示Khatri-Rao积操作，d＝vec(R_S)。由公式(16)可知，数据向量y可以看做是具有单快拍虚拟阵列输出的信号模型，因此，带有MMV问题的稀疏信号重构可以转变为SMV问题。另外，在此信号模型中，类似于公式(11)中的实值导向矩阵，公式(16)中的操作极大地增加了虚拟自由度，并且比公式(2)中的信号模型拥有更大的孔径。因此，角度分辨率得到很大的提高。

获得稀疏表示框架下的相应模型，如下所示

当整个空间的离散样本网格数比目标数目大很多时，目标的DOA是稀疏的，公式(2)中X的奇异值分解(SVD)可以表示为

$X=U_s{\mathrm{\Λ}}_s{U}_s^H+U_n{\mathrm{\Λ}}_n{U}_n^H---\left(17\right)$

其中和是由P个最大奇异值对应的奇异向量所组成。和是由剩下的MN-P个奇异值对应的奇异向量所组成。Λ_s和Λ_n分别是对角线元素为P个最大奇异值和对角线元素为剩下MN-P个奇异值的对角矩阵。将接收信号X与V_s相乘，有

X_SV＝AS_SV+N_SV(18)

其中X_SV＝XV_s，S_SV＝SV_s，N_SV＝NV_s。利用SVD技术后，接收数据X_SV比X具有更低的维数。利用样本网格公式(18)能够被稀疏表示，其中L≥P。然后发射完备字典和接收完备字典可以分别表示为 $A_t^{\widehat{\mathrm{\θ}}}=\[a_t\left({\widehat{\mathrm{\θ}}}_1\right),...,a_t\left({\widehat{\mathrm{\θ}}}_L\right)\]$ 和 $A_r^{\widehat{\mathrm{\θ}}}=\[a_r\left({\widehat{\mathrm{\θ}}}_1\right),...,a_r\left({\widehat{\mathrm{\θ}}}_L\right)\].$ 于是构造完备字典在稀疏表示框架下，公式(18)能够写成为

$X_{SV}=A_{\widehat{\mathrm{\θ}}}{S}_{SV}^{\widehat{\mathrm{\θ}}}+N_{SV}---\left(19\right)$

其中和S_SV有相同的行支持，即矩阵是稀疏的。为了估计可以将公式(19)构造成为l₁范数最小化问题，如下所示

$\min \left|\right|{\left(S_{SV}^{\widehat{\mathrm{\θ}}}\right)}^{\left(l_2\right)}|{|}_1,s.t.||X_{SV}-A_{\widehat{\mathrm{\θ}}}{S}_{SV}^{\widehat{\mathrm{\θ}}}|{|}_2^2\≤\stackrel{\‾}{\mathrm{\η}}---\left(20\right)$

其中||·||₁和||·||₂分别表示l₁范数和l₂范数。表示第q个元素等于中第q行l₂范数的列向量，是正则化参数。获得稀疏矩阵后，通过找到中的非零行能够对DOA进行估计。

类似公式(19)的稀疏表示框架，基于公式(16)中的信号模型，构造实值稀疏表示DOA估计方法。让表示整个空间域所有感兴趣DOA的离散采样网格。可能的DOA数将比目标数大很多，即L＞＞P。然后DOA估计的一维完备字典可以构造成其中且在稀疏表示框架下，公式(16)可以被完备字典代替，如下所示

$y={\underset{\‾}{B}}_{\widehat{\mathrm{\θ}}}{d}_{\widehat{\mathrm{\θ}}}+{\mathrm{\σ}}^{2}vec\left(I_Q\right)---\left(21\right)$

其中是稀疏向量，即其只有P个非零元素，冗杂元素为0。因此通过找到中的非零元素能够对DOA进行估计。

4、设计权值矩阵获得估计参数并构造实值l₁范数最小化框架，利用编程软件包SOC(二阶锥)计算方法，获得恢复矩阵，并寻找恢复矩阵中的非零行，实现对单基地MIMO雷达目标DOA的估计。

设计权值矩阵，如下所示

实际中，估计未知协方差矩阵可以通过得到。然后有和其中Δy是估计误差，接着结合公式(21)中的稀疏表示框架和估计误差Δy，具有SMV的实值l₁范数最小化问题构造如下

$\min \left|\right|d_{\widehat{\mathrm{\θ}}}|{|}_1,s.t.||\hat{y}-{\underset{\‾}{B}}_{\widehat{\mathrm{\θ}}}{d}_{\widehat{\mathrm{\θ}}}-{\mathrm{\σ}}^{2}vec\left(I_Q\right)|{|}_2^2\≤\mathrm{\η}---\left(22\right)$

其中η是正则化参数，η设置误差量并且在最后的DOA估计性能中起到重要的作用。在公式(22)的l₁范数最小化问题中，参数η的选取依靠Δy的分布。

根据W和U_Q的定义，他们是正交矩阵。因此，和具有零均值复高斯分布，这是因为高斯随机矩阵的正交不变性，其分布不受乘法限制，其中[·]_i,j表示矩阵的第(i,j)个元素。此外，由于和的实部和虚部是零均值实高斯分布，和均满足零均值高斯分布。因此，可以得出结论，估计误差满足

$\mathrm{\Δ}y=vec({\hat{R}}_{\stackrel{\‾}{Y}}-R_{\stackrel{\‾}{Y}})~AsN(0,\frac{1}{J}{R}_{\stackrel{\‾}{Y}}^T\⊗R_{\stackrel{\‾}{Y}})---\left(23\right)$

由公式(23)可知，利用降维转换和酉变换后，估计误差Δy满足其中AsN(μ,σ²)表示具有均值μ和协方差σ²的渐近正态分布。注意到Δy不是渐近标准正态分布，参数η不容易被估算出来。

构造权值矩阵使得估计误差Δy满足渐近标准正态分布，于是有

$\left|\right|{\stackrel{\‾}{W}}^{-\frac{1}{2}}\mathrm{\Δ}y|{|}_2^2~{\mathrm{As\χ}}^2(0,I_{Q^2})---\left(24\right)$

其中表示具有Q²自由度的渐近卡方分布。

构造实值l₁范数最小化框架，如下所示

结合公式(22)和公式(24)，构造具有SMV的实值l₁范数最小化问题，如下所示

$min\left|\right|d_{\widehat{\mathrm{\θ}}}|{|}_1,s.t.||{\hat{W}}^{-\frac{1}{2}}(\hat{y}-{\widehat{\mathrm{\σ}}}^{2}vec(I_Q\left)\right)-\left({\hat{W}}^{-\frac{1}{2}}{\underset{\‾}{B}}_{\widehat{\mathrm{\θ}}}\right)d_{\widehat{\mathrm{\θ}}}|{|}_2\≤\sqrt{{\mathrm{\η}}_1}---\left(25\right)$

其中是的估计并且通过能够计算出来。根据公式(24)，高概率1-ξ可以用来选择参数η₁，其中ξ是很小的。通常设置ξ＝0.001来估算η₁的值就足够了，通过执行Matlab软件中函数η₁＝chi2inv(1-ξ，Q²)即可。是σ²的估计值，并且通过Q-P个最小奇异值的平均值或者的最小特征值得到。公式(25)能够通过SOC(二阶锥)软件编程包计算，例如Sedumi和CVX。然后通过测绘获得目标的DOA估计。

本发明的方法的主要特点如下：

1、本发明通过降维转换，SNR增益得到提高，并且利用协方差矩阵向量化技术，使得虚拟孔径明显扩大，因此本发明比l₁-SVD和RVl₁-SVD算法具有更好的角度估计性能和更高的角度分辨率；

2、本发明由于将二维完备字典转换成一维的，将多测量矢量(MMV)问题转换成单测量矢量(SMV)问题，并且在l₁范数最小化问题中只涉及实值处理过程，因此本发明计算复杂度明显降低，比传统基于l₁范数的方法具有更低的计算复杂度，且计算时间对目标数具有更好的稳健性；

3、本发明由于实域转换等技术的应用，在低SNR区域角度估计性能优于l₁-SVD和RVl₁-SVD，且具有最低的SNR临界值。

4、表1是不同算法对于不同数目不相关目标的计算时间对比，可以看出图表一通过利用MATLAB软件中TIC和TOC指令，验证了l₁-SVD算法，RVl₁-SVD算法和本发明的计算复杂度。这里假设目标数目P是已知的。所有目标的SNR设置为10dB，快拍数为400。目标的DOA满足θ_i+1-θ_i＝10°(i＝1,2,…,P-1)。对于每个目标数，这些算法的计算时间为通过300次MonteCarlo试验的均值。图表一为这些算法的计算时间对比，可以看出，当目标数增加时l₁-SVD算法和RVl₁-SVD算法的计算时间明显增加，但是本发明的计算时间对目标数具有稳健性。另一方面，本发明的计算时间比l₁-SVD算法和RVl₁-SVD算法都小很多。

表1

附图说明

图1是本发明的整体框架图；

图2不同算法对于不相关目标均方根误差和信噪比的关系图；

图3不同算法对于不相关目标角度估计的均方根误差和快拍数的关系图；

图4不同算法对于不相关目标DOA估计误差和角度间隔的关系图；

图5不同算法对于不相关目标分辨率和信噪比的关系图。

具体实施方式

下面结合附图对本发明做更详细地描述：

结合图一，本发明方法的主要步骤是：1、建立单基地MIMO雷达的接收信号模型，并设计降维矩阵对接收数据进行降维处理；2、利用酉变换矩阵，将降维后的接收数据矩阵变成实域，构造实值扩展数据矩阵并获得其协方差矩阵；3、将实值协方差矩阵向量化以解决MMV问题，并获得稀疏表示框架下的相应模型；4、设计权值矩阵获得估计参数并构造实值l₁范数最小化框架，利用编程软件包SOC(二阶锥)计算方法，获得恢复矩阵，并寻找恢复矩阵中的非零行，实现对单基地MIMO雷达目标DOA的估计。

仿真中将本发明与l₁-SVD算法，实域l₁-SVD(表示为RVl₁-SVD)算法以及克拉美—罗界(CRB)算法进行了对比。在所有仿真中，假设一个单基地MIMO雷达系统，发射阵列和接收阵列的阵元数M＝N＝5，且都是阵元间空间距离为半个波长的均匀线阵。信噪比(SNR)定义为10log(σ_s ²/σ_n ²)，其中σ_s ²和σ_n ²分别表示信号和噪声功率。300次MonteCarlo试验的情况下，DOA估计的均方根误差(RMSE)定义为

$RMSE=\sqrt{\frac{1}{300P}\underset{i=1}{\overset{300}{\Σ}}\underset{p=1}{\overset{P}{\Σ}}{({\widehat{\mathrm{\θ}}}_{i,p}-{\mathrm{\θ}}_p)}^2}---\left(51\right)$

其中表示第i次试验波达角θ_p的估计。仿真中的所有算法，都假设目标数是已知的或者通过MDL准则已被估计，并且对l₁-SVD算法，RVl₁-SVD算法和本发明设置置信区间为0.999(即ξ＝0.001)。所有仿真中，空间网格是均匀的从-90°到90°，间隔为0.01°。

结合图二到图五

1、不同算法对于不相关目标均方根误差和信噪比的关系

图二展示了不同SNR下三种算法的角度估计性能，其中假设不相干目标数P＝3，其DOA分别为θ₁＝-10.13°，θ₂＝-0.13°，θ₃＝20.13°，保持快拍数固定为400，SNR从-10dB到20dB变换。可以看出，在低SNR区域RVl₁-SVD算法比l₁-SVD算法提供更好的角度估计性能并且在高SNR区域角度估计性能相近。与RVl₁-SVD算法和l₁-SVD算法相比，本发明的性能在所有SNR区域优于以上两种算法。原因是在本发明中SNR增益得到增强并且虚拟孔径扩大了。

2、不同算法对于不相关目标角度估计的均方根误差和快拍数的关系

图三展示了不同快拍数下三种算法的角度估计性能，其中假设不相干目标数P＝3，其DOA分别为θ₁＝-10.13°，θ₂＝-0.13°，θ₃＝20.13°，SNR固定为0dB并且快拍数从50到600变换。从图三可以得出结论，本发明比RVl₁-SVD算法和l₁-SVD算法都具有更好的角度估计性能，但是在很低快拍数情况下角度估计性能下降。这是因为在快拍数小的情况下，协方差向量被严重扰乱。因此，本发明在样本量足够大的情况下提供了更好的角度估计性能。

3、不同算法对于不相关目标DOA估计误差和角度间隔的关系

图四展示了这三种算法的DOA估计误差与角度间隔，SNR和快拍数分别固定为5dB和400。考虑两个不相关目标的DOA为θ₁＝0.13°，θ₂＝0.13°+Δθ，其中Δθ从2°到20°变化。从图四可以看出，与RVl₁-SVD算法和l₁-SVD算法相比，本发明对于邻近目标表现出最好的角度估计性能。这是因为在本发明中虚拟孔径显著扩大，即角度分辨率得以提高。

4、不同算法对于不相关目标分辨率和信噪比的关系

图五展示了三种算法的目标分辨率与SNR，其中快拍数固定为400。两个不相关邻近目标的DOA为θ₁＝0.52°，θ₂＝4.52°。在这一仿真中，如果在空间谱中至少有两个尖峰出现且它满足其中和是θ_i的估计，则这两个目标可以认为被成功分辨出来了。可以看出，在高SNR区域所有方法表现出100％正确的角度分辨率。当SNR降低时，各算法的目标分辨率都以某种程度开始下降，即所谓的SNR临界值。本发明与RVl₁-SVD和l₁-SVD算法相比有最低的SNR临界值，即在分辨邻近目标时本发明具有最好的能力。

Claims (5)

1.基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法，其特征在于包括以下步骤：

(1)发射阵列发射相互正交的相位编码信号，接收端进行匹配滤波处理后获得接收数据，建立单基地MIMO雷达的接收信号模型，并利用降维转换矩阵对接收数据进行降维处理；

(2)利用酉变换矩阵，将降维后的接收数据矩阵变成实域，构造实值扩展数据矩阵并获得其协方差矩阵；

(3)根据Khatri-Rao积，将实值协方差矩阵向量化以解决MMV问题，并获得稀疏表示框架下的相应模型；

(4)设计权值矩阵获得估计参数并构造实值l₁范数最小化框架，利用编程软件包二阶锥SOC计算方法，获得恢复矩阵，并寻找恢复矩阵中的非零行，实现对单基地MIMO雷达目标DOA的估计。

2.根据权利要求1所述的基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法，其特征在于，所述的对MIMO雷达系统的接收信号进行降维处理步骤是：

(1)根据单基地MIMO雷达的特殊结构，只有Q＝M+N-1个不同元素，导向矢量可以表示为：

$a_t\left(\mathrm{\θ}\right)\⊗a_r\left(\mathrm{\θ}\right)=G_{MN\×Q}b\left(\mathrm{\θ}\right)---\left(3\right)$

式中a_t(θ)和a_r(θ)分别为发射导向矢量和接收导向矢量，G_MN×Q和b(θ)是降维转换矩阵和导向矢量，它们分别为

$G_{MN\×Q}=\left[\begin{array}{c}{J}_0\\ J_1\\ .\\ .\\ .\\ J_{M-1}\end{array}\right]---\left(4\right)$

$b\left(\mathrm{\θ}\right)=\left\{\begin{array}{c}\[e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}},...,e^{-j\mathrm{\π}\frac{1}{2}\sin \mathrm{\θ}},e^{-j\mathrm{\π}\frac{1}{2}\sin \mathrm{\θ}},...,e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}}\],Q=2k\\ \[e^{-j\mathrm{\π}\frac{Q-1}{2}\sin \mathrm{\θ}},...,1,...,e^{-j\mathrm{\π}\frac{Q-1}{1}\sin \mathrm{\θ}}\],Q=2k+1\end{array}\right.---\left(5\right)$

其中J_m＝[0_N×m,I_N,0_N×(M-m-1)]，m＝0,1,...,M-1

(2)根据公式(4)，定义矩阵F＝G^HG，如下所示

(3)定义降维矩阵W＝F^-(1/2)G^H，满足WW^H＝I_M+N-1。利用W获得降维接收数据Y，有

$Y=F^{(1/2)}BS+WN=\hat{B}S+WN---\left(7\right)$

式中B＝[b(θ₁),...,b(θ_p)]，S和N是信号数据矩阵和噪声矩阵，

3.根据权利要求1所述的基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法，其特征在于，所述的构造实值扩展数据矩阵并获得其协方差矩阵如下步骤：

(1)定义酉变换矩阵为

$U_{2K}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_K& {\mathrm{jI}}_K\\ {\mathrm{\Π}}_K& -{\mathrm{j\Π}}_K\end{array}\right]---\left(9\right)$

其中Π_K表示具有反对角元素为1，其他元素为0的K×K交换矩阵，同时

$U_{2K+1}=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_K& 0& {\mathrm{jI}}_K\\ 0^T& \sqrt{2}& 0^T\\ {\mathrm{\Π}}_K& 0& -{\mathrm{j\Π}}_K\end{array}\right]---\left(10\right)$

利用酉变换矩阵U_Q，新的接收数据矩阵表示为

$\hat{Y}=U_Q^{H}Y=U_Q^H\hat{B}S+U_Q^{H}WN---\left(11\right)$

其中(·)^H表示共轭转置。根据公式(11)，复导向矩阵转换成为实值导向矩阵

(2)数据矩阵可以分为实数部分和虚数部分，即

$\mathrm{Re}\left(\hat{Y}\right)=\stackrel{\‾}{B}\mathrm{Re}\left(S\right)+\mathrm{Re}\left(U_Q^{H}WN\right)---\left(12\right)$

$\mathrm{Im}\left(\hat{Y}\right)=\stackrel{\‾}{B}\mathrm{Im}\left(S\right)+\mathrm{Im}\left(U_Q^{H}WN\right)---\left(13\right)$

其中Re(·)和Im(·)分别表示实数和虚数部分，构造实值扩展数据矩阵，如下所示

$\stackrel{\‾}{Y}=\[\mathrm{Re}\left(\hat{Y}\right),\mathrm{Im}\left(\hat{Y}\right)\]=\stackrel{\‾}{B}\stackrel{\‾}{S}+\stackrel{\‾}{N}---\left(14\right)$

其中 $\stackrel{\‾}{S}=\[\mathrm{Re}\left(S\right),\mathrm{Im}\left(S\right)\]$ 和 $\stackrel{\‾}{N}=\[\mathrm{Re}\left(U_Q^{H}WN\right),\mathrm{Im}\left(U_Q^{H}WN\right)\]$

(3)的实值协方差矩阵表示为

$R_{\stackrel{\‾}{Y}}=E\[\stackrel{\‾}{Y}{\stackrel{\‾}{Y}}^H\]=\stackrel{\‾}{B}{R}_S{\stackrel{\‾}{B}}^H+{\mathrm{\σ}}^2{I}_Q---\left(15\right)$

其中E[·]表示期望，σ²是噪声功率。

4.根据权利要求1所述的基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法，其特征在于，所述的获得稀疏表示框架下的相应模型按如下步骤：

(1)对进行向量化操作，有

其中表示Khatri-Rao积操作，d＝vec(R_S)，数据向量y可看做是具有单快拍虚拟阵列输出的信号模型，于是稀疏信号重构的MMV转变为SMV问题，同时，类似于公式(11)中的实值导向矩阵，

(2)构造DOA估计的一维完备字典为其中且在稀疏表示框架下，公式(16)可以被完备字典代替，如下所示

$y={\underset{\‾}{B}}_{\widehat{\mathrm{\θ}}}{d}_{\widehat{\mathrm{\θ}}}+{\mathrm{\σ}}^{2}vec\left(I_Q\right)---\left(21\right)$

其中是稀疏向量，即其只有P个非零元素，冗杂元素为0，因此通过找到中的非零元素能够对DOA进行估计。

5.根据权利要求1所述的基于稀疏表示的单基地多输入多输出雷达目标波达方向估计方法，其特征在于，所述的获得恢复矩阵，并寻找恢复矩阵中的非零行的步骤是：

(1)构造权值矩阵使得估计误差Δy满足渐近标准正态分布，

$\left|\right|{\stackrel{\‾}{W}}^{-\frac{1}{2}}\mathrm{\Δ}y|{|}_2^2~{\mathrm{As\χ}}^2(0,I_{Q^2})---\left(24\right)$

其中表示具有Q²自由度的渐近卡方分布

(2)构造具有SMV的实值l₁范数最小化问题，

$\min \left|\right|d_{\widehat{\mathrm{\θ}}}|{|}_1,s.t.||{\hat{W}}^{-\frac{1}{2}}(\hat{y}-{\widehat{\mathrm{\σ}}}^{2}vec(I_Q\left)\right)-\left({\hat{W}}^{-\frac{1}{2}}{\underset{\‾}{B}}_{\widehat{\mathrm{\θ}}}\right)d_{\widehat{\mathrm{\θ}}}|{|}_2\≤\sqrt{{\mathrm{\η}}_1}---\left(25\right)$

其中是的估计并且能够通过得到，通过执行Matlab软件中函数η₁＝chi2inv(1-ξ，Q²)，设置ξ＝0.001来估算η₁，是σ²的估计值，通过Q-P个最小奇异值的平均值或者的最小特征值得到，然后测绘获得目标的DOA估计。