CN108120967B - 一种平面阵列doa估计方法及设备 - Google Patents

Abstract

本发明公开了一种平面阵列DOA估计方法及设备。所述平面阵列DOA估计方法包括：步骤1：使用平面阵列接收来自至少一个信源的信号，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对；步骤2：基于所述互质线性子阵列对的实际接收信号，计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵、以及两个线性子阵列的实际接收信号之间的互协方差矩阵；步骤3：基于所述自协方差矩阵和互协方差矩阵估计所述至少一个信源的DOA的数值。解决了二维空间谱搜索及角度配对等导致算法复杂度高的问题，能够以较低的复杂度实现高精度的波达方向估计，适于实时性要求高的场合。

Description

一种平面阵列DOA估计方法及设备

技术领域

本发明涉及通信信号处理领域，尤其涉及一种平面阵列DOA估计方法及设备。

背景技术

波达方向(DOA)估计是雷达、声呐及无线通信等诸多领域的重要研究内容。利用多天线构成阵列可用来估计信号的来波方向。然而，传统的天线阵列通常受到相邻阵元间隔不大于半波长的制约，导致天线数一定时，阵列孔径受到限制。目前，一种称为互质阵列的天线阵列结构突破了半波长的限制而受到了广泛关注，可以有效提升角度的估计性能。但是目前的研究大多集中于一维线性阵列，无法满足二维DOA估计的需要。

另外，在传统二维DOA估计研究中，估计性能受到阵列孔径的限制，在尺寸受限时布设大孔径阵列天线非常困难，导致估计性能不好。同时二维空间谱的完全搜索和二维角度配对将导致计算复杂度过高，因此，难以应用至实时性要求高的场合。

因此，至少需要提出新的技术方案来对现有技术方案的不足之处进行改进。

发明内容

本发明的目的是通过以下技术方案实现的。

根据本发明的平面阵列DOA估计方法，包括：

步骤1：使用平面阵列接收来自至少一个信源的信号，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对；

步骤2：基于所述互质线性子阵列对的实际接收信号，计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵、以及两个线性子阵列的实际接收信号之间的互协方差矩阵；

步骤3：基于所述自协方差矩阵和互协方差矩阵估计所述至少一个信源的DOA的数值。

根据本发明的平面阵列DOA估计方法，所述相互平行的互质线性子阵列对包括：

第一线性子阵列和第二线性子阵列，所述第一线性子阵列包含M₁个阵元，所述M₁个阵元沿y轴方向以M₂λ/2为间隔进行布置，所述第二线性子阵列包含M₂个阵元，所述M₂个阵元沿y轴方向以M₁λ/2为间隔进行布置，所述第一线性子阵列的第一个阵元和所述第二线性子阵列的第一个阵元沿x轴方向对齐且间隔小于或等于λ/2，其中，M₁和M₂是互质的正整数，λ为信号波长。

根据本发明的平面阵列DOA估计方法，所述步骤2包括：

步骤2-1：获取所述互质线性子阵列对的实际接收信号x₁(t)和x₂(t)，其中，

x₁(t)和x₂(t)分别表示第一线性子阵列和第二线性子阵列在t时刻的实际接收信号，矩阵A₁＝[a₁(α₁),…,a₁(α_K)]和A₂＝[a₂(α₁),…,a₂(α_K)]分别表示第一线性子阵列和第二线性子阵列沿y轴的流型矩阵，

和

分别表示第一线性子阵列和第二线性子阵列对于第k个信源的导向矢量，矩阵

为对角矩阵，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为由入射角度分别为(α_k,β_k),k＝1,2,…,K,(K≥1)的K个信源信号组成的信号矢量，上标T表示转置运算，α_k表示入射方向与y轴之间的夹角，β_k表示入射方向与x轴之间的夹角，矢量n₁(t)和n₂(t)分别表示第一线性子阵列和第二线性子阵列所实际接收到的均值为零、方差为

的加性高斯白噪声矢量，矢量n₁(t)和n₂(t)与信号矢量s(t)不相关；

步骤2-2：计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵R₁和R₂、以及两个线性子阵列的实际接收信号之间的互协方差矩阵R₁₂：

其中

其中，H表示共轭转置运算，

和

分别表示M₁和M₂阶单位阵，Φ和R_s都是对角阵，

表示对角阵。

根据本发明的平面阵列DOA估计方法，所述步骤3包括：

步骤3-1：对所述互协方差矩阵R₁₂的表达式进行矢量化处理；

步骤3-2：基于所述自协方差矩阵R₁和R₂、以及经过矢量化处理的所述互协方差矩阵R₁₂的表达式，使用一维搜索来估计所述至少一个信源的DOA的数值。

根据本发明的平面阵列DOA估计方法，所述步骤3-1包括：通过以下公式对所述互协方差矩阵R₁₂的表达式进行矢量化处理，

其中，

其中，⊙表示Khatri-Rao积，*表示共轭运算，u表示矩阵ΦR_s对角线上的元素所构成的向量。

根据本发明的平面阵列DOA估计方法，所述步骤3-2包括：

步骤3-2-1：通过以下公式对所述自协方差矩阵R₁和R₂进行特征值分解，求取第k个信源与平行互质阵列之间的夹角α_k的估计值

其中，i＝1,2,，

和

分别表示自协方差矩阵经过特征值分解后的信号部分和噪声部分，表示第k个信源与两个线性子阵列各自的参考阵元之间的夹角a_i,k的估计值，

表示角度搜索区域，Δ_min,i表示最小搜索间隔；

步骤3-2-2：通过以下公式求取第k个信源与平行互质阵列之间的夹角β_k的估计值

其中，

其中，

分别表示矩阵A、A₁、A₂的估计，

表示矩阵

的伪逆。

根据本发明的平面阵列DOA估计设备，所述设备包括平面阵列、处理器和存储有可执行指令的存储器，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对，所述处理器执行所述可执行指令来完成根据上文所述的方法中的步骤。

根据本发明的平面阵列DOA估计设备，包括：

平面阵列模块，用于接收来自至少一个信源的信号，所述平面阵列模块包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对；

协方差矩阵计算模块，其与平面阵列模块连接，用于基于所述互质线性子阵列对的实际接收信号，计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵、以及两个线性子阵列的实际接收信号之间的互协方差矩阵；

DOA估计模块，其与协方差矩阵计算模块连接，用于基于所述自协方差矩阵和互协方差矩阵估计所述至少一个信源的DOA的数值。

本发明的优点在于：解决了二维空间谱搜索及角度配对等导致算法复杂度高的问题，能够以较低的复杂度实现高精度的波达方向估计，适于实时性要求高的场合。

附图说明

通过阅读下文具体实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出具体实施方式的目的，而并不认为是对本发明的限制。而且在整个附图中，用相同的参考符号表示相同的部件。在附图中：

图1示出了根据本发明实施方式的平面阵列DOA估计方法的示意流程图。

图2示出了根据本发明实施方式的相互平行的互质线性子阵列的结构示意图。

图3示出了包括根据本发明实施方式的平面阵列DOA估计方法在内的几种不同估计方法的均方根误差与信噪比之间的关系示意图。

图4示出了包括根据本发明实施方式的平面阵列DOA估计方法在内的几种不同估计方法的均方根误差与快拍数之间的关系示意图。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施方式。虽然附图中显示了本公开的示例性实施方式，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。

图1示出了根据本发明实施方式的平面阵列DOA估计方法100的示意流程图。

如图1所示，平面阵列DOA估计方法100包括以下步骤：

步骤S102：使用平面阵列接收来自至少一个信源的信号，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对。

步骤S104：基于所述互质线性子阵列对的实际接收信号，计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵、以及两个线性子阵列的实际接收信号之间的互协方差矩阵。

步骤S106：基于所述自协方差矩阵和互协方差矩阵估计所述至少一个信源的DOA的数值。

图2示出了根据本发明实施方式的相互平行的互质线性子阵列200的结构示意图。

如图2所示，平面阵列DOA估计方法100所使用的相互平行的互质线性子阵列对200包括第一线性子阵列(即，图2所示的子阵列1)和第二线性子阵列(即，图2所示的子阵列2)。

第一线性子阵列包含M₁个阵元，所述M₁个阵元沿y轴方向以M₂λ//2为间隔进行布置，所述第二线性子阵列包含M₂个阵元，所述M₂个阵元沿y轴方向以M₁λ/2为间隔进行布置，所述第一线性子阵列的第一个阵元和所述第二线性子阵列的第一个阵元沿x轴方向对齐且间隔小于或等于λ/2，其中，M₁和M₂是互质的正整数，λ为信号波长。

尽管在图2中仅仅示出了一个相互平行的互质线性子阵列对200，然而，平面阵列DOA估计方法100所使用的平面阵列可以包括多个相互平行的互质线性子阵列对200。

尽管在图1中未示出，然而，可选地，上述步骤S104可以包括以下步骤：

步骤2-1：获取所述互质线性子阵列对的实际接收信号x₁(t)和x₂(t)。

即，(1)构建二维接收阵列模型。

即，利用如图2所示的二维接收阵列模型接收来自至少一个信源的信号。

其中，

x₁(t)和x₂(t)分别表示第一线性子阵列和第二线性子阵列在t时刻的实际接收信号，矩阵A₁＝[a₁(α₁),…,a₁(α_K)]和A₂＝[a₂(α₁),…,a₂(α_K)]分别表示第一线性子阵列和第二线性子阵列沿y轴的流型矩阵，和分别表示第一线性子阵列和第二线性子阵列对于第k个信源的导向矢量，矩阵

为对角矩阵，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为由入射角度分别为(α_k,β_k),k＝1,2,…,K,(K≥1)的K个信源信号组成的信号矢量，上标T表示转置运算，α_k表示入射方向与y轴之间的夹角，β_k表示入射方向与x轴之间的夹角，矢量n₁(t)和n₂(t)分别表示第一线性子阵列和第二线性子阵列所实际接收到的均值为零、方差为

的加性高斯白噪声矢量，矢量n₁(t)和n₂(t)与信号矢量s(t)不相关。

步骤2-2：计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵R₁和R₂、以及两个线性子阵列的实际接收信号之间的互协方差矩阵R₁₂。

即，(2)二维估计问题降维中的互协方差矩阵R₁₂的计算、以及(3)最小一维空间谱搜索中的自协方差矩阵R₁和R₂的计算。

根据x₁(t)和x₂(t)的定义，互相关矩阵表示为

其中(·)^H表示共轭转置运算，

表示对角阵。由于Φ和R_s都是对角矩阵，ΦR_s亦是对角阵。

两个子阵对应的协方差矩阵分别为

和

其中，

和

分别表示M₁和M₂阶单位阵。

尽管在图1中未示出，然而，可选地，上述步骤S106可以包括以下步骤：

步骤3-1：对所述互协方差矩阵R₁₂的表达式进行矢量化处理。

可选地，所述步骤3-1包括：通过以下公式对所述互协方差矩阵R₁₂的表达式进行矢量化处理，

即，(2)二维估计问题降维中的具体降维方法。

其中，

其中，⊙表示Khatri-Rao积，*表示共轭运算，u表示矩阵ΦR_s对角线上的元素所构成的向量。

矩阵A仅与角度α_k,k＝1,2,…,K相关。同时R_s为实值对角矩阵，而Φ中对角元素与β_k,k＝1,2,…,K相关，所以向量u中元素与真实相位β_k相关。因此可实现二维角度的降维。

步骤3-2：基于所述自协方差矩阵R₁和R₂、以及经过矢量化处理的所述互协方差矩阵R₁₂的表达式，使用一维搜索来估计所述至少一个信源的DOA的数值。

可选地，所述步骤3-2包括：

步骤3-2-1：通过以下公式对所述自协方差矩阵R₁和R₂进行特征值分解，求取第k个信源与平行互质阵列之间的夹角α_k的估计值

即，(3)最小一维空间谱搜索中的一维空间谱搜索方法。

由于R₁和R₂仅与角度α_k,k＝1,2,…,K相关。因此，可将经典的多重信号分类算法应用至每个子阵来获得角度估计。具体为，协方差矩阵R_i(i＝1,2)的特征值分解为：

其中U_i,s和U_i,n分别表示信号子空间和噪声子空间。对每个子阵而言，DOA可通过最大化

来获得。

由于两个子阵均为稀疏阵，对每个真实相位而言都存在多个模糊相位。对第i个子阵而言，真实相位和模糊相位之间的关系为：

其中P_i为整数，此关系可用来减少空间谱的搜索区域，如定理1所示。

定理1：对子阵i而言，其阵元间隔为

最小频谱搜索间隔Δ_min,i为

搜索区域为

证明：令γ表示搜索区域的起点，Δ表示搜索间隔。为使Δ最小化，可建立如下优化问题

γ∈[0,π),Δ>0,

利用拉格朗日乘子算法，可得

其中λ≥0表示拉格朗日乘子。利用Karush-Kuhn-Tucker(KKT)条件，可得：

通过对上式求解，可得最优解为

和

根据定理1可知，一维搜索区域得到了最小化，从而可极大降低实现复杂度。通过最小区域搜索，可获得任意相关的相位，其他相关相位可通过真实相位与模糊相位之间的关系计算可得。最终利用两个子阵间的互质关系，即两个子阵只有在真实相位处才能获得公共谱峰，通过寻找两个子阵相位的公共部分，即

即可确定真实角度α_k,k＝1,2,…,K的估计值

步骤3-2-2：通过以下公式求取第k个信源与平行互质阵列之间的夹角β_k的估计值

即，(4)角度自动配对。

如上所述，可通过最小一维搜索获得角度α_k,k＝1,2,…,K的估计值

利用

可得两个子阵对应的阵列流型矩阵分别为

和

因此，r₁₂可以表示为：

其中

利用最小二乘准则，上述问题可解为：

其中

表示矩阵

的伪逆。根据u的定义可知，其相位揭示了相位β_k,k＝1,2,…,K的值，因此有：

注意到，β_k的估计依赖于α_k的估计值，因此估计过程中实现了角度α_k和β_k实现了自动配对，降低了实现复杂度。

DOA估计算法的复杂度主要由三部分构成，即协方差矩阵计算，特征值分解和空间谱搜索。由于空间谱搜索产生的复杂度远大于其他两部分，因此不同算法的复杂度可由空间谱搜索的复杂度表征，如表1所示，其中Δ_α和Δ_β分别表示两个角度的搜索间隔，K为信号源个数。可以看出，由于在空间谱搜索的估计算法中采用了合理的搜索区域、以及角度α_k和β_k之间的自动配对，本发明所提出的上述DOA估计方法100的计算复杂度远远低于其他两种算方法，这在M₁或M₂较大的场合变得更加明显。

表1

为了使本领域技术人员更直观地了解本文提出的平面阵列DOA估计方法100的技术效果，下文中给出了采用平面阵列DOA估计方法100所得到的部分结果。

图3示出了包括根据本发明实施方式的平面阵列DOA估计方法在内的几种不同估计方法的均方根误差与信噪比之间的关系示意图。

图3为快拍数为400时，本发明方法(即，图3中的“所提算法”)和完全空间谱搜索方法、MUSIC方法的均方根误差及CRB界随信噪比的变化情况。从图3可以看出，本发明方法的性能与完全空间谱搜索方法的性能基本一致，但是其计算复杂度却有极大的降低。

图4示出了包括根据本发明实施方式的平面阵列DOA估计方法在内的几种不同估计方法的均方根误差与快拍数之间的关系示意图。

图4为信噪比为5分贝时，本发明方法(即，图4中的“所提算法”)和完全空间谱搜索方法、MUSIC方法的均方根误差及CRB界随快拍数的变化情况。可以看出，本发明方法的性能与完全空间谱搜索方法的性能基本一致，但是其复杂度更低。

以上结果表明，根据本发明的上述技术方案通过包括(1)构建二维接收阵列模型；(2)二维估计问题降维；(3)最小一维空间谱搜索；(4)角度自动配对等步骤。能够在保证估计精度的同时，可极大降低波达方向估计的实现复杂度。针对二维空间谱搜索及角度配对等导致复杂度高的问题，本发明利用平行子阵之间的互质特性及相关性来实现降维与配对。本发明以较低的复杂度实现高精度的波达方向估计，适于实时性要求高的场合。

结合上文提出的平面阵列DOA估计方法100，还提出了一种平面阵列DOA估计设备，所述设备包括平面阵列、处理器和存储有可执行指令的存储器，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对，所述处理器执行所述可执行指令来完成根据上文所述的平面阵列DOA估计方法100中的步骤。

结合上文提出的平面阵列DOA估计方法100，还提出了另一种平面阵列DOA估计设备，包括：

平面阵列模块，用于接收来自至少一个信源的信号，所述平面阵列模块包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对；

协方差矩阵计算模块，其与平面阵列模块连接，用于基于所述互质线性子阵列对的实际接收信号，计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵、以及两个线性子阵列的实际接收信号之间的互协方差矩阵；

DOA估计模块，其与协方差矩阵计算模块连接，用于基于所述自协方差矩阵和互协方差矩阵估计所述至少一个信源的DOA的数值。

以上所述，仅为本发明示例性的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (5)

1.一种平面阵列DOA估计方法，其特征在于，包括：

步骤1：使用平面阵列接收来自至少一个信源的信号，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对；

步骤2：基于所述互质线性子阵列对的实际接收信号，计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵、以及两个线性子阵列的实际接收信号之间的互协方差矩阵；

步骤3：基于所述自协方差矩阵和互协方差矩阵估计所述至少一个信源的DOA的数值；所述步骤3包括：

步骤3-1：对所述互协方差矩阵R₁₂的表达式进行矢量化处理；

步骤3-2：基于所述自协方差矩阵R₁和R₂、以及经过矢量化处理的所述互协方差矩阵R₁₂的表达式，使用一维搜索来估计所述至少一个信源的DOA的数值；

所述步骤3-1包括：通过以下公式对所述互协方差矩阵R₁₂的表达式进行矢量化处理，

其中，

其中，⊙表示Khatri-Rao积，*表示共轭运算，u表示矩阵ΦR_s对角线上的元素所构成的向量；

所述步骤3-2包括：

步骤3-2-1：通过以下公式对所述自协方差矩阵R₁和R₂进行特征值分解，求取第k个信源与平行互质阵列之间的夹角α_k的估计值

其中，i＝1,2,

和

分别表示自协方差矩阵经过特征值分解后的信号部分和噪声部分，

表示第k个信源与两个线性子阵列各自的参考阵元之间的夹角a_i,k的估计值，

表示角度搜索区域，Δ_min,i表示最小搜索间隔；

步骤3-2-2：通过以下公式求取第k个信源与平行互质阵列之间的夹角β_k的估计值

其中，

其中，

分别表示矩阵A、A₁、A₂的估计，

表示矩阵

的伪逆。

2.根据权利要求1所述的平面阵列DOA估计方法，其特征在于，所述相互平行的互质线性子阵列对包括：

第一线性子阵列和第二线性子阵列，所述第一线性子阵列包含M₁个阵元，所述M₁个阵元沿y轴方向以M₂λ/2为间隔进行布置，所述第二线性子阵列包含M₂个阵元，所述M₂个阵元沿y轴方向以M₁λ/2为间隔进行布置，所述第一线性子阵列的第一个阵元和所述第二线性子阵列的第一个阵元沿x轴方向对齐且间隔小于或等于λ/2，其中，M₁和M₂是互质的正整数，λ为信号波长。

3.根据权利要求2所述的平面阵列DOA估计方法，其特征在于，所述步骤2 包括：

步骤2-1：获取所述互质线性子阵列对的实际接收信号x₁(t)和x₂(t)，其中，

x₁(t)和x₂(t)分别表示第一线性子阵列和第二线性子阵列在t时刻的实际接收信号，矩阵A₁＝[a₁(α₁),…,a₁(α_K)]和A₂＝[a₂(α₁),…,a₂(α_K)]分别表示第一线性子阵列和第二线性子阵列沿y轴的流型矩阵，

和

分别表示第一线性子阵列和第二线性子阵列对于第k个信源的导向矢量，矩阵

为对角矩阵，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为由入射角度分别为(α_k,β_k),k＝1,2,…,K,(K≥1)的K个信源信号组成的信号矢量，上标T表示转置运算，α_k表示入射方向与y轴之间的夹角，β_k表示入射方向与x轴之间的夹角，矢量n₁(t)和n₂(t)分别表示第一线性子阵列和第二线性子阵列所实际接收到的均值为零、方差为

的加性高斯白噪声矢量，矢量n₁(t)和n₂(t)与信号矢量s(t)不相关；

步骤2-2：计算所述互质线性子阵列对的两个线性子阵列各自的实际接收信号的自协方差矩阵R₁和R₂、以及两个线性子阵列的实际接收信号之间的互协方差矩阵R₁₂：

其中，

其中，H表示共轭转置运算，

和

分别是M₁和M₂阶单位阵，Φ和R_s都是对角阵，

表示对角阵。

4.根据权利要求1所述的平面阵列DOA估计方法，其特征在于，所述步骤3-1包括：通过以下公式对所述互协方差矩阵R₁₂的表达式进行矢量化处理，

其中，

其中，⊙表示Khatri-Rao积，*表示共轭运算，u表示矩阵ΦR_s对角线上的元素所构成的向量。

5.一种平面阵列DOA估计设备，所述设备包括平面阵列、处理器和存储有可执行指令的存储器，其特征在于，所述平面阵列包括至少一个平行互质阵列，所述平行互质阵列包括相互平行的互质线性子阵列对，所述处理器执行所述可执行指令来完成根据权利要求1至4中任一项所述的方法中的步骤。

Images