CN109507634A - 一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法 - Google Patents

Abstract

本发明公开了一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，可应用于雷达、声呐和移动通信领域。该方法包括通过接收信号的协方差或利用接收信号直接计算传播算子；利用传播算子估计信号子空间；基于秩亏准则构建一维空间谱函数，通过谱峰搜索或者多项式求根实现对DOA的估计。该发明首先利用线性的传播算子估计信号子空间，继而基于秩亏准则提出用于波达方向估计的一维空间谱函数。同时，该发明提供了一种当阵列构型满足特殊条件时的基于多项式求根的快速波达方向估计方法，降低了计算复杂度。运用该算法可以获得较传统算法更低的计算复杂度，同时估计精度保持接近。

Description

一种任意传感器阵列下的基于传感算子的盲远场信号波达方 向估计方法

技术领域

本发明属于信号处理领域，特别是涉及一种任意传感器阵列构型下基于秩亏准则的信号波达方向估计方法。

背景技术

利用传感器阵列实现盲信号波达方向(Direction-of-arrival,DOA)估计是雷达、声呐、无线通信、射电天文等领域的一项基本课题。DOA估计领域的经典算法有常规波束成形算法、最大似然算法、最小二乘算法、子空间算法和张量分解算法等。其中，以多重信号分类(Multiple Signal Classification,MUSIC)算法、基于旋转不变性的信号参数估计算法(Estimation of Signal Parameters via Rotational Invariance Techniques，ESPRIT)为代表的子空间算法由于其高分辨率、高精度、相对较低的复杂度等优点，成为学术界和工程界关注的焦点。MUSIC算法和ESPRIT算法分别代表了子空间算法的两种思路。MUSIC算法利用噪声子空间与阵列流型之间的正交性，通过对空间谱函数的搜索得到对DOA的估计，优点是适用于任意构型阵列，缺点是谱峰搜索过程复杂度较高。ESPRIT算法利用信号子空间和旋转不变阵列的性质，优点是对DOA估计具有闭式解，复杂度较低，缺点是只适用于旋转不变阵列。基于这两种算法，相应的改进算法层出不穷。例如针对MUSIC改进的求根MUSIC算法降低了复杂度，针对ESPRIT改进的广义ESPRIT算法利用秩亏准则解除了旋转不变阵列的限制。然而，这些算法均采用特征值分解运算对子空间进行估计，这种运算的计算复杂度与阵元数的立方成正比，当阵元数较大时，算法复杂度较高。

传播算子(Propagator Method,PM)算法同样利用子空间理论进行DOA估计，但其是利用传播算子这一可通过对接收信号简单矩阵运算得到的线性算子得到对子空间的估计，避免了复杂度较高的特征值分解运算。基于传播算子对MUSIC和ESPRIT算法的改进分别有N-PM算法和S-PM算法，相关的改进算法可用于实现电磁矢量传感器、稀疏L型阵列、双基地多输入多输出雷达等装置下的DOA估计。

发明内容

发明目的：提供一种任意传感器阵列构型下基于秩亏准则的信号波达方向估计方法，该方法可以降低传统广义ESPRIT算法的复杂度，更快速的获得估计结果，减少对硬件计算能力的依赖和对存储空间的占用，同时，估计精度与传统广义ESPRIT算法保持接近。

技术方案：为实现上述发明目的，本发明采用以下技术方案：

一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，包括以下步骤：

(1)通过接收信号的协方差或利用接收信号直接计算传播算子；

(2)利用传播算子估计信号子空间；

(3)基于秩亏准则构建一维空间谱函数，通过谱峰搜索或者多项式求根实现对DOA的估计。

进一步的，步骤(1)中基于信号独立假设，阵列流型矩阵A为满秩矩阵，其K行线性无关，剩余(M-K)行是前K行的线性组合，则将阵列流型矩阵A按如下形式分块：

其中，G₁是一个包含A前K行的非奇异K×K矩阵，G₂是包含A后(M-K)行的 (M-K)×K维矩阵，P^HG₁＝G₂，矩阵P即为传播算子；

将T组快拍采样数据记为Y＝[y(1),…,y(T)]，并计算采样协方差矩阵对接收数据Y和采样协方差矩阵按如下方式分块：

其中，分别包含Y的前K行和后(M-K)行，分别包含的左边K列和后边(M-K)列；对传播算子的估计由下列两式分别获得：

进一步的，步骤(2)中定义一个M×K维的新矩阵Q＝[I_K,P]^H，其中I_K是K×K维的单位矩阵，由传播算子P的定义，易得：

其中，A为阵列流型矩阵，G₁为一个包含A前K行的非奇异K×K矩阵，因为矩阵G₁非奇异，等式成立，且可得Q与阵列流型矩阵A张成相同的列空间，由此，利用传播算子实现对信号子空间的估计。

进一步的，步骤(3)中将阵列分成两阵元数均为(M-1)的子阵，假设两个子阵分别包含前(M-1)和后(M-1)个阵元，两子阵相应的选择矩阵定义为J₁和J₂，且其中I_M-1是(M-1)×(M-1)的单位矩阵，o是(M-1)×1的零矢量；

两子阵的阵列流型矩阵分别表示为A₁＝J₁A,A₂＝J₂A，令a₁(θ_k)、a₂(θ_k)分别表示 A₁和A₂的第k列，其表达式为：

阵列流型矢量a₁(θ_k)、a₂(θ_k)之间存在如下转移关系：a₂(θ_k)＝Γ(θ_k)a₁(θ_k)，其中δ_i＝c_i+1-c_i,i＝1,...,M-1；

两子阵接收信号的信号子空间分别表示为：

定义w(θ)＝[cosθ,sinθ]^T并构建一个(M-1)×(M-1)的对角矩阵和矩阵W(θ)＝Q₂-C(θ)Q₁，由传播算子定义和子阵转移关系得：

其中，H(θ)＝[(Γ(θ₁)-C(θ))a₁(θ₁),…,(Γ(θ_K)-C(θ))a₁(θ_K)]，

因为G₁非奇异，当θ＝θ_k时，H(θ)的第k列(Γ(θ_k)-C(θ))a₁(θ_k)的所有元素将为零；此时，如果K≤M-1，矩阵W(θ)将产生秩亏，方阵EW(θ)的行列式值将为零，其中E 是一个任意的(M-1)×K的满秩矩阵，E＝W^H(θ)；

由此，得到以下一维空间谱函数：

更进一步的，令δ_min表示集合{δ₁,…,δ_M-1}中2范数最小的元素，当阵列构型不满足δ_i＝α_iδ_min,时，则，通过对一维空间谱函数f(θ)估计DOA，在[-π/2,π/2]范围内对DOA进行搜索，K个最大的峰值即是DOA的估计值。

更进一步的，令δ_min表示集合{δ₁,…,δ_M-1}中2范数最小的元素，当阵列构型满足δ_i＝α_iδ_min,时，谱峰搜索问题转换为多项式求根问题：假设δ₁＝δ_min且α₁≤α₂≤…≤α_M-1，定义则C(θ)和W(θ)分别表示为：

W(z)＝Q₂-C(z)Q₁；

同时W^H(θ)被表示为空间谱函数的分母被表示为如下的多项式形式：

p(z)＝det{F(z)W(z)}；

因为α_i,i＝1,…,M-1均为正数，DOA将通过对多项式p(z)求根得到，单位圆内使空间谱函数最大的K个根即是DOA估计值。

有益效果：与现有技术相比，本发明首先利用线性的传播算子估计信号子空间，继而基于秩亏准则，提出了用于波达方向估计的一维空间谱函数，实现了对S-PM算法的改进，去除了传统的S-PM算法对于阵列构型的限制，将其推广到任意传感器阵列构型的情形。与传统的广义ESPRIT算法相比，其复杂度更低，而估计精度保持接近。同时，提出了一种当阵列构型满足特定条件时的基于求根多项式的快速波达方向估计方法，进一步降低了计算复杂度；运用该方法可以获得较传统算法更低的计算复杂度，同时估计精度保持接近。

附图说明

图1是本发明下阵列模型图；

图2是本发明下算法流程图；

图3是本发明与广义ESPRIT算法的复杂度对比图；

图4是本发明估计均方误差值与N-PM、广义ESPRIT算法和克拉美罗界对比图。

具体实施方式

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

以一阵元任意分布的传感器阵列为例，如图1所示，其阵元数为M，阵元坐标表示为(x_m,y_m),m＝1,…,M，且相邻阵元之间间隔不大于入射信号波长的一半。空间中存在K个远场、窄带、独立信源，其信号中心频率均为ω₀，且有M>K。假设信源与阵列之间的距离足够远，信号波前可被认为平面波。假设阵元相互独立，无互耦等干扰因素影响，入射信号可被表示成：

y(t)＝As(t)+n(t),t＝1,…,T(1)；

其中A表示阵列流型矩阵，θ＝[θ₁,…,θ_K]表示信号方向与y轴之间的夹角，且有θ_k∈[-π/2,π/2]，s(t)＝[s₁(t),…,s_K(t)]^T为K×1的复信号矢量，n(t)表示平稳零均值的加性高斯白噪声，方差为t表示采样时间点，T表示快拍数。阵列流型矩阵A的表达式为A＝[a(θ₁),…,a(θ_K)]，其中为对应于第k个信号的阵列流型矢量，c_m＝[x_m,y_m]^T,w(θ_k)＝[cosθ_k,sinθ_k]^T，λ＝c/ω₀为波长，c为光速。

如图2所示，一种任意传感器阵列构型下基于秩亏准则的信号波达方向估计方法，首先利用线性的传播算子估计信号子空间，继而基于秩亏准则提出用于波达方向估计的一维空间谱函数。同时，该发明提供了一种在特定阵列构型下的基于多项式求根的快速波达方向估计算法，替代了复杂度较高的谱峰搜索过程。运用该算法可以获得较传统算法更低的计算复杂度，同时估计精度保持接近。具体包括以下步骤：

(1)通过接收信号的协方差或利用接收信号直接计算传播算子；

基于信号独立假设，阵列流型矩阵A为满秩矩阵，其K行线性无关，而剩余(M-K) 行是前K行的线性组合。此时，将阵列流型矩阵A按如下形式分块：

其中G₁是一个包含A前K行的非奇异K×K矩阵，G₂是包含A后(M-K)行的(M-K)×K维矩阵。矩阵G₂可表示为G₁的线性变换：P^HG₁＝G₂，矩阵P即为传播算子。

传播算子P可简单得从接收信号或者采样协方差矩阵中计算得到。将T组快拍采样数据记为Y＝[y(1),…,y(T)]，并计算采样协方差矩阵对接收数据Y和采样协方差矩阵按如下方式分块：

其中分别包含Y的前K行和后(M-K)行，分别包含的左边K列和后边(M-K)列。对传播算子的估计可由下列两式分别获得：

(2)利用传播算子估计信号子空间；

定义一个M×K维的新矩阵Q＝[I_K,P]^H，其中I_K是K×K维的单位矩阵。由传播算子P的定义，易得：

因为矩阵G₁非奇异，等式成立，且可得Q与阵列流型A张成相同的列空间，也即span{Q}＝span{A}。由此，利用传播算子实现对信号子空间的估计。

(3)基于秩亏准则构建一维空间谱函数，通过谱峰搜索或者多项式求根实现对DOA的估计；

将阵列分成两阵元数均为(M-1)的子阵，在不失一般性的基础上，假设两个子阵分别包含前(M-1)和后(M-1)个阵元。两子阵相应的选择矩阵定义为J₁和J₂，且其中I_M-1是(M-1)×(M-1)的单位矩阵，o是(M-1)×1的零矢量。

两子阵的阵列流型矩阵可分别表示A₁＝J₁A,A₂＝J₂A。令a₁(θ_k),a₂(θ_k)分别表示 A₁和A₂的第k列，也即两个子阵对应于第k个信号的阵列流型矢量，其表达式为：

阵列流型矢量a₁(θ_k),a₂(θ_k)之间存在如下转移关系：a₂(θ_k)＝Γ(θ_k)a₁(θ_k)，其中δ_i＝c_i+1-c_i,i＝1,...,M-1

两子阵接收信号的信号子空间可分别表示为：

定义w(θ)＝[cosθ,sinθ]^T并构建一个(M-1)×(M-1)的对角矩阵和矩阵W(θ)＝Q₂-C(θ)Q₁。由传播算子定义和子阵转移关系可得：

其中H(θ)＝[(Γ(θ₁)-C(θ))a₁(θ₁),…,(Γ(θ_K)-C(θ))a₁(θ_K)]。

因为G₁非奇异，当θ＝θ_k时，H(θ)的第k列(Γ(θ_k)-C(θ))a₁(θ_k)的所有元素将为零。此时，如果K≤M-1，矩阵W(θ)将产生秩亏，方阵EW(θ)的行列式值将为零，其中E 是一个任意的(M-1)×K的满秩矩阵。本发明中，选择E＝W^H(θ)。

由此，可通过对如下一维空间谱函数估计DOA：

在[-π/2,π/2]范围内对DOA进行搜索，K个最大的峰值即是DOA的估计值。

令δ_min表示集合{δ₁,…,δ_M-1}中2范数最小的元素，当阵列构型满足δ_i＝α_iδ_min,时，谱峰搜索问题可转换为多项式求根问题。具体而言，不失一般性，假设δ₁＝δ_min且α₁≤α₂≤…≤α_M-1。定义则C(θ)和W(θ)可分别表示为：

W(z)＝Q₂-C(z)Q₁ (15)；

同时W^H(θ)可被表示为空间谱函数的分母可被表示为如下的多项式形式：

p(z)＝det{F(z)W(z)} (16)；

因为α_i,i＝1,…,M-1均为正数，DOA可通过对多项式p(z)求根得到，单位圆内使空间谱函数最大的K个根即是DOA估计值。

图3是本发明与广义ESPRIT算法的复杂度对比图，图中参数K＝2，T＝150，“Proposed-search”、“Proposed-root”分别表示采用空间谱搜索和多项式求根实现估计的本发明方法，“GESPRIT-search”、“GESPRIT-root”分别表示采用空间谱搜索和多项式求根实现估计的广义ESPRIT算法。如图3所示，无论是采用空间谱搜索还是多项式求根方法，本发明的复杂度在同等条件下均低于广义ESPRIT算法。

图4是本发明估计均方误差值与N-PM、广义ESPRIT算法和克拉美罗界对比图，“Proposed”、“N-PM”、“GESPRIT”、“CRB”分别表示本发明方法、N-PM算法、广义 ESPRIT算法和克拉美罗界。此处考虑M＝8阵元的任意构型阵列，阵元坐标分别为 (0,0)，(0.3λ,0.25λ)，(0.5λ,0)，(0.85λ,0.25λ)，(1.15λ,0)，(1.4λ,0.25λ)，(1.75λ,0)， (2.1λ,0.25λ)；空间中存在K＝3个信源，角度分别是(20°,40°,60°)，快拍数T＝100。如图4所示，相比N-PM算法，本发明方法的估计精度与广义ESPRIT算法接近，更接近克拉美罗界，而本发明方法的计算复杂度相比广义ESPRIT算法较小。

Claims (6)

1.一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，其特征在于，包括以下步骤：

(1)通过接收信号的协方差或利用接收信号直接计算传播算子；

(2)利用传播算子估计信号子空间；

(3)基于秩亏准则构建一维空间谱函数，通过谱峰搜索或者多项式求根实现对DOA的估计。

2.根据权利要求1所述的一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，其特征在于，步骤(1)中基于信号独立假设，阵列流型矩阵A为满秩矩阵，其K行线性无关，剩余(M-K)行是前K行的线性组合，则将阵列流型矩阵A按如下形式分块：

其中，G₁是一个包含A前K行的非奇异K×K矩阵，G₂是包含A后(M-K)行的(M-K)×K维矩阵，P^HG₁＝G₂，矩阵P即为传播算子；

将T组快拍采样数据记为Y＝[y(1),…,y(T)]，并计算采样协方差矩阵对接收数据Y和采样协方差矩阵按如下方式分块：

其中，分别包含Y的前K行和后(M-K)行，分别包含的左边K列和后边(M-K)列；对传播算子的估计由下列两式分别获得：

3.根据权利要求1所述的一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，其特征在于，步骤(2)中定义一个M×K维的新矩阵Q＝[I_K,P]^H，其中I_K是K×K维的单位矩阵，由传播算子P的定义，易得：

其中，A为阵列流型矩阵，G₁为一个包含A前K行的非奇异K×K矩阵，因为矩阵G₁非奇异，等式成立，且可得Q与阵列流型矩阵A张成相同的列空间，由此，利用传播算子实现对信号子空间的估计。

4.根据权利要求1所述的一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，其特征在于，步骤(3)中将阵列分成两阵元数均为(M-1)的子阵，假设两个子阵分别包含前(M-1)和后(M-1)个阵元，两子阵相应的选择矩阵定义为J₁和J₂，且其中I_M-1是(M-1)×(M-1)的单位矩阵，o是(M-1)×1的零矢量；

两子阵的阵列流型矩阵分别表示为A₁＝J₁A,A₂＝J₂A，令a₁(θ_k)、a₂(θ_k)分别表示A₁和A₂的第k列，其表达式为：

阵列流型矢量a₁(θ_k)、a₂(θ_k)之间存在如下转移关系：a₂(θ_k)＝Γ(θ_k)a₁(θ_k)，其中

两子阵接收信号的信号子空间分别表示为：

定义w(θ)＝[cosθ,sinθ]^T并构建一个(M-1)×(M-1)的对角矩阵和矩阵W(θ)＝Q₂-C(θ)Q₁，由传播算子定义和子阵转移关系得：

其中，H(θ)＝[(Γ(θ₁)-C(θ))a₁(θ₁),…,(Γ(θ_K)-C(θ))a₁(θ_K)]，

因为G₁非奇异，当θ＝θ_k时，H(θ)的第k列(Γ(θ_k)-C(θ))a₁(θ_k)的所有元素将为零；此时，如果K≤M-1，矩阵W(θ)将产生秩亏，方阵EW(θ)的行列式值将为零，其中E是一个任意的(M-1)×K的满秩矩阵，E＝W^H(θ)；

由此，得到以下一维空间谱函数：

5.根据权利要求4所述的一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，其特征在于：令δ_min表示集合{δ₁,…,δ_M-1}中2范数最小的元素，当阵列构型不满足时，则，通过对一维空间谱函数f(θ)估计DOA，在[-π/2,π/2]范围内对DOA进行搜索，K个最大的峰值即是DOA的估计值。

6.根据权利要求4所述的一种任意传感器阵列下的基于传感算子的盲远场信号波达方向估计方法，其特征在于：令δ_min表示集合{δ₁,…,δ_M-1}中2范数最小的元素，当阵列构型满足时，谱峰搜索问题转换为多项式求根问题：假设δ₁＝δ_min且α₁≤α₂≤…≤α_M-1，定义则C(θ)和W(θ)分别表示为：

W(z)＝Q₂-C(z)Q₁；

同时W^H(θ)被表示为空间谱函数的分母被表示为如下的多项式形式：

p(z)＝det{F(z)W(z)}；

因为α_i,i＝1,…,M-1均为正数，DOA将通过对多项式p(z)求根得到，单位圆内使空间谱函数最大的K个根即是DOA估计值。