CN108732549A - 一种基于协方差矩阵重构的阵元缺损mimo雷达doa估计方法 - Google Patents

Abstract

本发明公开了一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，包括首先接收阵元缺损下的MIMO雷达接收信号经匹配滤波处理，获得虚拟阵列在k时刻输出信号计算阵元缺损MIMO雷达虚拟阵列的协方差矩阵其次以接收阵元数N作为子方块矩阵的大小，对阵元缺损MIMO雷达虚拟阵列的协方差矩阵进行划分；然后利用Toeplitz矩阵的差分性质重构具有Toeplitz特性的子方块矩阵根据各个子方块矩阵的重构值获得重构后的数据协方差矩阵最后根据重构后的协方差矩阵采用ESPRIT算法进行目标DOA估计。本发明能有效解决接收阵元缺损下的MIMO雷达协方差矩阵秩亏的问题，提高MIMO雷达的DOA估计性能。

Description

一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法

技术领域

本发明涉及多输入多输出(MIMO)雷达目标参数估计方法，特别是涉及一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法。

背景技术

多输入多输出(Multiple-Input Multiple-Output，MIMO)雷达的灵感来源于通信系统中的MIMO技术，由于其存在诸多优点，得到众多学者的广泛关注。与传统的相控阵雷达相比，MIMO雷达能灵活地根据需求来设计发射波形，还可通过虚拟阵元来扩展阵列的孔径，以提升目标可识别数目的上限，并提高目标分辨力以及参数估计性能。因此，对MIMO雷达的研究具有重要的意义和实用价值。

波达方向(Direction of Arrival，DOA)估计是雷达信号处理中的一项重要研究内容，针对MIMO雷达的DOA估计，人们先后提出了多重信号分类(multiple signalclassification，MUSIC)算法和旋转不变子空间(estimation of signal parameters viarotational invariance technique，ESPRIT)算法等子空间类算法以及其他算法。在众多DOA估计算法中，子空间类算法因计算方便、估计性能优越而获得广泛的应用。子空间类算法的思想是对协方差矩阵进行特征分解求取特征向量，分离出信号子空间和噪声子空间，然后利用信号子空间与噪声子空间的正交性或信号子空间本身的旋转不变性进行目标DOA估计。随着MIMO雷达天线阵列规模的不断增加，以及受恶劣环境或人为干扰等因素的影响，大大增加了天线阵元损坏的概率。阵元缺损导致部分阵元的目标接收信号丢失，使协方差矩阵的秩发生退化，从而导致子空间类算法性能下降甚至完全失效。针对传统线性阵列阵元缺损下的DOA估计问题，现有技术一利用天线阵元间的相关特性，将单次快拍下的信号矢量变换为等效低秩矩阵，通过矩阵填充算法估计缺失数据，但是该方法需逐个快拍填充缺失数据，具有较大的计算量。针对MIMO雷达阵元缺损下的DOA估计问题，现有技术二提出了一种差分阵处理方法，该方法根据KR积变换理论，建立MIMO雷达虚拟阵元差分阵列的接收数据模型，通过差分阵列中正常工作阵元的接收数据恢复相应的缺失数据，最后采用空间平滑算法进行解相干来提高阵元缺损MIMO雷达的DOA估计性能，但该方法需通过空间平滑算法来获得满秩协方差矩阵，当MIMO雷达的天线阵元较多时，计算量过大。在对运动目标的单快拍成像时，人们通过利用稀疏阵列MIMO雷达采集数据来降低硬件成本，其被稀疏掉的天线阵元可等价于相应位置上的阵元损坏，现有技术三采用矩阵填充算法恢复被稀疏掉阵元的接收数据，当面阵中出现整行或整列接收天线阵元同时被稀疏掉时，稀疏MIMO雷达的接收数据无法满足矩阵填充所要求的的强不相干性条件，因此该文通过Hankel矩阵变换将其转换成二重块Hankel矩阵，再采用矩阵填充算法对二重块Hankel矩阵进行填补，从而重构出被稀疏掉阵元的接收数据，以提高MIMO雷达的成像性能。现有技术四针对接收阵元缺损下的MIMO声呐/雷达的DOA估计问题，提出了一种基于差分阵列的数据协方差矩阵重构的方法，该方法利用MIMO声呐/雷达能形成虚拟均匀线性满阵时其协方差矩阵满足Toeplitz的性质，通过正常的协方差矩阵元素来恢复相应的缺失元素，然而该方法仅适用于发射和接收阵元间距满足特定关系时阵元缺损MIMO声呐/雷达的DOA估计问题。在大多数情况下，MIMO雷达虚拟阵列的数据协方差矩阵为非Toeplitz矩阵，此时利用Toeplitz矩阵的差分性质对协方差矩阵进行重构会破坏原协方差矩阵的结构，从而导致DOA估计性能下降。因此，研究如何恢复MIMO雷达缺损阵元的数据以消除阵元缺损对目标DOA估计带来的不利影响是非常有必要的。

发明内容

发明目的：为有效解决接收阵元缺损下的MIMO雷达协方差矩阵秩亏的问题，从而提高MIMO雷达的DOA估计性能，本发明提供一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法。

技术方案：本发明提供了一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，包括以下步骤：

(1)接收阵元缺损下的MIMO雷达接收信号经匹配滤波处理，从而获得虚拟阵列在k时刻输出信号计算阵元缺损MIMO雷达虚拟阵列的协方差矩阵其中，K为快拍数，(·)^H表示共轭转置运算；

(2)以接收阵元数N作为子方块矩阵的大小，对阵元缺损MIMO雷达虚拟阵列的协方差矩阵进行划分，即其中，M为发射阵元数，为协方差矩阵中第(i,j)位置处的N×N子方块矩阵；

(3)为了恢复各子方块矩阵中的缺失元素，利用Toeplitz矩阵的差分性质重构具有Toeplitz特性的子方块矩阵根据各个子方块矩阵的重构值获得重构后的数据协方差矩阵

(4)根据重构后的协方差矩阵采用ESPRIT算法进行目标DOA估计。

进一步的，所述步骤(1)中假设一个阵元正常下的MIMO雷达系统具有M个发射阵元，N个接收阵元，发射与接收阵列均为均匀线阵，阵元间距分别为d_t和d_r，假设有L个远场窄带非相干目标，第l(l＝1,...,L)个目标的方位角度为θ_l，则k时刻阵列接收信号经匹配滤波后的接收数据向量为：

x(k)＝As(k)+n(k) (1)；

其中，为阵元正常下的MIMO雷达的虚拟阵列流型矩阵，表示Kronecker积；为发射阵列的导向矢量，为接收阵列的导向矢量，其中，l＝1,...,L，(·)^T表示转置，λ为载波波长；发射信号为s(k)＝[s₁(k),…,s_L(k)]^T；n(k)为高斯白噪声且与信号彼此独立；

根据上述假设，阵元正常下的MIMO雷达虚拟阵列的协方差矩阵为：

其中，E(·)表示取数学期望；(·)^H表示共轭转置；为信源协方差矩阵，其中，为第l个信源的功率；为噪声功率；I为单位矩阵；

若阵元缺损MIMO雷达中第q(1≤q≤N)个接收阵元缺损，则阵元缺损MIMO雷达虚拟阵列流型中exp{-j2π[(q-1)d_r+md_t]sin(θ_l)/λ},m＝0,…,M-1这些项均被0取代，从而得到阵元缺损MIMO雷达的虚拟阵列流型矩阵则阵元缺损MIMO雷达虚拟阵列的协方差矩阵为：

为了分析阵元正常下的MIMO雷达协方差矩阵R的结构特性，将协方差矩阵R中的任意元素可表示为：

其中，m,m'＝1,2,…,M，n,n'＝1,2,…,N；[(m-1)N+n,(m'-1)N+n']表示该元素在矩阵R中的位置；

假设c＝(m'-1)N+n'-[(m-1)N+n]为元素位置的差分值，则：

定义C＝{(m'-1)N+n'-[(m-1)N+n]}为协方差矩阵中的元素位置的差分集合，定义R(c)＝{R[(m-1)N+n,(m'-1)N+n']|(m'-1)N+n'-[(m-1)N+n]＝c∈C}为差分值等于c(c∈C)的协方差矩阵元素；假设阵元正常下的MIMO雷达收发阵列阵元间距之间的关系为d_t＝Wd_r,由式(4)可得：

通过式(6)可以看出，仅当d_t＝Nd_r时，属于同一差分值的阵元正常下的MIMO雷达虚拟阵列的协方差矩阵元素固定不变，此时该协方差矩阵为Toeplitz矩阵；当d_t≠Nd_r时，阵元正常下的MIMO雷达虚拟阵列的协方差矩阵不再是Toeplitz矩阵；

实际接收数据矩阵长度是有限长的，则阵元缺损MIMO雷达的数据协方差矩阵的近似估计为其中，为阵元缺损MIMO雷达虚拟阵列在k时刻的输出信号；K为快拍数。

进一步的，所述步骤(2)中阵元正常下的MIMO雷达虚拟阵列流型矩阵表示为：

其中，A_t＝[a_t(θ₁),…,a_t(θ_L)]为发射阵列流型矩阵；A_r＝[a_r(θ₁),…,a_r(θ_L)]为接收阵列流型矩阵；D_m(·)为取矩阵的第m行为对角元素构造对角矩阵；

由阵元正常下的MIMO雷达虚拟阵列的协方差矩阵的表达式可得MIMO雷达虚拟阵列的协方差矩阵表示为：

以接收阵元数划分阵元正常下的MIMO雷达虚拟阵列的协方差矩阵R，从而获得矩阵行数与列数均等于接收阵元数的子方块矩阵；以划分后的子方块矩阵来表示协方差矩阵R，即：

结合式(8)和式(9)，第(i,j)位置处的子方块矩阵表示为：

其中，

由式(10)看出，子方块矩阵r_ij可等价于阵元数为N的均匀线阵的接收数据的协方差矩阵，其中，A_r可认为是均匀线阵的阵列流型矩阵，可认为是信源协方差矩阵；

按照式(9)划分阵元缺损MIMO雷达虚拟阵列的协方差矩阵获得矩阵行数与列数均等于接收阵元数的子方块矩阵以划分后的子方块矩阵来表示协方差矩阵即：

进一步的，所述步骤(3)中定义S＝{u-v|u,v＝1,…,N}为子方块矩阵中元素位置的差分集合，定义P(s)＝{(u,v)|u-v＝s∈S,且为差分集合中值等于s(s∈S)且的(u,v)对的个数，则得到子方块矩阵中差分值等于s的非零元素的均值为：

利用Toeplitz矩阵的差分性质重构子方块矩阵得：

依次对协方差矩阵中的各子方块矩阵进行重构，获得重构后的数据协方差矩阵为：

进一步的，所述步骤(4)中根据重构后的协方差矩阵采用ESPRIT算法进行DOA估计，对进行特征分解得到：

其中，Λ_s为L个大特征值构成的对角矩阵；U_s为对应于L个大特征值的特征向量；Λ_n为MN-L个小特征值构成的对角矩阵；U_n为对应于MN-L个小特征值的特征向量；

根据ESPRIT算法可知，信号子空间与阵列流型间的关系为U_s＝AT，其中T为唯一的非奇异矩阵；分别取出阵列流型矩阵A的前N(M-1)行和后N(M-1)行为A₁和A₂，由子阵阵列流型间的旋转不变性可知A₂＝A₁Φ_t，再取出信号子空间U_s的前N(M-1)行和后N(M-1)行为U₁和U₂；

由信号子空间与阵列流型之间的关系可知U₁＝A₁T，U₂＝A₂T，则矩阵U₁和U₂的关系为U₂＝U₁Ψ_t，其中Ψ_t＝T^-1Φ_tT；通过矩阵U₁和U₂可求得然后对Ψ_t进行特征值分解得到特征值χ_l(l＝1,…,L)，则第l个目标的DOA估计值为：

其中，angle(·)表示取相角。

有益效果：与现有技术相比，本发明具有以下优点：

(1)在实际应用中，受恶劣环境和人为干扰等因素的影响，MIMO雷达天线阵元损坏的概率大大增加，阵元损坏会导致协方差矩阵的秩退化，从而导致子空间类算法性能下降甚至完全失效。本发明根据MIMO雷达的协方差矩阵是由子方块Toeplitz矩阵组成的特性，利用Toeplitz矩阵的差分性质，从而重构出满秩的协方差矩阵，达到了恢复协方差矩阵中缺失数据的目的，有效地解决了接收阵元缺损下的MIMO雷达DOA估计问题。

(2)本发明方法利用子方块矩阵中属于同一差分值的非零元素均值来重构协方差矩阵，其计算量较低、实时性较高，同时起到协方差矩阵的降噪效果，从而有效提高接收阵元缺损MIMO雷达的DOA估计性能。

(3)本发明方法不要求MIMO雷达的发射和接收阵元间距满足某种特定的关系，能有效适用于任意阵列结构下的接收阵元缺损MIMO雷达的DOA估计问题，具有普遍适用性和通用性。

附图说明

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

图2是收发阵元间距满足d_t＝Nd_r时的目标DOA估计结果图；

图3是收发阵元间距不满足d_t＝Nd_r时的目标DOA估计结果图；

图4是DOA估计均方根误差随信噪比的变化关系；

图5是DOA估计均方根误差随快拍数的变化关系；

图6是DOA估计均方根误差随缺损接收阵元数的变化曲线。

具体实施方式

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

在实际应用中由于恶劣环境或人为干扰等因素而导致MIMO雷达部分阵元缺损，使得其接收数据缺失及其协方差矩阵秩亏，从而导致子空间类算法的DOA估计性能恶化甚至完全失效。本发明提出了一种接收阵元缺损下基于协方差矩阵重构的MIMO雷达DOA估计方法，该方法根据MIMO雷达协方差矩阵中以接收阵元数划分的子方块矩阵具有Toeplitz特性，利用正常工作接收阵元的协方差矩阵元素来恢复相应的缺失元素，从而重构出完整的数据协方差矩阵，提高阵元缺损MIMO雷达的DOA估计性能。本发明方法能有效解决接收阵元缺损下的MIMO雷达协方差矩阵秩亏的问题，从而提高MIMO雷达的DOA估计性能。

如图1所示，本发明的一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，包括以下步骤：

(1)接收阵元缺损下的MIMO雷达接收信号经匹配滤波处理，从而获得虚拟阵列在k时刻输出信号计算阵元缺损MIMO雷达虚拟阵列的协方差矩阵其中，K为快拍数，(·)^H表示共轭转置运算。

假设一个阵元正常下的MIMO雷达系统具有M个发射阵元，N个接收阵元，发射与接收阵列均为均匀线阵，阵元间距分别为d_t和d_r。假设有L个远场窄带非相干目标，第l(l＝1,...,L)个目标的方位角度为θ_l，则k时刻阵列接收信号经匹配滤波后的接收数据向量为：

x(k)＝As(k)+n(k) (1)；

其中，为阵元正常下的MIMO雷达的虚拟阵列流型矩阵，表示Kronecker积；为发射阵列的导向矢量，为接收阵列的导向矢量，其中，l＝1,...,L，(·)^T表示转置，λ为载波波长；发射信号为s(k)＝[s₁(k),…,s_L(k)]^T；n(k)为高斯白噪声且与信号彼此独立。

根据上述假设，阵元正常下的MIMO雷达虚拟阵列的协方差矩阵为：

其中，E(·)表示取数学期望；(·)^H表示共轭转置；为信源协方差矩阵，其中，为第l个信源的功率；为噪声功率；I为单位矩阵。

若阵元缺损MIMO雷达中第q(1≤q≤N)个接收阵元缺损，则阵元缺损MIMO雷达虚拟阵列流型中exp{-j2π[(q-1)d_r+md_t]sin(θ_l)/λ}，m＝0,…,M-1这些项均被0取代，从而得到阵元缺损MIMO雷达的虚拟阵列流型矩阵则阵元缺损MIMO雷达虚拟阵列的协方差矩阵为：

为了分析阵元正常下的MIMO雷达协方差矩阵R的结构特性，将协方差矩阵R中的任意元素可表示为：

其中，m,m'＝1,2,…,M，n,n'＝1,2,…,N；[(m-1)N+n,(m'-1)N+n']表示该元素在矩阵R中的位置；

假设c＝(m'-1)N+n'-[(m-1)N+n]为元素位置的差分值，则：

定义C＝{(m'-1)N+n'-[(m-1)N+n]}为协方差矩阵中的元素位置的差分集合，定义R(c)＝{R[(m-1)N+n,(m'-1)N+n']|(m'-1)N+n'-[(m-1)N+n]＝c∈C}为差分值等于c(c∈C)的协方差矩阵元素。假设阵元正常下的MIMO雷达收发阵列阵元间距之间的关系为d_t＝Wd_r,由式(4)可得：

通过式(6)可以看出，仅当d_t＝Nd_r时，属于同一差分值的阵元正常下的MIMO雷达虚拟阵列的协方差矩阵元素固定不变，此时该协方差矩阵为Toeplitz矩阵。当d_t≠Nd_r时，阵元正常下的MIMO雷达虚拟阵列的协方差矩阵不再是Toeplitz矩阵，此时采用现有技术四的方法重构协方差矩阵时，会将本来不具有Toeplitz特性的协方差矩阵重构为具有Toeplitz特性的协方差矩阵，这样会破坏MIMO雷达协方差矩阵的结构，从而导致DOA估计的性能不理想。

实际接收数据矩阵长度是有限长的，则阵元缺损MIMO雷达虚拟阵列的协方差矩阵的近似估计为其中，为阵元缺损MIMO雷达虚拟阵列在k时刻的输出信号；K为快拍数。

步骤2：以接收阵元数N作为子方块矩阵的大小，对阵元缺损MIMO雷达虚拟阵列的协方差矩阵进行划分，即其中，M为发射阵元数，为协方差矩阵中第(i,j)位置处的N×N子方块矩阵。

为了能重构阵元缺损下MIMO雷达的协方差矩阵，本发明分析了在阵元正常下的MIMO雷达协方差矩阵中以接收阵元数划分的子方块矩阵的结构特性。阵元正常下的MIMO雷达虚拟阵列流型矩阵表示为：

其中，A_t＝[a_t(θ₁),…,a_t(θ_L)]为发射阵列流型矩阵；A_r＝[a_r(θ₁),…,a_r(θ_L)]为接收阵列流型矩阵；D_m(·)为取矩阵的第m行为对角元素构造对角矩阵。

由式(2)可得阵元正常下的MIMO雷达虚拟阵列的协方差矩阵表示为：

以接收阵元数划分阵元正常下的MIMO雷达虚拟阵列的协方差矩阵R，从而获得矩阵行数与列数均等于接收阵元数的子方块矩阵。以划分后的子方块矩阵来表示协方差矩阵R，即：

结合式(8)和式(9)，第(i,j)位置处的子方块矩阵表示为：

其中，

由式(10)看出，子方块矩阵r_ij可等价于阵元数为N的均匀线阵的接收数据的协方差矩阵，其中，A_r可认为是均匀线阵的阵列流型矩阵，可认为是信源协方差矩阵。均匀线阵的接收数据的协方差矩阵为Toeplitz矩阵。因此阵元正常下的MIMO雷达虚拟阵列协方差矩阵R中的各子方块矩阵r_ij具有Toeplitz特性，从而可以根据Toeplitz矩阵的差分性质，在各子方块矩阵中利用正常矩阵元素来恢复相应的缺失元素，重构出完整的阵元缺损MIMO雷达虚拟阵列的协方差矩阵。

按照式(9)划分阵元缺损MIMO雷达虚拟阵列的协方差矩阵获得矩阵行数与列数均等于接收阵元数的子方块矩阵以划分后的子方块矩阵来表示协方差矩阵即：

步骤3：为了恢复各子方块矩阵中的缺失元素，利用Toeplitz矩阵的差分性质重构具有Toeplitz特性的子方块矩阵根据各个子方块矩阵的重构值获得重构后的数据协方差矩阵

定义S＝{u-v|u,v＝1,…,N}为子方块矩阵中元素位置的差分集合，定义P(s)＝{(u,v)|u-v＝s∈S,且为差分集合中值等于s，s为差分集合S中一个元素，即s∈S，且的(u,v)对的个数，则可以得到子方块矩阵中差分值等于s的非零元素的均值为：

利用Toeplitz矩阵的差分性质重构子方块矩阵得：

依次对协方差矩阵中的各子方块矩阵进行重构，获得重构后的数据协方差矩阵为：

步骤4：根据重构后的协方差矩阵采用ESPRIT算法进行目标DOA估计。

根据重构后的协方差矩阵采用ESPRIT算法进行DOA估计，对进行特征分解得到：

其中，Λ_s为L个大特征值构成的对角矩阵；U_s为对应于L个大特征值的特征向量；Λ_n为MN-L个小特征值构成的对角矩阵；U_n为对应于MN-L个小特征值的特征向量。

根据ESPRIT算法的原理可知，信号子空间与阵列流型间的关系为U_s＝AT，其中T为唯一的非奇异矩阵。分别取出阵列流型矩阵A的前N(M-1)行和后N(M-1)行为A₁和A₂，由子阵阵列流型间的旋转不变性可知A₂＝A₁Φ_t，再取出信号子空间U_s的前N(M-1)行和后N(M-1)行为U₁和U₂。

由信号子空间与阵列流型之间的关系可知U₁＝A₁T，U₂＝A₂T，则矩阵U₁和U₂的关系为U₂＝U₁Ψ_t，其中Ψ_t＝T^-1Φ_tT。通过矩阵U₁和U₂可求得然后对Ψ_t进行特征值分解得到特征值χ_l(l＝1,…,L)，则第l个目标的DOA估计值为：

其中，angle(·)表示取相角。

本发明的技术效果可以通过以下仿真结果进一步说明。为了验证本发明方法在接收阵元缺损MIMO雷达DOA估计方面的有效性，设计了以下几组仿真实验，将现有技术四的方法和本发明方法的阵元缺损MIMO雷达DOA估计性能进行对比，并以阵元正常工作时和阵元缺损时MIMO雷达直接利用ESPRIT算法来估计目标DOA的性能作为参照。为了比较一致，在现有技术四的方法中，应用ESPRIT算法从重构后协方差矩阵中估计目标DOA。仿真中，MIMO雷达的发射阵元数为5，接收阵元数为25，回波噪声选取零均值的加性高斯白噪声。

定义DOA估计的均方根误差(Root Mean Square Error,RMSE)为：

式中，M_T为蒙特卡洛实验次数，为第m_t次蒙特卡洛实验中第l个目标的DOA估计值。

仿真实验1：假设接收阵列中存在8个随机的天线阵元缺损，接收阵元间距d_r＝0.5λ，发射阵元间距d_t＝Nd_r＝12.5λ，选取方位角度为10°、13°和20°的三个目标，快拍数K＝300，信噪比(Signal-Noise Ratio，SNR)为-5dB，进行50次蒙特卡洛实验，仿真结果如图2所示。由图2可知，阵元缺损下直接采用ESPRIT算法估计的目标DOA与真实值相比存在较大波动。收发阵元间距满足d_t＝Nd_r时，MIMO雷达虚拟阵列的协方差矩阵为Toeplitz矩阵，因此本发明方法和现有技术四的方法估计的目标DOA基本与真实值重合。

仿真实验2：假设接收阵列中存在8个随机的天线阵元缺损，接收阵元间距d_r＝0.5λ，发射阵元间距d_t＝d_r＝0.5λ，选取方位角度为10°、13°和20°的三个目标，快拍数K＝300，信噪比为-5dB，进行50次蒙特卡洛实验，仿真结果如图3所示。由图3可知，阵元缺损下直接采用ESPRIT算法估计的目标DOA存在较大的误差。收发阵元间距不满足d_t＝Nd_r时，MIMO雷达虚拟阵列的协方差矩阵不再是Toeplitz矩阵，现有技术四的方法将本来不具有Toeplitz特性的协方差矩阵重构为具有Toeplitz特性的协方差矩阵，破坏了原协方差矩阵的结构，因此现有技术四的方法的目标DOA估计误差明显高于本发明方法。

仿真实验3：假设接收阵列中存在8个随机的天线阵元缺损，选取方位角度为-10°、0°和15°的三个目标，快拍数K＝300，信噪比变化范围为-20dB～0dB，进行300次蒙特卡洛实验，仿真结果如图4所示。由图4可知，随着信噪比的不断增加，本发明方法的DOA估计精度逐渐提高，其精度明显优于现有技术四的方法的估计精度。此外，本发明方法利用子方块矩阵中非零元素的均值来恢复其缺失元素，具有降噪的效果，因此其DOA估计精度也优于阵元正常工作下的DOA估计精度。

仿真实验4：假设接收阵列中存在8个随机的天线阵元缺损，选取方位角度为-10°、0°和15°的三个目标，信噪比为-5dB，令快拍数由50～350变化，进行300次蒙特卡洛实验，仿真结果如图5所示。由图5可知，随着快拍数的不断增多，估计的协方差矩阵更加趋近于理想协方差矩阵，因此不同方法的DOA估计精度随着快拍数增加而逐渐提高。在不同的快拍数下，本发明方法始终保持最佳的DOA估计精度。

仿真实验5：为了检验本发明方法在不同缺损接收阵元数下的稳健性，选取方位角度为-10°、0°和15°的三个目标，快拍数K＝300，信噪比为-5dB，缺损接收阵元数由0～16依次增加，进行300次蒙特卡洛实验，仿真结果如图6所示。由图6可知，在缺损接收阵元数不断增加的过程中，本发明方法的DOA估计性能始终最优，具有良好的稳健性。

Claims (5)

1.一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，其特征在于，包括以下步骤：

(1)接收阵元缺损下的MIMO雷达接收信号经匹配滤波处理，从而获得虚拟阵列在k时刻输出信号计算阵元缺损MIMO雷达虚拟阵列的协方差矩阵其中，K为快拍数，(·)^H表示共轭转置运算；

(2)以接收阵元数N作为子方块矩阵的大小，对阵元缺损MIMO雷达虚拟阵列的协方差矩阵进行划分，即其中，M为发射阵元数，为协方差矩阵中第(i,j)位置处的N×N子方块矩阵；

(3)为了恢复各子方块矩阵中的缺失元素，利用Toeplitz矩阵的差分性质重构具有Toeplitz特性的子方块矩阵根据各个子方块矩阵的重构值获得重构后的数据协方差矩阵

(4)根据重构后的协方差矩阵采用ESPRIT算法进行目标DOA估计。

2.根据权利要求1所述的一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，其特征在于：所述步骤(1)中假设一个阵元正常下的MIMO雷达系统具有M个发射阵元，N个接收阵元，发射与接收阵列均为均匀线阵，阵元间距分别为d_t和d_r，假设有L个远场窄带非相干目标，第l(l＝1,...,L)个目标的方位角度为θ_l，则k时刻阵列接收信号经匹配滤波后的接收数据向量为：

x(k)＝As(k)+n(k) (1)；

其中，为阵元正常下的MIMO雷达的虚拟阵列流型矩阵，表示Kronecker积；为发射阵列的导向矢量，为接收阵列的导向矢量，其中，l＝1,...,L，(·)^T表示转置，λ为载波波长；发射信号为s(k)＝[s₁(k),…,s_L(k)]^T；n(k)为高斯白噪声且与信号彼此独立；

根据上述假设，阵元正常下的MIMO雷达虚拟阵列的协方差矩阵为：

其中，E(·)表示取数学期望；(·)^H表示共轭转置；为信源协方差矩阵，其中，为第l个信源的功率；为噪声功率；I为单位矩阵；

若阵元缺损MIMO雷达中第q(1≤q≤N)个接收阵元缺损，则阵元缺损MIMO雷达虚拟阵列流型中exp{-j2π[(q-1)d_r+md_t]sin(θ_l)/λ},m＝0,…,M-1这些项均被0取代，从而得到阵元缺损MIMO雷达的虚拟阵列流型矩阵则阵元缺损MIMO雷达虚拟阵列的协方差矩阵为：

为了分析阵元正常下的MIMO雷达协方差矩阵R的结构特性，将协方差矩阵R中的任意元素可表示为：

其中，m,m'＝1,2,…,M，n,n'＝1,2,…,N；[(m-1)N+n,(m'-1)N+n']表示该元素在矩阵R中的位置；

假设c＝(m'-1)N+n'-[(m-1)N+n]为元素位置的差分值，则：

定义C＝{(m'-1)N+n'-[(m-1)N+n]}为协方差矩阵中的元素位置的差分集合，定义R(c)＝{R[(m-1)N+n,(m'-1)N+n']|(m'-1)N+n'-[(m-1)N+n]＝c∈C}为差分值等于c(c∈C)的协方差矩阵元素；假设阵元正常下的MIMO雷达收发阵列阵元间距之间的关系为由式(4)可得：

通过式(6)可以看出，仅当d_t＝Nd_r时，属于同一差分值的阵元正常下的MIMO雷达虚拟阵列的协方差矩阵元素固定不变，此时该协方差矩阵为Toeplitz矩阵；当d_t≠Nd_r时，阵元正常下的MIMO雷达虚拟阵列的协方差矩阵不再是Toeplitz矩阵；

实际接收数据矩阵长度是有限长的，则阵元缺损MIMO雷达的数据协方差矩阵的近似估计为其中，为阵元缺损MIMO雷达虚拟阵列在k时刻的输出信号；K为快拍数。

3.根据权利要求1所述的一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，其特征在于：所述步骤(2)中阵元正常下的MIMO雷达虚拟阵列流型矩阵表示为：

其中，A_t＝[a_t(θ₁),…,a_t(θ_L)]为发射阵列流型矩阵；A_r＝[a_r(θ₁),…,a_r(θ_L)]为接收阵列流型矩阵；D_m(·)为取矩阵的第m行为对角元素构造对角矩阵；

由阵元正常下的MIMO雷达虚拟阵列的协方差矩阵的表达式可得MIMO雷达虚拟阵列的协方差矩阵表示为：

以接收阵元数划分阵元正常下的MIMO雷达虚拟阵列的协方差矩阵R，从而获得矩阵行数与列数均等于接收阵元数的子方块矩阵；以划分后的子方块矩阵来表示协方差矩阵R，即：

结合式(8)和式(9)，第(i,j)位置处的子方块矩阵表示为：

其中，

由式(10)看出，子方块矩阵r_ij可等价于阵元数为N的均匀线阵的接收数据的协方差矩阵，其中，A_r可认为是均匀线阵的阵列流型矩阵，可认为是信源协方差矩阵；

按照式(9)划分阵元缺损MIMO雷达虚拟阵列的协方差矩阵获得矩阵行数与列数均等于接收阵元数的子方块矩阵以划分后的子方块矩阵来表示协方差矩阵即：

4.根据权利要求1所述的一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，其特征在于：所述步骤(3)中定义S＝{u-v|u,v＝1,…,N}为子方块矩阵中元素位置的差分集合，定义P(s)＝{(u,v)|u-v＝s∈S,且为差分集合中值等于s(s∈S)且的(u,v)对的个数，则得到子方块矩阵中差分值等于s的非零元素的均值为：

利用Toeplitz矩阵的差分性质重构子方块矩阵得：

依次对协方差矩阵中的各子方块矩阵进行重构，获得重构后的数据协方差矩阵为：

5.根据权利要求1所述的一种基于协方差矩阵重构的阵元缺损MIMO雷达DOA估计方法，其特征在于：所述步骤(4)中根据重构后的协方差矩阵采用ESPRIT算法进行DOA估计，对进行特征分解得到：

其中，Λ_s为L个大特征值构成的对角矩阵；U_s为对应于L个大特征值的特征向量；Λ_n为MN-L个小特征值构成的对角矩阵；U_n为对应于MN-L个小特征值的特征向量；

根据ESPRIT算法可知，信号子空间与阵列流型间的关系为U_s＝AT，其中T为唯一的非奇异矩阵；分别取出阵列流型矩阵A的前N(M-1)行和后N(M-1)行为A₁和A₂，由子阵阵列流型间的旋转不变性可知A₂＝A₁Φ_t，再取出信号子空间U_s的前N(M-1)行和后N(M-1)行为U₁和U₂；

由信号子空间与阵列流型之间的关系可知U₁＝A₁T，U₂＝A₂T，则矩阵U₁和U₂的关系为U₂＝U₁Ψ_t，其中Ψ_t＝T^-1Φ_tT；通过矩阵U₁和U₂可求得然后对Ψ_t进行特征值分解得到特征值χ_l(l＝1,…,L)，则第l个目标的DOA估计值为：

其中，angle(·)表示取相角。