CN104392114A - 一种基于空时数据的高分辨目标方位估计方法 - Google Patents

Abstract

本发明公开了一种基于空时数据的高分辨目标方位估计方法，属于雷达技术领域。本发明通过采样记录天线各传感器通道接收的数据，将采样数据排成列向量的形式，在时间域、空间域和延迟域形成相关矩阵，对相关矩阵进行去噪处理，利用循环特征分解方法计算信号子空间，从信号和噪声子空间的关系求解噪声子空间，得到噪声子空间的投影矩阵，基于投影矩阵反解目标方位和俯仰信息，避免了经典MUSIC方法在区分信号子空间和噪声子空间时性能出现急剧下降的问题，提高了确定信源的方位角时的准确性。

Description

一种基于空时数据的高分辨目标方位估计方法

技术领域

本发明涉及雷达技术领域，特别涉及一种基于空时数据的高分辨目标方位估计方法。

背景技术

空间定位功能是天线应用的关键技术之一。为了追求更高的定位精度，根据方位估计与时间信号的频率估计的相似性，许多时域非线性定位方法推广到空间方位估计中便产生了所谓的高分辨方位估计方法。

美国的Schmidt提出了著名的多重信号分类(MUSIC)方法，实现了现代高分辨方位估计技术的飞跃，从而促进了特征子空间类方法的兴起，该方法通过数据矩阵的奇异值分解(SVD)或空间协方差矩阵的特征值分解(EVD)来获得相互正交的信号子空间和噪声子空间，然后利用天线阵列导向矢量和子空间的关系来构造信号参数的估计算法，得到方位估计的渐近无偏估计，突破了以往方位估计算法中阵列孔径对参数估计性能的瑞利限制，目标分辨能力可以达到波束宽度的1/3～1/5。

在实现本发明的过程中，发明人发现现有技术至少存在以下问题：

在实际应用中，信噪比低于零以及采样数较小时，经典MUSIC方法在区分信号子空间和噪声子空间时性能出现急剧下降的问题，导致在确定信源的方位角时准确性下降。

发明内容

为了解决现有技术的问题，本发明提供了一种基于空时数据的高分辨目标方位估计方法，所述方法包括：

步骤一，通过传感器阵列对信号进行采样，获取到第一信号矢量x₁(t)、第二信号矢量x₂(t)，其中x₁(t)＝A₁s(t)+n₁(t)，x₂(t)＝A₁Φs(t)+n₂(t)；

步骤二，在开始采样后的第一个时刻，构造与所述第一信号矢量、所述第二信号矢量的对应的矩阵

$R_{11}^{\left(1\right)}=A_1{R}_s\left(1\right)A_1^H+{\mathrm{\σ}}_1^2\left(1\right)I,$

$R_{22}^{\left(1\right)}=A_1{\mathrm{\ΦR}}_s\left(1\right){\mathrm{\Φ}}^H{A}_1^H+{\mathrm{\σ}}_2^2\left(1\right)I,$

其中R_s为自相关矩阵，为所述采样过程中的噪声，在开始采样后的第二个时刻，构造与所述第一信号矢量、所述第二信号矢量对应的矩阵

$R_{11}^{\left(2\right)}=A_1{R}_s\left(2\right)A_1^H+{\mathrm{\σ}}_1^2\left(2\right)I,$

$R_{22}^{\left(2\right)}=A_1{\mathrm{\ΦR}}_s\left(2\right){\mathrm{\Φ}}^H{A}_1^H+{\mathrm{\σ}}_2^2\left(2\right)I,$

依次在每个时刻均构造两个相关矩阵，并根据时间顺序，将全部相关矩阵组合为相关矩阵组其中k代表时刻数，Λ_m＝Φ^m-1， $B_m^{\left(k\right)}=A_m{R}_s^H\left(k\right);$

步骤三，将所述全部相关矩阵进行特征分解，对所述全部相关矩阵中的噪声进行估计，得到每个噪声的估计值 对所述全部特征矩阵进行去噪，得到去噪相关矩阵组

步骤四，构造初始矩阵U(0)，建立第一代价函数

$\tilde{C}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{m=1}{\overset{2}{\mathrm{\Σ}}}\left\{{\left[R_{\mathrm{mm}}^{\left(k\right)}\right]}^{H}U(l-1)U^H(l-1)R_{\mathrm{mm}}^{\left(k\right)}\right\},$

其中，l代表循环次数，对所述第一代价函数进行特征值分解，得到第一分解式

${\tilde{C}}_{\mathrm{eig}}=\stackrel{\‾}{V}(l-1)\stackrel{\‾}{D}(l-1){\stackrel{\‾}{V}}^H(l-1),$

其中，为特征值矩阵,为特征向量矩阵,令矩阵V(l-1)为矩阵的前P个大特征值对应的特征向量形成的矩阵,建立第二代价函数

其中，l代表循环次数，对所述第二代价函数进行特征值分解，得到第二分解式

其中，为特征值矩阵,为特征向量矩阵,所述第一代价函数、第二代价函数为U(l-1)到U(l)的一次循环，对所述U(l)循环计算，直至满足

||U(l)U^H(l)-U(l-1)U^H(l-1)||_F＜ε

时停止循环，获取停止循环时的l；

步骤五，根据所述停止循环时的l，确定信号子空间U(l)，令得到投影矩阵根据所述投影矩阵得到角度估计

$P\left(\widehat{\mathrm{\θ}}\right)=\frac{1}{\left|a^H\left(\widehat{\mathrm{\θ}}\right)\right[I-{\tilde{U}}_s{\tilde{U}}_s^H\left]a\left(\widehat{\mathrm{\θ}}\right)\right|}.$

本发明提供的技术方案带来的有益效果是：

通过采样记录天线各传感器通道接收的数据，将采样数据排成列向量的形式，在时间域、空间域和延迟域形成相关矩阵，对相关矩阵进行去噪处理，利用循环特征分解方法计算信号子空间，从信号和噪声子空间的关系求解噪声子空间，得到噪声子空间的投影矩阵，基于投影矩阵反解目标方位和俯仰信息，避免了经典MUSIC方法在区分信号子空间和噪声子空间时性能出现急剧下降的问题，提高了确定信源的方位角、俯仰角时准确性。

附图说明

为了更清楚地说明本发明的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种基于空时数据的高分辨目标方位估计方法示意图；

图2为本发明提供的传感器的结构示意图；

图3为本发明提供的检测平面示意图；

图4(a)为本发明提供的子空间误差随信噪比变化情况示意图；

图4(b)为本发明提供的子空间误差随快拍数变化情况示意图；

图4(c)为本发明提供的子空间误差随传感器个数变化情况示意图；

图5(a)为本发明提供的MUSIC方法的空间谱随信噪比变化示意图；

图5(b)为本发明提供的本发明方法的空间谱随信噪比变化示意图；

图5(c)为本发明提供的本发明的谱峰和MUSIC方法的谱峰之差随信噪比变化示意图；

图6(a)为本发明提供的MUSIC方法的空间谱随快拍数变化示意图；

图6(b)为本发明提供的空间谱随快拍数变化示意图；

图6(c)为本发明提供的谱峰和MUSIC方法的谱峰之差随快拍数变化示意图；

图7(a)为本发明提供的MUSIC方法的空间谱随传感器个数变化示意图；

图7(b)为本发明提供的空间谱随传感器个数变化示意图；

图7(c)为本发明提供的谱峰和MUSIC方法的谱峰之差随传感器个数变化示意图。

具体实施方式

为使本发明的结构和优点更加清楚，下面将结合附图对本发明的结构作进一步地描述。

实施例一

本发明提供了一种基于空时数据的高分辨目标方位估计方法，如图1所示，该装置包括：

步骤一，通过传感器阵列对信号进行采样，获取到第一信号矢量x₁(t)、第二信号矢量x₂(t)，其中x₁(t)＝A₁s(t)+n₁(t)，x₂(t)＝A₁Φs(t)+n₂(t)。

在实施中，天线由M+1个均匀线性传感器组成,它由两个具有相同结构的子阵组成，两个子阵列均含M个天线阵元，具体结构如图2所示，第1个子阵含阵元1到M，第2个子阵含阵元2到M+1。设第p个信源的方位角为θ_p。在第t(t＝1,2,…,T)次采样下所有P个信源在第m个阵元上的总接收信号为

$x_m\left(t\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{a}_{\mathrm{mp}}{s}_p\left(t\right)+n_m\left(t\right),$

其中d为相邻阵元的间距，λ为信源的波长，s_p第p个信源复包络，n_m为第m个振元上的噪声。第t次采样下第1个子阵的接收数据矢量可表示为

x₁(t)＝A₁s(t)+n₁(t)，

其中A₁＝[a₁,a₂,…,a_P]，a_p为列向量，其第i个元素为exp(j2π(i-1)θ_p/λ)，s(t)＝[s₁(t),s₂(t),…,s_P(t)]^T，n₁(t)＝[n₁(t),n₂(t),…,n_M(t)]^T，第t次采样子阵2的接收数据矢量可表示为

x₂(t)＝A₁Φs(t)+n₂(t)，

其中Φ为对角矩阵， $\mathrm{\Φ}=\mathrm{diag}[e^{{\mathrm{j\ω}}_1},e^{{\mathrm{j\ω}}_2}...,e^{{\mathrm{j\ω}}_P}],{\mathrm{\ω}}_p=\frac{2\mathrm{\π}d\sin {\mathrm{\θ}}_p}{\mathrm{\λ}},$ 噪声为

n₂(t)＝[n_M+1(t),n_M+2(t),…,n_2M(t)]^T。

步骤二，在开始采样后的第一个时刻，构造与第一信号矢量、第二信号矢量的对应的矩阵

$R_{11}^{\left(1\right)}=A_1{R}_s\left(1\right)A_1^H+{\mathrm{\σ}}_1^2\left(1\right)I,$

$R_{22}^{\left(1\right)}=A_1{\mathrm{\ΦR}}_s\left(1\right){\mathrm{\Φ}}^H{A}_1^H+{\mathrm{\σ}}_2^2\left(1\right)I,$

其中R_s为自相关矩阵，为所述采样过程中的噪声，在开始采样后的第二个时刻，构造与第一信号矢量、第二信号矢量对应的矩阵

$R_{11}^{\left(2\right)}=A_1{R}_s\left(2\right)A_1^H+{\mathrm{\σ}}_1^2\left(2\right)I,$

$R_{22}^{\left(2\right)}=A_1{\mathrm{\ΦR}}_s\left(2\right){\mathrm{\Φ}}^H{A}_1^H+{\mathrm{\σ}}_2^2\left(2\right)I,$

依次在每个时刻均构造两个相关矩阵，并根据时间顺序，将全部相关矩阵组合为相关矩阵组其中k代表时刻数，Λ_m＝Φ^m-1， $B_m^{\left(k\right)}=A_m{R}_s^H\left(k\right).$

在实施中，因为在每个采样时刻，第一子阵和第二子阵分别获取采样信号值，因此在同一采样时刻可以生成两个信号矢量矩阵，当采样过程结束后，将生成的多个信号矢量矩阵重新组合成相关矩阵组

步骤三，将全部相关矩阵进行特征分解，对全部相关矩阵中的噪声进行估计，得到每个噪声的估计值对全部特征矩阵进行去噪，得到去噪相关矩阵组

在实施中，对相关矩阵组中的每个矩阵进行特征分解，在分解后的矩阵选取M-P个小特征值，根据上述特征值的平均值对噪声进行估计，得到每个噪声的估计值从相关矩阵组中的每个矩阵中去除对应的噪声后，得到的部分去噪矩阵如下所示：

${\hat{R}}_{11}^{\left(1\right)}=R_{11}^{\left(1\right)}-{\widehat{\mathrm{\σ}}}_1^2\left(1\right)I,$

${\hat{R}}_{22}^{\left(1\right)}=R_{22}^{\left(1\right)}-{\widehat{\mathrm{\σ}}}_2^2\left(1\right)I,$

${\hat{R}}_{11}^{\left(2\right)}=R_{11}^{\left(2\right)}+{\widehat{\mathrm{\σ}}}_1^2\left(2\right)I,$

${\hat{R}}_{22}^{\left(2\right)}=R_{22}^{\left(2\right)}-{\widehat{\mathrm{\σ}}}_2^2\left(2\right)I,$

.

.

.

${\hat{R}}_{11}^{\left(K\right)}=R_{11}^{\left(K\right)}+{\widehat{\mathrm{\σ}}}_1^2\left(K\right)I,$

${\hat{R}}_{22}^{\left(K\right)}=R_{22}^{\left(K\right)}-{\widehat{\mathrm{\σ}}}_2^2\left(K\right)I,$

将上述去噪矩阵组合为去噪矩阵组

步骤四，构造初始矩阵U(0)，建立第一代价函数

$\tilde{C}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{m=1}{\overset{2}{\mathrm{\Σ}}}\left\{{\left[R_{\mathrm{mm}}^{\left(k\right)}\right]}^{H}U(l-1)U^H(l-1)R_{\mathrm{mm}}^{\left(k\right)}\right\},---\left(1\right)$

其中，l代表循环次数，对第一代价函数进行特征值分解，得到第一分解式

${\tilde{C}}_{\mathrm{eig}}=\stackrel{\‾}{V}(l-1)\stackrel{\‾}{D}(l-1){\stackrel{\‾}{V}}^H(l-1),---\left(2\right)$

其中，为特征值矩阵,为特征向量矩阵,令矩阵V(l-1)为矩阵的前P个大特征值对应的特征向量形成的矩阵,建立第二代价函数

其中，l代表循环次数，对第二代价函数进行特征值分解，得到第二分解式

其中，为特征值矩阵,为特征向量矩阵,第一代价函数、第二代价函数为U(l-1)到U(l)的一次循环，对U(l)循环计算，直至满足

||U(l)U^H(l)-U(l-1)U^H(l-1)||_F＜ε，(5)

时停止循环，获取停止循环时的l。

在实施中，对初始矩阵U(0)构造两个代价函数和对进行特征值分解后得到与相关的多项式，对进行特征分解后得到与相关的多项式，此时容易发现公式(1)至公式(2)，是从U(l-1)多项式到多项式的过程，公式(3)至(4)是从V(l-1)多项式到多项式的过程，因此，公式(1)至公式(2)、公式(3)至公式(4)为一次U(l-1)到U(l)的循环，每完成一次循环，()内的数值加1，直至公式(5)成立时，根据此时l的数值，确定矩阵U(L)。

步骤五，根据停止循环时的l，确定信号子空间U(l)，令得到投影矩阵根据投影矩阵得到角度估计

$P\left(\widehat{\mathrm{\θ}}\right)=\frac{1}{\left|a^H\left(\widehat{\mathrm{\θ}}\right)\right[I-{\tilde{U}}_s{\tilde{U}}_s^H\left]a\left(\widehat{\mathrm{\θ}}\right)\right|}.$

在实施中，为了便于理解，这里将本发明所使用的传感器阵列简化为一个质点，该质点固定在平面二维坐标系的原点位置，如图3所示，在XOY平面内，该信源相对于该质点的角度范围为-90°～90°，在XOY平面内，通过计算采样结果中峰值最高的采样信号对应的角度值，确定信源的方位。公式中的表示在步骤二中生成的第一信号矢量、第二信号矢量中矩阵A₁中的向量，在第一信号矢量、第二信号矢量中仅含有方位角一个角度，向量中也仅含有方位角这一个角度，其中对方位角θ的角度估计，而代表向量的共轭转置向量。

为进一步说明本发明的基于空时数据的目标方位估计方法较经典高分辨定位方(如MUSIC)方法的优越性，做如下仿真实验，各实验均取100次独立实验平均结果。

实验1：信号子空间误差性能分析。

信号子空间的精度是保证定位准确性的重要前提，本实验对发明中提出的多循环双向特征分解方法和经典的MUSIC使用的特征分解类方法得到的信号子空间性能进行对比。实验中两个信源来自θ₁＝-3°和θ₂＝3°，图4(a)为子空间误差随信噪比变化情况，试验中取快拍为200和500两种不同情况进行实验，由图可见，在相同快拍下，本发明所得到的子空间误差小于经典MUSIC的特征分解类方法，本发明在200次快拍下的子空间误差在信噪比为-10dB到10dB的大部分范围内甚至低于MUSIC方法在500次快拍下的误差。

图4(b)为子空间误差随快拍数变化情况，快拍数从100到1000范围内变化。试验中取信噪比为0dB和5dB两种不同情况进行实验，由图可见，在相同信噪比下，本发明所得到的子空间误差小于经典MUSIC的特征分解类方法。

图4(c)为子空间误差随传感器个数变化情况，传感器个数从8到20范围内变化。试验中取信噪比为0dB和5dB两种不同情况进行实验，由图可见，在相同的传感器个数下，本发明所得到的子空间误差小于经典MUSIC的特征分解类方法。

实验2：方位估计性能随信噪比变化。

本实验对比发明方法和经典MUSIC方法定位性能随信噪比变化情况，定位性能通常体现于谱峰函数。图5(a)为MUSIC方法的空间谱，图5(b)为本发明方法的空间谱，通过5(a)和5(b)的比较可见，MUSIC方法在-10dB信噪比下分辨不出两个信号源，而本发明方法在-10dB信噪比下可完全分辨两个信源，且得到的谱峰比较尖锐，分辨率较高。图5(c)将本文得到的谱峰和MUSIC方法的谱峰进行相减，可见得到的均是正数值，说明在-10dB到10dB范围内，本发明得到的谱峰高于MUSIC方法的谱峰，具有更好的分辨性能。高的谱峰值同时说明本发明的多循环双向特征分解方法得到的子空间更接近于真实的信号子空间，和噪声子空间具有更好的正交性。

实验3：方位估计性能随快拍数变化。

本实验对比发明方法和经典MUSIC方法定位性能随快拍数变化情况。图6(a)为MUSIC方法的空间谱，图6(b)为本发明方法的空间谱，通过6(a)和6(b)的比较可见，MUSIC方法在200和300次快拍下分辨不出两个信号源，而本发明方法在200和300次快拍下可完全分辨两个信源，且得到的谱峰比较尖锐，分辨率较高。图6(c)将本文得到的谱峰和MUSIC方法的谱峰进行相减，可见得到的均是正数值，说明在100到1000次快拍范围内，本发明得到的谱峰高于MUSIC方法的谱峰，具有更好的分辨性能。

实验4：方位估计性能随传感器个数变化。

本实验对比发明方法和经典MUSIC方法定位性能随传感器个数变化。图7(a)为MUSIC方法的空间谱，图7(b)为本发明方法的空间谱，通过7(a)和7(b)的比较可见，MUSIC方法在8个传感器下完全不能分辨两个目标，在10个传感器下分辨性能很差，本发明方法在10个传感器下可清晰地分辨两个源信号，在8个传感器下的分辨能力优于MUSIC方法10个传感器的分辨性能。图7(c)将本文得到的谱峰和MUSIC方法的谱峰进行相减，可见得到的均是正数值，说明本发明方法得到的谱峰比较尖锐，分辨率较高。在8到20个传感器变化范围内，本发明得到的谱峰高于MUSIC方法的谱峰，具有更好的分辨性能。

本实施例中提出的一种基于空时数据的高分辨目标方位估计方法，通过采样记录天线各传感器通道接收的数据，将采样数据排成列向量的形式，在时间域、空间域和延迟域形成相关矩阵，对相关矩阵进行去噪处理，利用循环特征分解方法计算信号子空间，从信号和噪声子空间的关系求解噪声子空间，得到噪声子空间的投影矩阵，基于投影矩阵反解目标方位和俯仰信息，避免了经典MUSIC方法在区分信号子空间和噪声子空间时性能出现急剧下降的问题，提高了确定信源的方位角、俯仰角时准确性。

需要说明的是：上述实施例提供的基于空时数据的高分辨目标方位估计方法进行信源定位的实施例，仅作为该方法中在实际应用中的说明，还可以根据实际需要而将上述方法在其他应用场景中使用，其具体实现过程类似于上述实施例，这里不再赘述。

以上所述仅为本发明的实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于空时数据的高分辨目标方位估计方法，其特征在于，所述包括：

步骤一，通过传感器阵列对信号进行采样，获取到第一信号矢量x₁(t)、第二信号矢量x₂(t)，其中x₁(t)＝A₁s(t)+n₁(t)，x₂(t)＝A₁Φs(t)+n₂(t)；

步骤二，在开始采样后的第一个时刻，构造与所述第一信号矢量、所述第二信号矢量的对应的矩阵

$R_{11}^{\left(1\right)}=A_1{R}_s\left(1\right)A_1^H+{\mathrm{\σ}}_1^2\left(1\right)I,$

$R_{22}^{\left(1\right)}=A_1\mathrm{\Φ}{R}_s\left(1\right){\mathrm{\Φ}}^H{A}_1^H+{\mathrm{\σ}}_2^2\left(1\right)I,$

其中R_s为自相关矩阵， 为所述采样过程中的噪声，在开始采样后的第二个时刻，构造与所述第一信号矢量、所述第二信号矢量对应的矩阵

$R_{11}^{\left(2\right)}=A_1{R}_s\left(2\right)A_1^H+{\mathrm{\σ}}_1^2\left(2\right)I,$

$R_{22}^{\left(2\right)}=A_1\mathrm{\Φ}{R}_s\left(2\right){\mathrm{\Φ}}^H{A}_1^H+{\mathrm{\σ}}_2^2\left(2\right)I,$

依次在每个时刻均构造两个相关矩阵，并根据时间顺序，将全部相关矩阵组合为相关矩阵组其中k代表时刻数，Λ_m＝Φ^m-1， $B_m^{\left(k\right)}=A_m{R}_s^H\left(k\right);$

步骤三，将所述全部相关矩阵进行特征分解，对所述全部相关矩阵中的噪声进行估计，得到每个噪声的估计值 对所述全部特征矩阵进行去噪，得到去噪相关矩阵组

步骤四，构造初始矩阵U(0)，建立第一代价函数

$\tilde{C}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{m=1}{\overset{2}{\mathrm{\Σ}}}\left\{{\left[R_{\mathrm{mm}}^{\left(k\right)}\right]}^{H}U(l-1)U^H(l-1)R_{\mathrm{mm}}^{\left(k\right)}\right\},$

其中，l代表循环次数，对所述第一代价函数进行特征值分解，得到第一分解式

${\tilde{C}}_{\mathrm{eig}}=\stackrel{\‾}{V}(l-1)\stackrel{\‾}{D}(l-1){\stackrel{\‾}{V}}^H(l-1),$

其中，为特征值矩阵,为特征向量矩阵,令矩阵V(l-1)为矩阵的前P个大特征值对应的特征向量形成的矩阵,建立第二代价函数

其中，l代表循环次数，对所述第二代价函数进行特征值分解，得到第二分解式

其中，为特征值矩阵,为特征向量矩阵,所述第一代价函数、第二代价函数为U(l-1)到U(l)的一次循环，对所述U(l)循环计算，直至满足

${\left|\right|U\left(l\right)U^H\left(l\right)-U(l-1)U^H(l-1)\left|\right|}_F<\mathrm{\&epsiv;}$

时停止循环，获取停止循环时的l；

步骤五，根据所述停止循环时的l，确定信号子空间U(l)，令得到投影矩阵根据所述投影矩阵得到角度估计

$P\left(\widehat{\mathrm{\&theta;}}\right)=\frac{1}{\left|a^H\left(\widehat{\mathrm{\&theta;}}\right)\right[I-{\tilde{U}}_s{\tilde{U}}_s^H\left]a\left(\widehat{\mathrm{\&theta;}}\right)\right|}.$