CN107315161B - 基于压缩感知的非圆信号波达方向角估计方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知的非圆信号波达方向角估计方法，主要解决现有技术中阵元利用率低，信号识别数量少的问题，其实现步骤是：1)获取嵌套阵列输出信号，计算协方差矩阵和椭圆协方差矩阵，构造相应的协方差向量，计算这两个向量中元素的维数，构造相应的虚拟阵列向量，构造稀疏重构向量；2)根据稀疏重构向量构造超完备基并定义第一、第二空域稀疏向量，波达方向角范围以及常数矩阵；3)通过稀疏重构方法获得第一空域稀疏向量的最优估计，绘制幅度谱图，并寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求角度值。本发明在阵元数量有限的情况下大大提高了阵列可识别的信源数目，适用于目标侦察和无源定位。

Description

基于压缩感知的非圆信号波达方向角估计方法

技术领域

本发明属于信号处理技术领域，特别涉及一种电磁信号的阵列信号波达方向角估计方法，可用于对飞机、舰船运动目标的侦察与无源定位。

背景技术

信号的波达方向角DOA估计是阵列信号处理领域的一个重要分支，它是指利用天线阵列对空间声学信号、电磁信号进行感应接收，再运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。

在现代通信中，二相相移键控以及M进制幅移键控等非圆信号的应用越来越多，因此有关非圆信号的DOA估计受到了越来越多的关注。目前已有关于利用阵列天线处理非圆信号的一些方法被提出，比较有代表性的是Zhangmeng Liu等人发表的论文“Direction-of-arrival estimation of noncircular signals via sparse representation”(《IEEETransactions on aerospace and electronic systems》,VOL 48,No.3,pp.2690-26982012.07)中公开了一种利用压缩感知求解进行非圆信号DOA估计的方法，该方法是基于均匀阵列的。

另一方面，为了在较少的阵元条件下得到尽量大的角度自由度，检测更多的信源，一些新的阵列结构被提出，比较有代表性的是嵌套阵列以及互质阵列。P Piya等人在其发表的论文“Nested Arrays:A Novel Approach to Array Processing With EnhancedDegrees of Freedom”(《IEEE transactions on signal processing》,VOL 58,NO.8,August 2010)中公开了一种基于嵌套阵列的DOA估计方法，该方法能够使用M+N个阵元，生成2MN+2N-1个虚拟阵元，可检测MN+N-1个信号。该方法具有估计多于阵元数目的信号数的能力，但是，该阵列的讨论都集中在接收信号为圆信号的条件下，对于如何利用该阵列进行非圆信号的处理目前还没有研究。

在实际应用中，对于具有非圆信号，在给定一定数量阵元的环境下，如果不能合理利用这些阵元以及信号的非圆特性，就不能估计足够多的信号，造成侦察和定位资源的浪费。

发明内容

本发明的目的在于针对上述现有技术存在的不足，提出一种基于压缩感知的非圆信号波达方向角估计方法，以在非圆信号环境下，利用嵌套阵列进行信号处理算法设计，避免因不能合理利用阵元和信号特性造成的资源浪费。

为实现上述目的，本发明技术方案包括如下：

(1)用M+N个天线接收机形成嵌套阵列，其中M、N分别表示两个天线接收阵列的阵元数，其取值范围为M≥1，N≥1；

(2)假设空间中有K个非圆目标信号，通过采样和滤波得到嵌套阵列输出信号：Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]^T，其中，y_i(t)表示嵌套阵列的第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)^T表示矩阵转置运算；

(3)利用嵌套阵列输出信号Y(t)，计算虚拟均匀阵列协方差向量和椭圆协方差向量

(4)根据虚拟均匀阵列协方差向量和椭圆协方差向量构造一个稀疏重构向量r：

(5)构造超完备基A_ds(θ)并定义第一空域稀疏向量u与第二空域稀疏向量δ，波达方向角范围θ以及常数矩阵T：

(5.1)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,将其定义为波达方向角的范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ_q表示第q个刻度值，1≤q≤Q，Q＞＞M；

(5.2)构造一个信号稀疏化后对应的C_d×Q维的差阵列超完备基A_d(θ)和一个C_s×Q维的和阵列超完备基A_s(θ)：

A_d(θ)＝[a_d(θ₁),…,a_d(θ_q),…,a_d(θ_Q)]，

A_s(θ)＝[a_s(θ₁),…,a_s(θ_q),…,a_s(θ_Q)]，

其中，α_d(θ_q)表示第q个刻度值对应的差阵列阵列流型向量，α_s(θ_q)表示第q个刻度值对应的和阵列阵列流型向量，

α_d(θ_q)＝(exp(j2πd·(C_d-1)/2·sinθ_q/λ),.…,1,…exp[-j2πd·(C_d-1)/2·sinθ_q/λ])^T，

α_s(θ_q)＝(exp(-j2πd·2·sinθ_q/λ),……,exp[-j2πd·C_s·sinθ_q/λ])^T，1≤q≤Q，

(5.3)根据差阵列超完备基A_d(θ)与和阵列超完备基A_s(θ)得到超完备基A_ds(θ)：

(5.4)定义一个2Q×1维的空域稀疏向量：u＝[u₁,u₂,…,u_q,…,u_Q,u₁,u₂,…,u_q,…,u_Q]^T，定义一个Q×1维的空域稀疏向量：δ＝[δ₁,δ₂,…,δ_q,…,δ_Q]^T，其中u_q和δ_q分别表示第一未知变量和第二未知变量，1≤q≤Q；

(5.5)定义一个单位矩阵T_Q为Q×Q的对角阵，即对角元素全为1，其它元素全为0，并定义一个常数矩阵T＝[T_Q,T_Q]；

(6)通过稀疏重构方法获得第一空域稀疏向量δ的最优估计

(7)以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以第一空域稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

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

1)本发明采用了嵌套阵列模型进行波达方向角度估计，克服了现有技术中采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目的缺点，提高了在阵元数目相同的条件下的阵列可识别信源数目。

2)本发明利用信号的协方差矩阵R_d，以及椭圆协方差矩阵R_s对信号进行估计，增多了可估计的非圆信号数目。

3)本发明利用了嵌套阵列的特点和非圆信号的特性，将非圆信号在嵌套阵列上进行信号处理，在M+N个阵元上能估计出(MN+M+N-1)/2+MN+N-1个信号，大大提高了阵列的利用率，增加了阵列可识别信源的数目。

4)在波达角度估计中引入压缩感知算法，提高了角度估计的精度以及在较低的信噪比条件下的估计性能。

附图说明

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

图2是本发明中嵌套阵列的结构示意图。

具体实施方式

参照图1，本实例的实现步骤如下：

步骤1：用M+N个天线接收机形成嵌套阵列。

(1a)将每个天线接收机称为一个阵元，用M个天线接收机形成第一均匀线性阵列a，其阵元间距为d，定义第一均匀线性阵列a的第一个阵元为起始阵元，定义起始阵元位置D(1)＝1，第一均匀线性阵列a的其它阵元位置依次为D(2)＝2，D(3)＝3，D(4)＝4，…，D(M)＝M；其中，M的取值范围为M≥1，d的取值范围为0＜d≤λ/2，λ为入射到阵列的窄带信号波长；

(1b)用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+1)d，第二均匀线性阵列b的阵元位置依次设置为D(M+1)＝M+1，D(M+2)＝2(M+1)，D(M+2)＝3(M+1)，…，D(M+N)＝N(M+1)，其中，N的取值范围为N≥1；

(1c)将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为Md的位置，将第二均匀线性阵列b的所有阵元依次放置于第一均匀线性阵列a之后，形成嵌套阵列，如图2所示。

步骤2：获得嵌套阵列输出信号Y(t)。

假设空间中有K个非圆目标信号，使用嵌套阵列天线接收机，对空间目标信号进行快拍采样和匹配滤波操作，得到嵌套阵列输出信号：Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]^T，其中，K的取值范围是K＜MN+M+N-1，y_i(t)表示嵌套阵列第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)^T表示矩阵转置运算。

步骤3：计算协方差矩阵R_d和椭圆协方差矩阵R_s。

利用嵌套阵列输出信号Y(t)，计算协方差矩阵R_d和椭圆协方差矩阵R_s：

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

步骤4：构造等效协方差向量r_d和等效椭圆协方差向量r_s。

根据协方差矩阵R_d和椭圆协方差矩阵R_s中的元素，分别构造等效协方差向量r_d和等效椭圆协方差向量r_s：

其中，R_d(i,j)表示协方差矩阵R_d中位于第i行、第j列的元素，i的取值范围为1≤i≤M+N，j的取值范围为1≤j≤M+N，R_s(i,j)表示椭圆协方差矩阵R_s中位于第i行，第j列的元素。

步骤5：计算等效协方差向量和等效椭圆协方差向量中所有元素的维数。

根据等效协方差向量r_d和等效椭圆协方差向量r_s中每一个元素所在的行和列在嵌套阵列中对应的阵元位置，计算等效协方差向量r_d中所有元素的维数E_i,j和等效椭圆协方差向量r_s中所有元素的维数F_i,j：

E_i,j＝D(j)-D(i)，

F_i,j＝D(j)+D(i)，

其中，D(i)表示嵌套阵列中第i个阵元的位置，D(j)表示嵌套阵列中第j个阵元的位置。

步骤6：获得虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量

(6a)根据等效协方差向量r_d中所有元素的维数，删除等效协方差向量r_d中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量

(6b)根据等效椭圆协方差向量r_s中所有元素的维数，删除等效椭圆协方差向量r_s中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量

步骤7：根据虚拟均匀阵列协方差向量和椭圆协方差向量构造一个稀疏重构向量r：

步骤8：构造超完备基A_ds(θ)并定义第一空域稀疏向量u与第二空域稀疏向量δ，波达方向角范围θ以及常数矩阵T：

(8a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,将其定义为波达方向角的范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ_q表示第q个刻度值，1≤q≤Q，Q＞＞M；

(8b)构造一个信号稀疏化后对应的C_d×Q维的差阵列超完备基A_d(θ)和一个C_s×Q维的和阵列超完备基A_s(θ)：

A_d(θ)＝[a_d(θ₁),…,a_d(θ_q),…,a_d(θ_Q)]，

A_s(θ)＝[a_s(θ₁),…,a_s(θ_q),…,a_s(θ_Q)]，

其中，α_d(θ_q)表示第q个刻度值对应的差阵列阵列流型向量，α_s(θ_q)表示第q个刻度值对应的和阵列阵列流型向量，其表示公式如下：

α_d(θ_q)＝(exp(j2πd·(C_d-1)/2·sinθ_q/λ),.…,1,…exp[-j2πd·(C_d-1)/2·sinθ_q/λ])^T，

α_s(θ_q)＝(exp(-j2πd·2·sinθ_q/λ),……,exp[-j2πd·C_s·sinθ_q/λ])^T，1≤q≤Q；

(8c)根据差阵列超完备基A_d(θ)与和阵列超完备基A_s(θ)得到超完备基A_ds(θ)：

(8d)定义一个2Q×1维的空域稀疏向量：u＝[u₁,u₂,…,u_q,…,u_Q,u₁,u₂,…,u_q,…,u_Q]^T，定义一个Q×1维的空域稀疏向量：δ＝[δ₁,δ₂,…,δ_q,…,δ_Q]^T，其中u_q和δ_q分别表示第q个第一未知变量和第q个第二未知变量，1≤q≤Q；

(8e)定义一个单位矩阵T_Q为Q×Q的对角阵，即对角元素全为1，其它元素全为0，并定义一个常数矩阵T＝[T_Q,T_Q]。

步骤9：通过稀疏重构方法获得第一空域稀疏向量δ的最优估计

(9a)根据稀疏表示的思想，将波达方向角估计问题转化为求解约束优化方程式：

min||δ||₁

s.t.δ(i)≥||[u(i),u(N+i)]||₂

其中，||·||₁和||·||₂分别表示求向量的1-范数和向量的2-范数，s.t.表示约束关系，β表示误差的允许值，β的取值为现实环境中的噪声标准差,根据求得u，u(i)表示u中第i个元素，θ_i∈θ；

(9b)利用凸优化方法求解上述优化方程，得到第一空域稀疏向量的最优解

步骤10：以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以第一空域稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

以上描述仅是本发明的一个具体实例，不构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不违背本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

Claims (4)

1.基于压缩感知的非圆信号波达方向角估计方法，其特征在于，包括：

(1)用M+N个天线接收机形成嵌套阵列，其中M、N分别表示两个天线接收阵列的阵元数，其取值范围为M≥1，N≥1；

(2)假设空间中有K个非圆目标信号，通过采样和滤波得到嵌套阵列输出信号：Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]^T，其中，y_i(t)表示嵌套阵列的第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)^T表示矩阵转置运算；

(3)利用嵌套阵列输出信号Y(t)，计算虚拟均匀阵列协方差向量和椭圆协方差向量

(4)根据虚拟均匀阵列协方差向量和椭圆协方差向量构造一个稀疏重构向量r：

(5)构造超完备基A_ds(θ)并定义第一空域稀疏向量u与第二空域稀疏向量δ，波达方向角范围θ以及常数矩阵T：

(5.1)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,将其定义为波达方向角的取值空间θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ_q表示第q个刻度值，1≤q≤Q，Q＞＞M；

(5.2)构造一个信号稀疏化后对应的C_d×Q维的差阵列超完备基A_d(θ)和一个C_s×Q维的和阵列超完备基A_s(θ)：

A_d(θ)＝[a_d(θ₁),…,a_d(θ_q),…,a_d(θ_Q)]，

A_s(θ)＝[a_s(θ₁),…,a_s(θ_q),…,a_s(θ_Q)]，

其中，a_d(θ_q)表示第q个刻度值对应的差阵列流型向量，a_s(θ_q)表示第q个刻度值对应的和阵列流型向量，

a_d(θ_q)＝(exp(j2πd·(C_d-1)/2·sinθ_q/λ),…,1,…exp[-j2πd·(C_d-1)/2·sinθ_q/λ])^T

α_s(θ_q)＝(exp(-j2πd·2·sinθ_q/λ),……,exp[-j2πd·C_s·sinθ_q/λ])^T，1≤q≤Q，

(5.3)根据差阵列超完备基A_d(θ)与和阵列超完备基A_s(θ)得到超完备基A_ds(θ)：

(5.4)定义一个2Q×1维的空域稀疏向量：u＝[u₁,u₂,…,u_q,…,u_Q,u₁,u₂,…,u_q,…,u_Q]^T，定义一个Q×1维的空域稀疏向量：δ＝[δ₁,δ₂,…,δ_q,…,δ_Q]^T，其中u_q和δ_q分别表示第一未知变量和第二未知变量，1≤q≤Q；

(5.5)定义一个单位矩阵T_Q为Q×Q的对角阵，即对角元素全为1，其它元素全为0，并定义一个常数矩阵T＝[T_Q,T_Q]；

(6)通过稀疏重构方法获得第一空域稀疏向量δ的最优估计

(7)以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以第一空域稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

2.根据权利要求1所述的方法，其中步骤(1)中用M+N个天线接收机形成嵌套阵列，按如下步骤进行：

(1a)将每个天线接收机称为一个阵元，用M个天线接收机形成第一均匀线性阵列a，其阵元间距为d，定义第一均匀线性阵列a的第一个阵元为起始阵元，定义起始阵元位置D(1)＝1，第一均匀线性阵列a的其它阵元位置依次为D(2)＝2，D(3)＝3，D(4)＝4，…，D(M)＝M；

(1b)用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+1)d，第二均匀线性阵列b的阵元位置依次设置为D(M+1)＝M+1，D(M+2)＝2(M+1)，D(M+2)＝3(M+1)，…，D(M+N)＝N(M+1)，其中，M≥1，N≥1，0＜d≤λ/2，λ为入射到阵列的窄带信号波长；

(1c)将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为Md的位置，将第二均匀线性阵列b的所有阵元依次放置于第一均匀线性阵列a之后，形成嵌套阵列。

3.根据权利要求1所述的方法，其中步骤(3)中计算虚拟均匀阵列协方差向量和椭圆协方差向量按如下步骤进行：

(3a)根据输出信号Y(t)计算协方差矩阵R_d和椭圆协方差矩阵R_s：

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

(3b)将等效协方差矩阵R_d和椭圆协方差矩阵R_s中的元素分别进行排列，得到等效协方差向量r_d和等效椭圆协方差向量r_s：

r_d＝[R_d(1,1),R_d(2,1),…,R_d(M+N,1),R_d(1,2),…,R_d(M+N,2),…,R_d(1,M+N),…,R_d(M+N,M+N)]^T

r_s＝[R_s(1,1),R_s(2,1),…,R_s(M+N,1),R_s(1,2),…,R_s(M+N,2),…,R_s(1,M+N),…,R_s(M+N,M+N)]^T

其中，R_d(i,j)表示协方差矩阵R_d中位于第i行，第j列的元素，i的取值范围为1≤i≤M+N，j的取值范围为1≤j≤M+N；R_s(i,j)表示椭圆协方差矩阵R_s中位于第i行，第j列的元素；

(3c)计算等效协方差向量r_d中所有元素的维数E_i,j和等效椭圆协方差向量r_s中所有元素的维数F_i,j：

E_i,j＝D(j)-D(i)

F_i,j＝D(j)+D(i)

其中，D(i)表示嵌套阵列中第i个阵元的位置，D(j)表示嵌套阵列中第j个阵元的位置；

(3d)删除等效协方差向量r_d中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量删除等效椭圆协方差向量r_s中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量其中中元素个数为C_d＝2MN+2N-1，r_s中元素个数为C_s＝MN+M+N。

4.根据权利要求1所述的方法，其中步骤(6)中计算第一空域稀疏向量的最优解按如下步骤进行：

(6a)根据稀疏表示的思想，将波达方向角估计问题转化为求解约束优化方程式：

min||δ||₁

s.t.δ(i)≥||[u(i),u(N+i)]||₂

其中，||·||₁和||·||₂分别表示求向量的1-范数和向量的2-范数，s.t.表示约束关系，β表示误差的允许值，β的取值为现实环境中的噪声标准差,根据求得u，u(i)表示u中第i个元素，θ_i∈θ；

(6b)利用凸优化方法求解上述优化方程，得到第一空域稀疏向量的最优解