CN103760547B - 基于互相关矩阵的双基mimo雷达角度估计方法 - Google Patents

Abstract

本发明公开了一种基于互相关矩阵的双基MIMO雷达角度估计方法，主要解决双基MIMO雷达角度估计运算量大、计算复杂的问题。其实现步骤是：1）对雷达回波信号进行匹配滤波，并按发射阵和接收阵形成数据；2）用形成数据分别构造互协方差矩阵；3）用互协方差矩阵行向量的线性无关性分别求解发射阵和接收阵导向矢量零空间的正交投影算子，并构造发射阵和接收阵求根多项式；4）求解目标相对于发射阵和接收阵位置；5）用匹配滤波后数据自相关协方差矩阵行向量的线性无关性求解合成阵列导向矢量零空间的正交投影算子，并构造代价函数进行角度配对。本发明以小的运算量实现了高精度的MIMO雷达目标角度估计，可用于雷达、通信中的目标定位。

Description

基于互相关矩阵的双基MIMO雷达角度估计方法

技术领域

本发明属于雷达技术领域，具体涉及到双基地MIMO雷达角度估计方法，可用于雷达目标参数估计与检测。

背景技术

双基地MIMO雷达目标相对于发射阵位置DOD和接收阵位置DOA估计引起了很多学者的研究兴趣。大部分的研究方法基于特征分解类算法。如K.T.Wong等人在IEEETransactionOnAntennaandPropagation期刊第48卷第8期第1235页至1245页提出了二维多重信号分类算法MUSIC。X.Zhang等人在IEEECommunicationsLetters期刊第14卷第12期第1161页至1163页提出了降维MUSIC算法。这两种算法性能相近。陈多芳等人2008年在ElectronicsLetters期刊第44卷第12期第770页至771页提出的双基地MIMO雷达ESPRIT算法，该算法利用发射阵和接收阵的旋转不变性估计双基地MIMO雷达到达角，并获得了较好的估计性能。M.L.Bencheikh等人2010年在ElectronicsLetters期刊第46卷第15期第1081页至1083页提出了联合MUSIC-ESPRIT算法，相比多项式求根算法，该算法运算量有一定的减少，但仍然无法满足工程需要。

特征分解类算法的复杂度受阵列孔径影响较大，为了降低运算量及加性噪声对算法的影响，JingminXin等人2004年在IEEETransactiononSignalProcessing期刊第52卷第4期第876页至893页提出了无特征值分解的子空间算法来估计一维DOA。尽管该算法估计性能相比特征分解类算法较差，但是运算量小，并且在低信噪比和小快拍条件下，能够较好的估计两个离得很近的相干和非相干目标。GuangminWang等人2011年在IEEETransactiononSignalProcessing期刊第59卷第7期第3197页至3212页提出了无特征值分解的子空间算法来估计二维DOA。该算法估计性能优于CCM-ESPRIT和JSVD算法，但是由于限定方位维和俯仰维阵元数相等，大大限制了阵列排布；此外，该算法涉及的L型阵列接收的方位维和俯仰维数据相互独立，而双基地MIMO雷达接收的DOD和DOA数据包含在同一数据中，因此，该算法无法直接应用于双基地MIMO雷达角度估计。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于互相关矩阵的双基MIMO雷达角度估计方法，以避免大的运算量，减弱加性噪声的影响，提高角度估计精度，实现对双基地MIMO雷达角度的快速估计。

实现本发明目的技术思路是：从接收数据分离DOD和DOA数据并构造新的互相关协方差矩阵，根据互相关矩阵行向量的线性无关性求解导向矢量零空间的投影算子，利用求根MUSIC算法求解目标角度，以合成阵列导向矢量零空间的投影算子构造代价函数实现DOD和DOA的配对。具体实现步骤包含如下：

1)M个发射天线发射电磁波照射空间P个不相关目标，N个接收天线接收目标散射回波，并对回波信号依次进行混频、低通滤波、通道分离，得到MN×1维发射阵排列数据X＝As(t)+n(t)其中A＝[b₁,b₂,…，b_i,…,b_P]，为合成阵列导向矢量，i＝1,…,P， 为Kronecker积；a_r(φ_i)和a_t(θ_i)分别为第i个目标的接收阵和发射阵导向矢量，θ_i为第i个目标相对于发射阵位置，φ_i为第i个目标相对于接收阵位置；s(t)为发射信号矢量，n(t)是复高斯白噪声矢量；M、N为大于零的整数；

2)将MN×1维发射阵排列数据X按照接收阵排列，得到NM×1维接收阵排列数据Y；

3)构造发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy和接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx；

3a)用MN×1维发射阵排列数据X中的M个元素和NM×1维接收阵排列数据Y中的N个元素构成第一互相关矩阵R_xykj,k＝1,…,N,j＝1,…,M；用MN×1维发射阵排列数据X中的M个元素共轭反转和NM×1维接收阵排列数据Y中的N个元素构成第二互相关矩阵将这两个互相关矩阵按照行排列，构成发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy的一个子阵，该互相关协方差矩阵R_xy的子阵数小于等于M和N中的最小值；

3b)用NM×1维接收阵排列数据Y中的N个元素和MN×1维发射阵排列数据X中的M个元素构成第三互相关矩阵R_yxjk；用NM×1维接收阵排列数据Y中的N个元素共轭反转和MN×1维发射阵排列数据X中的M个元素构成互第四互相关矩阵将这两个互相关矩阵按照行排列，构成接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx的一个子阵，该互相关协方差矩阵R_yx的子阵数小于等于M和N中的最小值；

4)将所述发射阵排列数据X与接收阵排列数据Y的互相关协方差矩阵R_xy划分为P×2NΓ维子阵和(M-P)×2NΓ维子阵将所述接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx划分为P×2MΓ维子阵和(N-P)×2MΓ维子阵

5）利用所述R_xy行向量的线性无关性和所述R_yx行向量的线性无关性，计算发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ：

6)利用发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ，分别构造发射阵求根多项式f(z_t)和接收阵求根多项式f(z_r)；

7）分别求解发射阵求根多项式单位圆上的根z_t和接收阵求根多项式单位圆上的根z_r，得到目标相对于发射阵位置估计值和目标相对于接收阵位置估计值

8）利用MN×1维发射阵排列数据X的自相关协方差矩阵R_xx行向量的线性无关性，计算合成阵列导向矢量A零空间的正交投影算子Π_θ,φ；

9）利用合成阵列导向矢量A零空间的正交投影算子Π_θ,φ构造代价函数f(θ_m,φ_n)＝1/b_mn ^HΠ_θ,φb_mn，其中 $b_{\mathrm{mn}}=a_r\left({\mathrm{\φ}}_n\right)\⊗a_t\left({\mathrm{\θ}}_m\right),m,n=1,\·\·\·,P;$

10）搜索目标相对于发射阵位置估计值和目标相对于接收阵位置估计值取代价函数f(θ_m,φ_n)模值最大时所对应的目标相对于发射阵位置估计值和目标相对于接收阵位置估计值即为同一目标相对于发射阵位置和接收阵位置。

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

(1)现有的双基地MIMO雷达角度估计技术大多采用接收天线接收的数据直接构建自相关协方差矩阵，噪声抑制效果较差，角度估计精度不高。本发明分别将接收天线接收的数据按照发射阵和接收阵排列，构建互相关协方差矩阵，大大减弱了加性噪声对角度估计的影响，提高估计精度。

(2)现有双基地MIMO雷达角度估计技术是通过对协方差矩阵进行特征分解或奇异值分解来获取噪声子空间和信号子空间，其实现复杂度高、运算量大。本发明根据互相关协方差矩阵行向量的线性无关性，求解噪声子空间，避免了特征分解带来的大运算量，加快了双基地MIMO雷达角度估计速度。

附图说明

图1是本发明的实现流程图；

图2是本发明中发射阵导向矢量零空间的正交投影算子的子流程图；

图3是本发明中接收阵导向矢量零空间的正交投影算子的子流程图；

图4是本发明中合成阵列导向矢量零空间的正交投影算子的子流程图；

图5是用本发明进行DOD，DOA的仿真结果；

图6是本发明与现有ESPRIT算法的估计性能比较图；

图7是本发明与现有ESPRIT算法复杂度比较图。

具体实施方式

以下参照附图对本发明的实现步骤及效果作进一步说明：

参照图1，本发明的具体实现步骤如下：

步骤1，对雷达回波信号进行匹配滤波，实现通道分离。

本发明的雷达系统为双基地雷达，发射站与接收站分别是由M、N个阵元构成阵元间距为半波长的均匀线阵，其中M、N为大于零的整数；各个发射子阵发射同载频的正交信号，且多普勒频率对信号正交性无影响，目标与天线间的距离远远大于发射和接收天线阵元间距。

对同一个距离单元内P个不相关的信源，第i个信源的位置为(θ_i,φ_i)，其中θ_i为目标相对于发射阵位置DOD，φ_i为目标相对于接收阵位置DOA，接收端信号依次经过混频、低通滤波和通道分离后输出匹配滤波后的MN×1维发射阵排列数据X为：

X＝[x₁₁,x₁₂,…,x_1j,…,x_1M,x₂₁,x₂₂,…,x_2j,…,x_2M,…,x_k1,x_k2,…,x_kj,…,x_kM,…,x_N1,x_N2,…,x_Nj,…,x_NM]^T〈1〉

＝As(t)+n(t)

其中，x_kj为第j个发射天线发射电磁波经目标散射后第k个接收天线接收到的数据，k＝1,…,N,j＝1,…,M；A＝[b₁,b₂,…,b_i,…,b_P]为合成阵列导向矢量，是MN×P维矩阵，i＝1,…，P，a_r(φ_i)＝[1,exp(jπsin(φ_i)),…,exp(jπ(N-1)sin(φ_i))]为接收阵导向矢量,a_t(θ_i)＝[1,exp(jπsin(θ_i)),…,exp(jπ(M-1)sin(θ_i))]为发射阵导向矢量，为Kronecker积，s(t)＝[s₁(t),s₂(t),…,s_P(t)]为发射信号矢量,s_i(t)＝α_iexp(j2πf_di),α_p和f_dp为第i个目标的幅度和多普勒频率，n(t)是均值为零方差为σ²I_MM复高斯白噪声矢量，I_NM是MN×MN维单位阵，T为转置。

步骤2，将匹配滤波后的MN×1维发射阵排列数据X按照接收阵形成数据。

2a)将匹配滤波后的MN×1维发射阵排列数据X，按照接收阵排列形成NM×1维接收阵排列数据Y：Y＝[x₁₁，x₂₁，…，x_k1，…x_N1，x₁₂，x₂₂，…，x_k2，…，x_N2，…，x_1j，x_2j，…x_kj，…，x_Nj，…，x_1M，x_2M，…，x_kM，…，x_NM]^T〈2〉

2b)对式〈1〉和〈2〉进一步化简，得到化简后的MN×1维发射阵排列数据X和NM×1维接收阵排列数据Y：

X＝[x′₁,x′₂,…,x'_k,…,x'_N]^T

Y＝[y′₁,y₂′,…,y_j',…,y_M]′^T〈3〉

其中，x'_k＝x_k1,x_k2,…,x_kj,…,x_kM,y'_j＝x_1j,x_2j,…,x_kj,…,x_Nj,k＝1,…,N,j＝1,…,M，x_kj为第j个发射天线发射电磁波经目标散射后第k个接收天线接收到的数据。

步骤3，用MN×1维发射阵排列数据X和NM×1维接收阵排列数据Y分别构造发射阵排列数据X与接收阵排列数据Y互协方差矩阵和接收阵排列数据Y与发射阵排列数据X协方差矩阵。

3a)用MN×1维发射阵排列数据X中的M个元素和NM×1维接收阵排列数据Y中的N个元素构成第一互相关矩阵R_xykj，用NM×1维接收阵排列数据Y中的N个元素和MN×1维发射阵排列数据X中的M个元素构成第三互相关矩阵R_yxjk：

R_xykj＝E{x_'ky'_j ^H}＝a_rk(φ)a_t(θ)R_s(a_tj(θ)a_r(φ))^H

R_yxjk＝E{y'_jx'_k ^H}＝a_tj(θ)a_r(φ)R_s(a_rk(φ)a_t(θ))^H〈4〉其中k＝1,…，N,j＝1,…,MR_s＝E{s(t)s(t)^H},a_rk(φ)＝exp(jπ(k-1)sin(φ)),a_tj(θ)＝exp(jπ(j-1)sin(θ))，H为共轭转置，E为数学期望；

3b)用MN×1维发射阵排列数据X中的M个元素共轭反转构成MN×1维发射阵排列共轭反转数据用NM×1维接收阵排列数据Y中的N个元素共轭反转构成NM×1维接收阵排列共轭反转数据

$\tilde{X}={[x_N^{\′},\·\·\·,x_k^{\′},\·\·\·x_2^{\′},x_1^{\′}]}^H=J_{\mathrm{NM}}{X}^{*}$

$\tilde{Y}={[y_M^{\&prime;},\&CenterDot;\&CenterDot;\&CenterDot;,y_j^{\&prime;},\&CenterDot;\&CenterDot;\&CenterDot;,y_2^{\&prime;},y_1^{\&prime;}]}^H=J_{\mathrm{NM}}{Y}^{*}---<5>$

其中J_NM为NMxNM维置换矩阵，()^*为共轭算子。

3c)用MN×1维发射阵排列共轭反转数据中的M个元素和NM×1维接收阵排列数据Y中的N个元素构成第二互相关矩阵用NM×1维接收阵排列共轭反转数据中的N个元素和MN×1维发射阵排列数据X中的M个元素构成互第四互相关矩阵

$R_{\tilde{x}\mathrm{ykj}}=E\left\{{\tilde{x}}_k^{\&prime;}{{y^{\&prime;}}_j}^T\right\}=J_{\mathrm{NM}}{R}_{\mathrm{xykj}}^{*}=a_{\mathrm{rk}}\left(\mathrm{\&phi;}\right)a_t\left(\mathrm{\&theta;}\right)D_{\mathrm{\&theta;}}^{-(M-1)}{R}_s^{*}{\left[a_{\mathrm{tj}}\left(\mathrm{\&theta;}\right)a_r\left(\mathrm{\&phi;}\right)\right]}^T$

$R_{\tilde{y}\mathrm{xjk}}=E\left\{{{\tilde{y}}^{\&prime;}}_j{x_k^{\&prime;}}^T\right\}=J_{\mathrm{NM}}{R}_{\mathrm{yxjk}}^{*}=a_{\mathrm{tj}}\left(\mathrm{\&theta;}\right)a_r\left(\mathrm{\&phi;}\right)D_{\mathrm{\&phi;}}^{-(N-1)}{R}_s^{*}{\left[a_{\mathrm{rk}}\left(\mathrm{\&phi;}\right)a_t\left(\mathrm{\&theta;}\right)\right]}^T$

其中D_θ＝diag{exp[jπsin(θ₁)],…,exp[jπsin(θ_P)]},T为转置，

${\tilde{x}}_k^{\&prime;}={x^{*}}_{(N-k+1)M},{x^{*}}_{(N-k+1)(M-1),}\&CenterDot;\&CenterDot;\&CenterDot;,{x^{*}}_{(N-k+1)j},\&CenterDot;\&CenterDot;\&CenterDot;,{x^{*}}_{(N-k+1)1},{x^{*}}_{(N-k+1)j}$ 为第j个发射天线发射电磁波经目标散射后第N-k+1个接收天线接收到数据的共轭，D_φ＝diag{exp[jπsin(φ₁)],…,exp[jπsin(φ_P)]}， ${\tilde{y}}_j^{\&prime;}={x^{*}}_{N(M-j+1)},{x^{*}}_{(N-1)(M-j+1)},\&CenterDot;\&CenterDot;\&CenterDot;,{x^{*}}_{k(M-j+1)},\&CenterDot;\&CenterDot;\&CenterDot;,{x^{*}}_{1(M-j+1)},{x^{*}}_{k(M-j+1)}$ 为第M-j+1个发射天线发射电磁波经目标散射后第k个接收天线接收到数据的共轭；

3d)分别构造发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy和接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx：

$\begin{array}{c}{R}_{\mathrm{xy}}=[R_{\mathrm{xy}11},R_{\tilde{x}y11},R_{\mathrm{xy}22},R_{\tilde{x}y22},\&CenterDot;\&CenterDot;\&CenterDot;,R_{\mathrm{xy\&Gamma;\&Gamma;}},R_{\tilde{x}\mathrm{y\&Gamma;\&Gamma;}}]\\ =a_t\left(\mathrm{\&theta;}\right)\{a_{r1}\left(\mathrm{\&phi;}\right)R_s{\left(a_{t1}\left(\mathrm{\&theta;}\right)a_r\left(\mathrm{\&phi;}\right)\right)}^H{a}_{r1}\left(\mathrm{\&phi;}\right){D_{\mathrm{\&theta;}}}^{-(M-1)}{R}_s^{*}{\left[a_{t1}\left(\mathrm{\&theta;}\right)a_r\left(\mathrm{\&phi;}\right)\right]}^T,\&CenterDot;\&CenterDot;\&CenterDot;\\ ,a_{\mathrm{r\&Gamma;}}\left(\mathrm{\&phi;}\right)R_s{\left(a_{\mathrm{t\&Gamma;}}\left(\mathrm{\&theta;}\right)a_r\left(\mathrm{\&phi;}\right)\right)}^H,a_{\mathrm{r\&Gamma;}}\left(\mathrm{\&phi;}\right){D_{\mathrm{\&theta;}}}^{-(M-1)}{R}_s^{*}{\left[a_{\mathrm{t\&Gamma;}}\left(\mathrm{\&theta;}\right)a_r\left(\mathrm{\&phi;}\right)\right]}^T\}\end{array}$

其中Γ＝min(M,N)。

步骤4，利用R_xy行向量的线性无关性，计算发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ。

现有计算导向矢量零空间正交投影算子的算法主要包括两大类：一类是特征值分解或奇异值分解算法；另一类是无特征值分解算法。本发明采用无特征值分解算法，利用R_xy行向量的线性无关性计算发射阵导向矢量a_t(θ)零空间的正交投影算子。

参照图2，本步骤的具体实现如下：

4a)假定P<M，将发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy分为两个子阵：

$R_{\mathrm{xy}}={\left[\begin{array}{c}{R}_{\mathrm{xy}}^1\\ R_{\mathrm{xy}}^2\end{array}\right]}_{\}M-P}^{\}P}---<8>$

其中和分别是P×2NΓ维子阵和(M-P)×2NΓ维子阵,由于发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy具有满秩，故R_xy行向量线性无关，存在子阵和子阵之间的线性算子P_θ，其关系式表示为：

$R_{\mathrm{xy}}^2=P_{\mathrm{\&theta;}}^H{R}_{\mathrm{xy}}^1---<9>$

进而得到线性算子P_θ为：

$P_{\mathrm{\&theta;}}\text{=}{\left[R_{\mathrm{xy}}^1{\left(R_{\mathrm{xy}}^1\right)}^H\right]}^{-1}{R}_{\mathrm{xy}}^1{R}_{\mathrm{xy}}^2---<10>$

其中()^-1为求逆算子，H表示共轭转置；

4b)利用线性算子P_θ，构建发射投影矩阵Q_θ：

$Q_{\mathrm{\&theta;}}=[P_{\mathrm{\&theta;}}^T,-{I_{M-P}}^T],---<11>$

其中,I_M-P为(M-P)×(M-P)维单位阵，T表示转置；

4c)根据发射投影矩阵Q_θ与发射阵导向矢量a_t(θ)的正交性，得到如下关系式：

$Q_{\mathrm{\&theta;}}^H{a}_t\left(\mathrm{\&theta;}\right)=O_{(M-P)\&times;P},---<12>$

其中,Ο_(M-P)×P为(M-P)×P维零矩阵；

4d)根据式〈12〉，利用正交投影公式得到发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ：

${\mathrm{\&Pi;}}_{\mathrm{\&theta;}}=Q_{\mathrm{\&theta;}}{\left(Q_{\mathrm{\&theta;}}^H{Q}_{\mathrm{\&theta;}}\right)}^{-1}{Q}_{\mathrm{\&theta;}}^H.---<13>$

步骤5，利用R_yx行向量的线性无关性计算接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ。

现有计算导向矢量零空间正交投影算子的算法主要包括两大类：一类是特征值分解或奇异值分解算法；另一类是无特征值分解算法。本发明采用无特征值分解算法，利用R_yx行向量的线性无关性计算接收阵导向矢量a_r(φ)零空间的正交投影算子。

参照图3，本步骤的具体实现如下：

5a)假定P<N，将接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx分为两个子阵：

$R_{\mathrm{yx}}={\left[\begin{array}{c}{R}_{\mathrm{yx}}^1\\ R_{\mathrm{yx}}^2\end{array}\right]}_{\}N-P}^{\}P}---<14>$

其中和分别是P×2MΓ维子阵和(N-P)×2MΓ维子阵,由于接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx具有满秩，故R_yx行向量线性无关，存在子阵和子阵之间的线性算子P_φ，其关系式表示为：

$R_{\mathrm{yx}}^2=P_{\mathrm{\&phi;}}^H{R}_{\mathrm{yx}}^1---<15>$

进而得到线性算子P_φ为：

$P_{\mathrm{\&phi;}}={\left[R_{\mathrm{yx}}^1{\left(R_{\mathrm{yx}}^1\right)}^H\right]}^{-1}{R}_{\mathrm{yx}}^1{R}_{\mathrm{yx}}^2---<16>$

其中()^-1为求逆算子，H表示共轭转置；

5b)利用线性算子P_φ，构建接收投影矩阵Q_φ：

$Q_{\mathrm{\&phi;}}=[P_{\mathrm{\&phi;}}^T,-{I_{N-P}}^T],---<17>$

其中,I_N-P为(N-P)×(N-P)维单位阵，T表示转置。

5c)根据接收投影矩阵Q_φ与接收阵导向矢量a_r(φ)的正交性，可得如下关系式：

$Q_{\mathrm{\&phi;}}^H{a}_r\left(\mathrm{\&phi;}\right)=O_{(N-P)\&times;P},---<18>$

其中,Ο_(N-P)×P为(N-P)×P维零矩阵；

5d)根据式〈18〉，利用正交投影公式，得到接收阵导向矢量a_r(φ)零空间的正交投影算子Πφ：

${\mathrm{\&Pi;}}_{\mathrm{\&phi;}}=Q_{\mathrm{\&phi;}}{\left(Q_{\mathrm{\&phi;}}^H{Q}_{\mathrm{\&phi;}}\right)}^{-1}{Q}_{\mathrm{\&phi;}}^H.---<19>$

步骤6，利用发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ，分别构造发射阵求根多项式f(z_t)和接收阵求根多项式f(z_r)：

$f\left(z_t\right)=a_t{\left(z_t^{-1}\right)}^T{\mathrm{\&Pi;}}_{\mathrm{\&theta;}}{a}_t\left(z_t\right),$ 其中 $a_t\left(z_t\right)={[1,z_t,z_t^2,\&CenterDot;\&CenterDot;\&CenterDot;,z_t^{M-1}]}^T,z_t=e^{\mathrm{j\&pi;}\sin \left(\mathrm{\&theta;}\right)};$

$f\left(z_r\right)=a_r{\left(z_r^{-1}\right)}^T{\mathrm{\&Pi;}}_{\mathrm{\&theta;}}{a}_r\left(z_r\right),$ 其中 $a_r\left(z_r\right)={[1,z_r{z}_r^2,\&CenterDot;\&CenterDot;\&CenterDot;,z_r^{M-1}]}^T,z_r=e^{\mathrm{j\&pi;}\sin \left(\mathrm{\&phi;}\right)}.---<20>$

步骤7，分别求解发射阵求根多项式单位圆上的根z_t和接收阵求根多项式单位圆上的根z_r，得到目标相对于发射阵位置估计值和目标相对于接收阵位置估计值分别表述如下：

$\widehat{\mathrm{\&theta;}}=\arcsin (\arg \left(z_t\right)/\mathrm{\&pi;}),$

步骤8，利用MN×1维发射阵排列数据X的自相关协方差矩阵R_xx行向量的线性无关性，计算合成阵列导向矢量A零空间的正交投影算子Π_θ,φ。

现有计算导向矢量零空间正交投影算子的算法主要包括两大类：一类是特征值分解或奇异值分解算法；另一类是无特征值分解算法。本发明采用无特征值分解算法，利用R_xx行向量的线性无关性计算合成阵列导向矢量A零空间的正交投影算子。

参照图4，本步骤的具体实现如下：

8a)利用MN×1维发射阵排列数据X构造自相关协方差矩阵R_xx，R_xx＝E{xx^H}，E为数学期望；

8b)假定P<MN，将自相关协方差矩阵R_xx划分为PxP维子阵和(NM-P)×P维子阵

$R_{\mathrm{xx}}={\left[\begin{array}{c}{R}_{\mathrm{xx}}^1\\ R_{\mathrm{xx}}^2\end{array}\right]}_{\}\mathrm{NM}-P}^{\}P}---<22>$

由于自相关协方差矩阵R_xx具有满秩，故R_xx行向量线性无关，得到的关系式，存在子阵与子阵之间的线性算子P_θ,φ，其关系式表示为：

$R_{\mathrm{xx}}^2=P_{\mathrm{\&theta;},\mathrm{\&phi;}}^H{R}_{\mathrm{xx}}^1---<23>$

进而得到线性算子P_θ,φ为：

$P_{\mathrm{\&theta;},\mathrm{\&phi;}}={\left[R_{\mathrm{xx}}^1{\left(R_{\mathrm{xx}}^1\right)}^H\right]}^{-1}{R}_{\mathrm{xx}}^1{R}_{\mathrm{xx}}^2---<24>$

其中()^-1为求逆算子，H表示共轭转置；

8c)利用子阵与子阵之间的线性算子P_θ,φ，构建合成阵列投影矩阵Q_θ,φ：

$Q_{\mathrm{\&theta;},\mathrm{\&phi;}}={[P_{\mathrm{\&theta;},\mathrm{\&phi;}}^T,-I_{\mathrm{NM}-P}]}^T---<25>$

其中I_NM-P为(NM-P)×(NM-P)维单位阵；

8d)根据合成阵列投影矩阵Q_θ,φ和合成阵列导向矢量A的正交性，得到如下关系式：

$Q_{\mathrm{\&theta;},\mathrm{\&phi;}}^{H}A=O_{(\mathrm{NM}-P)\&times;P}---<26>$

其中Ο_(NM-P)×P为(NM-P)×P维零矩阵；

8e)根据步骤8d)得到的关系式和正交投影公式，得到合成阵列导向矢量A零空间的正交投影算子Π_θ,φ：

${\mathrm{\&Pi;}}_{\mathrm{\&theta;},\mathrm{\&phi;}}=Q_{\mathrm{\&theta;},\mathrm{\&phi;}}{\left(Q_{\mathrm{\&theta;},\mathrm{\&phi;}}^H{Q}_{\mathrm{\&theta;},\mathrm{\&phi;}}\right)}^{-1}{Q}_{\mathrm{\&theta;},\mathrm{\&phi;}}^H.---<27>$

步骤9，利用合成阵列导向矢量A零空间的正交投影算子Π_θ,φ构造代价函数：

f(θ_m,φ_n)＝1/b_mn ^HΠ_θ,φb_mn〈28〉

其中 $b_{\mathrm{mn}}=a_r\left({\mathrm{\&phi;}}_n\right)\&CircleTimes;a_t\left({\mathrm{\&theta;}}_m\right),m,n=1,\&CenterDot;\&CenterDot;\&CenterDot;,P.$

步骤10，角度配对。

搜索目标相对于发射阵位置估计值和目标相对于接收阵位置估计值取代价函数f(θ_m,φ_n)模值最大时所对应的目标相对于发射阵位置估计值和目标相对于接收阵位置估计值即为同一目标相对于发射阵位置和接收阵位置。

至此完成基于互相关矩阵的双基MIMO雷达角度估计。

本发明的效果通过以下仿真试验进一步说明：

1.仿真条件

仿真参数如下：发射阵及接收阵均为中心频率对应的半波长均匀线阵，两个不相干目标位于(θ₁,φ₁)=(5°,-5°)和(θ₂,φ₂)=(15°,5°)，快拍数为100，发射阵天线个数M=6,接收阵天线个数N=8。角度估计的均方根误差RMSE采用公式：其中和分别为第i个目标相对于发射阵位置估计值和目标相对于接收阵位置估计值，θ_i和φ_i分别为第i个目标相对于发射阵位置和目标相对于接收阵位置，P为目标个数，E为数学期望。

2.仿真内容：

仿真1，采用本发明在信噪比为10dB时对目标角度进行100次角度估计，仿真结果如图5所示。

仿真2，采用本发明和现有ESPRIT方法对仿真条件下的两个目标进行角度估计，仿真目标角度估计均方根误差随信噪比变化，每个信噪比下进行100次蒙特卡洛仿真试验，仿真结果如图6所示。

仿真3，采用本发明和现有ESPRIT方法对仿真条件下的两个目标进行角度估计，仿真目标角度估计所需时间随发射接收阵元数变化，信噪比为10dB时每个阵元数下进行100次蒙特卡洛仿真实验，仿真结果如图7所示。

3.仿真分析

从图5可以看出，本发明能够正确估计目标角度，并正确完成配对。

从图6可以看出，采用本发明估计目标角度的均方误差随信噪比增加而降低，与现有ESPRIT算法相比，采用本发明可以更高的精度实现双基地MIMO雷达对目标的角度的估计。

从图7可以看出，采用本发明估计目标角度所消耗时间随阵元数增加而增加，相比现有ESPRIT算法，当阵元数较小时本发明角度估计所需时间较现有ESPRIT算法长，但随着阵元数的增加，现有ESPRIT算法消耗时间明显增长，即本发明运算量明显小于现有ESPRIT算法。

Claims (7)

1.一种基于互相关矩阵的双基MIMO雷达角度估计方法，包括如下步骤：

1)M个发射天线发射电磁波照射空间P个不相关目标，N个接收天线接收目标散射回波，并对回波信号依次进行混频、低通滤波、通道分离，得到MN×1维发射阵排列数据X＝As(t)+n(t)其中A＝[b₁,b₂,…,b_i,…,b_P]，为合成阵列导向矢量，i＝1,…,P， 为Kronecker积；a_r(φ_i)和a_t(θ_i)分别为第i个目标的接收阵和发射阵导向矢量，θ_i为第i个目标相对于发射阵位置，φ_i为第i个目标相对于接收阵位置；s(t)为发射信号矢量，n(t)是复高斯白噪声矢量；M、N为大于零的整数；

2)将MN×1维发射阵排列数据X按照接收阵排列，得到NM×1维接收阵排列数据Y；

3)构造发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy和接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx；

3a)用MN×1维发射阵排列数据X中的M个元素和NM×1维接收阵排列数据Y中的N个元素构成第一互相关矩阵R_xykj,k＝1,…,N,j＝1,…,M；用MN×1维发射阵排列数据X中的M个元素共轭反转和NM×1维接收阵排列数据Y中的N个元素构成第二互相关矩阵将这两个互相关矩阵按照行排列，构成发射阵排列数据X与接收阵排列数据Y互相关协方差矩阵R_xy的一个子阵，该互相关协方差矩阵R_xy的子阵数小于等于M和N中的最小值；

3b)用NM×1维接收阵排列数据Y中的N个元素和MN×1维发射阵排列数据X中的M个元素构成第三互相关矩阵R_yxjk；用NM×1维接收阵排列数据Y中的N个元素共轭反转和MN×1维发射阵排列数据X中的M个元素构成第四互相关矩阵将这两个互相关矩阵按照行排列，构成接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx的一个子阵，该互相关协方差矩阵R_yx的子阵数小于等于M和N中的最小值；

4)将所述发射阵排列数据X与接收阵排列数据Y的互相关协方差矩阵R_xy划分为P×2NΓ维子阵和(M-P)×2NΓ维子阵将所述接收阵排列数据Y与发射阵排列数据X互相关协方差矩阵R_yx划分为P×2MΓ维子阵和(N-P)×2MΓ维子阵

其中Γ＝min(M,N)，M、N分别为发射阵和接收阵天线个数；

5)利用所述R_xy行向量的线性无关性和所述R_yx行向量的线性无关性，计算发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ：

6)利用发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ，分别构造发射阵求根多项式f(z_t)和接收阵求根多项式f(z_r)；

7)分别求解发射阵求根多项式单位圆上的根z_t和接收阵求根多项式单位圆上的根z_r，得到目标相对于发射阵位置估计值和目标相对于接收阵位置估计值

8)利用MN×1维发射阵排列数据X的自相关协方差矩阵R_xx行向量的线性无关性，计算合成阵列导向矢量A零空间的正交投影算子Π_θ,φ；

9)利用合成阵列导向矢量A零空间的正交投影算子Π_θ,φ构造代价函数f(θ_m,φ_n)＝1/b_mn ^HΠ_θ,φb_mn，其中

10)搜索目标相对于发射阵位置估计值和目标相对于接收阵位置估计值取代价函数f(θ_m,φ_n)模值最大时所对应的目标相对于发射阵位置估计值和目标相对于接收阵位置估计值即为同一目标相对于发射阵位置和接收阵位置。

2.根据权利要求1所述的基于互相关矩阵的双基MIMO雷达角度估计方法，其中所述步骤3a)中的发射阵排列数据X与接收阵排列数据Y的互相关协方差矩阵R_xy，其公式如下：

其中，R_xykj＝E{x'_ky'_j ^H}，x′_k＝x_k1，x_k2，…，x_kj，…，x_kM，k＝1,…,N,j＝1,…,M，y'_j＝x_1j,x_2j,…,x_kj,…,x_Nj，x_kj为第j个发射天线发射电磁波经目标散射后第k个接收天线接收到的数据， 为第j个发射天线发射电磁波经目标散射后第N-k+1个接收天线接收到数据的共轭，Γ＝min(M,N)，M、N分别为发射阵和接收阵天线个数，E为数学期望，H为共轭转置，T为转置，*为共轭。

3.根据权利要求1所述的基于互相关矩阵的双基MIMO雷达角度估计方法，其中所述步骤3b)中的接收阵排列数据Y与发射阵排列数据X的互相关协方差矩阵R_yx，其公式如下：

其中,R_yxjk＝E{y'_jx'_k ^H}，y'_j＝x_1j,x_2j,…,x_kj,…,x_Nj，k＝1,…,N,j＝1,…,M，x'_k＝x_k1,x_k2,…,x_kj,…,x_kM，x_kj为第j个发射天线发射电磁波经目标散射后第k个接收天线接收到的数据，x^* _(N-1)(M-j+1),…,x^* _k(M-j+1),…,x^* _1(M-j+1)，x^* _k(M-j+1)为第M-j+1个发射天线发射电磁波经目标散射后第k个接收天线接收到数据的共轭，Γ＝min(M,N)，M、N分别为发射阵和接收阵天线个数，E为数学期望，H为共轭转置，T为转置，*为共轭。

4.根据权利要求1所述的基于互相关矩阵的双基MIMO雷达角度估计方法，其中步骤5)所述的利用所述R_xy行向量的线性无关性和所述R_yx行向量的线性无关性，计算发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ，按如下步骤进行：

5a)利用所述R_xy行向量的线性无关性和所述R_yx行向量的线性无关性得到 这两个关系式，进而得到：

其中()^-1为求逆算子，H表示共轭转置，P_θ为所述子阵和子阵之间的线性算子，P_φ为所述子阵和子阵之间的线性算子；

5b)利用步骤5a)得到的两个算子P_θ和P_φ，分别构建发射投影矩阵Q_θ和接收投影矩阵Q_φ：

其中I_M-P为(M-P)×(M-P)维单位阵，I_N-P为(N-P)×(N-P)维单位阵，T表示转置；

5c)根据发射投影矩阵Q_θ与发射阵导向矢量a_t(θ)的正交性，以及接收投影矩阵Q_φ与接收阵导向矢量a_r(φ)的正交性，可得如下关系式：

其中，Ο_(N-P)×P为(N-P)×P维零矩阵，Ο_(M-P)×P为(M-P)×P维零矩阵；

5d)根据所述步骤5c)的两个关系式得到发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ：

5.根据权利要求1所述的基于互相关矩阵的双基MIMO雷达角度估计方法，其中步骤6)所述的用发射阵导向矢量a_t(θ)零空间的正交投影算子Π_θ，以及接收阵导向矢量a_r(φ)零空间的正交投影算子Π_φ，分别构造发射阵求根多项式f(z_t)和接收阵求根多项式f(z_r)，其公式如下：

其中z_t＝e^jπsin(θ)；

其中z_r＝e^jπsin(φ)。

6.根据权利要求1所述的基于互相关矩阵的双基MIMO雷达角度估计方法，其中所述步骤7)中目标相对于发射阵位置估计值和目标相对于接收阵位置估计值分别表述如下：

7.根据权利要求1所述的基于互相关矩阵的双基MIMO雷达角度估计方法，其中步骤8)所述的利用MN×1维发射阵排列数据X的自相关协方差矩阵R_xx行向量的线性无关性，计算合成阵列导向矢量A零空间的正交投影算子Π_θ,φ，按如下步骤进行：

8a)利用MN×1维发射阵排列数据X构造自相关协方差矩阵R_xx，R_xx＝E{xx^H}，并将其划分为PxP维子阵和(NM-P)×P维子阵其中E为数学期望；

8b)利用所述R_xx行向量的线性无关性，得到的关系式，进而得到所述子阵与所述子阵之间的线性算子P_θ,φ：

8c)利用步骤8b)得到的子阵与子阵之间的线性算子P_θ,φ，构建合成阵列投影矩阵Q_θ,φ：

其中I_NM-P为(NM-P)×(NM-P)维单位阵；

8d)根据合成阵列投影矩阵Q_θ,φ和合成阵列导向矢量A的正交性，得到的关系式，其中Ο_(NM-P)×P为(NM-P)×P维零矩阵；

8e)根据步骤8d)得到的关系式，得到合成阵列导向矢量A零空间的正交投影算子：