本发明公开了一种基于压缩感知的非圆信号波达方向角估计方法，主要解决现有技术中阵元利用率低，信号识别数量少的问题，其实现步骤是：1)获取嵌套阵列输出信号，计算协方差矩阵和椭圆协方差矩阵，构造相应的协方差向量，计算这两个向量中元素的维数，构造相应的虚拟阵列向量，构造稀疏重构向量；2)根据稀疏重构向量构造超完备基并定义第一、第二空域稀疏向量，波达方向角范围以及常数矩阵；3)通过稀疏重构方法获得第一空域稀疏向量的最优估计，绘制幅度谱图，并寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求角度值。本发明在阵元数量有限的情况下大大提高了阵列可识别的信源数目，适用于目标侦察和无源定位。
1.基于压缩感知的非圆信号波达方向角估计方法，其特征在于，包括：
(1)用M+N个天线接收机形成嵌套阵列，其中M、N分别表示两个天线接收阵列的阵元数，其取值范围为M≥1，N≥1；
(2)假设空间中有K个非圆目标信号，通过采样和滤波得到嵌套阵列输出信号：Y(t)＝[y1(t),…,yi(t),…,yM+N(t)]T，其中，yi(t)表示嵌套阵列的第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)T表示矩阵转置运算；
(3)利用嵌套阵列输出信号Y(t)，计算虚拟均匀阵列协方差向量 和椭圆协方差向量(4)根据虚拟均匀阵列协方差向量 和椭圆协方差向量 构造一个稀疏重构向量r：(5)构造超完备基Ads(θ)并定义第一空域稀疏向量u与第二空域稀疏向量δ，波达方向角范围θ以及常数矩阵T：(5.1)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,将其定义为波达方向角的取值空间θ＝[θ1,θ2,…,θq,…,θQ]，θq表示第q个刻度值，1≤q≤Q，Q＞＞M；
(5.2)构造一个信号稀疏化后对应的Cd×Q维的差阵列超完备基Ad(θ)和一个Cs×Q维的和阵列超完备基As(θ)：Ad(θ)＝[ad(θ1),…,ad(θq),…,ad(θQ)]，
As(θ)＝[as(θ1),…,as(θq),…,as(θQ)]，
其中，ad(θq)表示第q个刻度值对应的差阵列流型向量，as(θq)表示第q个刻度值对应的和阵列流型向量，ad(θq)＝(exp(j2πd·(Cd-1)/2·sinθq/λ),…,1,…exp[-j2πd·(Cd-1)/2·sinθq/λ])Tαs(θq)＝(exp(-j2πd·2·sinθq/λ),……,exp[-j2πd·Cs·sinθq/λ])T，1≤q≤Q，(5.3)根据差阵列超完备基Ad(θ)与和阵列超完备基As(θ)得到超完备基Ads(θ)：(5.4)定义一个2Q×1维的空域稀疏向量：u＝[u1,u2,…,uq,…,uQ,u1,u2,…,uq,…,uQ]T，定义一个Q×1维的空域稀疏向量：δ＝[δ1,δ2,…,δq,…,δQ]T，其中uq和δq分别表示第一未知变量和第二未知变量，1≤q≤Q；
(5.5)定义一个单位矩阵TQ为Q×Q的对角阵，即对角元素全为1，其它元素全为0，并定义一个常数矩阵T＝[TQ,TQ]；
(6)通过稀疏重构方法获得第一空域稀疏向量δ的最优估计
(7)以波达方向角范围θ＝[θ1,θ2,…,θ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)计算协方差矩阵Rd和椭圆协方差矩阵Rs：其中，(·)H表示共轭转置运算；
(3b)将等效协方差矩阵Rd和椭圆协方差矩阵Rs中的元素分别进行排列，得到等效协方差向量rd和等效椭圆协方差向量rs：rd＝[Rd(1,1),Rd(2,1),…,Rd(M+N,1),Rd(1,2),…,Rd(M+N,2),…,Rd(1,M+N),…,Rd(M+N,M+N)]Trs＝[Rs(1,1),Rs(2,1),…,Rs(M+N,1),Rs(1,2),…,Rs(M+N,2),…,Rs(1,M+N),…,Rs(M+N,M+N)]T其中，Rd(i,j)表示协方差矩阵Rd中位于第i行，第j列的元素，i的取值范围为1≤i≤M+N，j的取值范围为1≤j≤M+N；Rs(i,j)表示椭圆协方差矩阵Rs中位于第i行，第j列的元素；
(3c)计算等效协方差向量rd中所有元素的维数Ei,j和等效椭圆协方差向量rs中所有元素的维数Fi,j：Ei,j＝D(j)-D(i)
其中，D(i)表示嵌套阵列中第i个阵元的位置，D(j)表示嵌套阵列中第j个阵元的位置；
(3d)删除等效协方差向量rd中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量 删除等效椭圆协方差向量rs中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量 其中 中元素个数为Cd＝2MN+2N-1，rs中元素个数为Cs＝MN+M+N。
4.根据权利要求1所述的方法，其中步骤(6)中计算第一空域稀疏向量的最优解 按如下步骤进行：(6a)根据稀疏表示的思想，将波达方向角估计问题转化为求解约束优化方程式：min||δ||1
s.t.δ(i)≥||[u(i),u(N+i)]||2
其中，||·||1和||·||2分别表示求向量的1-范数和向量的2-范数，s.t.表示约束关系，β表示误差的允许值，β的取值为现实环境中的噪声标准差,根据 求得u，u(i)表示u中第i个元素，θi∈θ；
(6b)利用凸优化方法求解上述优化方程，得到第一空域稀疏向量的最优解
基于压缩感知的非圆信号波达方向角估计方法
技术领域
[0001] 本发明属于信号处理技术领域，特别涉及一种电磁信号的阵列信号波达方向角估计方法，可用于对飞机、舰船运动目标的侦察与无源定位。
背景技术
[0002] 信号的波达方向角DOA估计是阵列信号处理领域的一个重要分支，它是指利用天线阵列对空间声学信号、电磁信号进行感应接收，再运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。
[0003] 在现代通信中，二相相移键控以及M进制幅移键控等非圆信号的应用越来越多，因此有关非圆信号的DOA估计受到了越来越多的关注。目前已有关于利用阵列天线处理非圆信号的一些方法被提出，比较有代表性的是Zhangmeng Liu等人发表的论文“Direction-of-arrival estimation of noncircular signals via sparse representation”(《IEEE Transactions on aerospace and electronic systems》,VOL 48,No.3,pp.2690-2698 2012.07)中公开了一种利用压缩感知求解进行非圆信号DOA估计的方法，该方法是基于均匀阵列的。
[0004] 另一方面，为了在较少的阵元条件下得到尽量大的角度自由度，检测更多的信源，一些新的阵列结构被提出，比较有代表性的是嵌套阵列以及互质阵列。P Piya等人在其发表的论文“Nested Arrays:A Novel Approach to Array Processing With Enhanced Degrees of Freedom”(《IEEE transactions on signal processing》,VOL 58,NO.8,August 2010)中公开了一种基于嵌套阵列的DOA估计方法，该方法能够使用M+N个阵元，生成2MN+2N-1个虚拟阵元，可检测MN+N-1个信号。该方法具有估计多于阵元数目的信号数的能力，但是，该阵列的讨论都集中在接收信号为圆信号的条件下，对于如何利用该阵列进行非圆信号的处理目前还没有研究。
[0005] 在实际应用中，对于具有非圆信号，在给定一定数量阵元的环境下，如果不能合理利用这些阵元以及信号的非圆特性，就不能估计足够多的信号，造成侦察和定位资源的浪费。
发明内容
[0006] 本发明的目的在于针对上述现有技术存在的不足，提出一种基于压缩感知的非圆信号波达方向角估计方法，以在非圆信号环境下，利用嵌套阵列进行信号处理算法设计，避免因不能合理利用阵元和信号特性造成的资源浪费。
[0007] 为实现上述目的，本发明技术方案包括如下：
[0008] (1)用M+N个天线接收机形成嵌套阵列，其中M、N分别表示两个天线接收阵列的阵元数，其取值范围为M≥1，N≥1；
[0009] (2)假设空间中有K个非圆目标信号，通过采样和滤波得到嵌套阵列输出信号：Y(t)＝[y1(t),…,yi(t),…,yM+N(t)]T，其中，yi(t)表示嵌套阵列的第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)T表示矩阵转置运算；
[0010] (3)利用嵌套阵列输出信号Y(t)，计算虚拟均匀阵列协方差向量 和椭圆协方差向量
[0011] (4)根据虚拟均匀阵列协方差向量 和椭圆协方差向量 构造一个稀疏重构向量r：
[0012]
[0013] (5)构造超完备基Ads(θ)并定义第一空域稀疏向量u与第二空域稀疏向量δ，波达方向角范围θ以及常数矩阵T：
[0014] (5.1)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,将其定义为波达方向角的范围θ＝[θ1,θ2,…,θq,…,θQ]，θq表示第q个刻度值，1≤q≤Q，Q＞＞M；
[0015] (5.2)构造一个信号稀疏化后对应的Cd×Q维的差阵列超完备基Ad(θ)和一个Cs×Q维的和阵列超完备基As(θ)：
[0016] Ad(θ)＝[ad(θ1),…,ad(θq),…,ad(θQ)]，
[0017] As(θ)＝[as(θ1),…,as(θq),…,as(θQ)]，
[0018] 其中，αd(θq)表示第q个刻度值对应的差阵列阵列流型向量，αs(θq)表示第q个刻度值对应的和阵列阵列流型向量，
[0019] αd(θq)＝(exp(j2πd·(Cd-1)/2·sinθq/λ),.…,1,…exp[-j2πd·(Cd-1)/2·sinθq/λ])T，
[0020] αs(θq)＝(exp(-j2πd·2·sinθq/λ),……,exp[-j2πd·Cs·sinθq/λ])T，1≤q≤Q，[0021] (5.3)根据差阵列超完备基Ad(θ)与和阵列超完备基As(θ)得到超完备基Ads(θ)：
[0022]
[0023] (5.4)定义一个2Q×1维的空域稀疏向量：u＝[u1,u2,…,uq,…,uQ,u1,u2,…,uq,…,uQ]T，定义一个Q×1维的空域稀疏向量：δ＝[δ1,δ2,…,δq,…,δQ]T，其中uq和δq分别表示第一未知变量和第二未知变量，1≤q≤Q；
[0024] (5.5)定义一个单位矩阵TQ为Q×Q的对角阵，即对角元素全为1，其它元素全为0，并定义一个常数矩阵T＝[TQ,TQ]；
[0025] (6)通过稀疏重构方法获得第一空域稀疏向量δ的最优估计
[0026] (7)以波达方向角范围θ＝[θ1,θ2,…,θq,…,θQ]的值为x轴坐标，以第一空域稀疏向量的最优解 的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。
[0027] 本发明与现有技术相比具有以下优点：
[0028] 1)本发明采用了嵌套阵列模型进行波达方向角度估计，克服了现有技术中采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目的缺点，提高了在阵元数目相同的条件下的阵列可识别信源数目。
[0029] 2)本发明利用信号的协方差矩阵Rd，以及椭圆协方差矩阵Rs对信号进行估计，增多了可估计的非圆信号数目。
[0030] 3)本发明利用了嵌套阵列的特点和非圆信号的特性，将非圆信号在嵌套阵列上进行信号处理，在M+N个阵元上能估计出(MN+M+N-1)/2+MN+N-1个信号，大大提高了阵列的利用率，增加了阵列可识别信源的数目。
[0031] 4)在波达角度估计中引入压缩感知算法，提高了角度估计的精度以及在较低的信噪比条件下的估计性能。
附图说明
[0032] 图1是本发明的实现流程图；
[0033] 图2是本发明中嵌套阵列的结构示意图。
具体实施方式
[0034] 参照图1，本实例的实现步骤如下：
[0035] 步骤1：用M+N个天线接收机形成嵌套阵列。
[0036] (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，λ为入射到阵列的窄带信号波长；
[0037] (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；
[0038] (1c)将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为Md的位置，将第二均匀线性阵列b的所有阵元依次放置于第一均匀线性阵列a之后，形成嵌套阵列，如图2所示。
[0039] 步骤2：获得嵌套阵列输出信号Y(t)。
[0040] 假设空间中有K个非圆目标信号，使用嵌套阵列天线接收机，对空间目标信号进行快拍采样和匹配滤波操作，得到嵌套阵列输出信号：Y(t)＝[y1(t),…,yi(t),…,yM+N(t)]T，其中，K的取值范围是K＜MN+M+N-1，yi(t)表示嵌套阵列第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)T表示矩阵转置运算。
[0041] 步骤3：计算协方差矩阵Rd和椭圆协方差矩阵Rs。
[0042] 利用嵌套阵列输出信号Y(t)，计算协方差矩阵Rd和椭圆协方差矩阵Rs：
[0043]
[0044]
[0045] 其中，(·)H表示矩阵共轭转置运算。
[0046] 步骤4：构造等效协方差向量rd和等效椭圆协方差向量rs。
[0047] 根据协方差矩阵Rd和椭圆协方差矩阵Rs中的元素，分别构造等效协方差向量rd和等效椭圆协方差向量rs：
[0048]
[0049]
[0050] 其中，Rd(i,j)表示协方差矩阵Rd中位于第i行、第j列的元素，i的取值范围为1≤i≤M+N，j的取值范围为1≤j≤M+N，Rs(i,j)表示椭圆协方差矩阵Rs中位于第i行，第j列的元素。
[0051] 步骤5：计算等效协方差向量和等效椭圆协方差向量中所有元素的维数。
[0052] 根据等效协方差向量rd和等效椭圆协方差向量rs中每一个元素所在的行和列在嵌套阵列中对应的阵元位置，计算等效协方差向量rd中所有元素的维数Ei,j和等效椭圆协方差向量rs中所有元素的维数Fi,j：
[0053] Ei,j＝D(j)-D(i)，
[0054] Fi,j＝D(j)+D(i)，
[0055] 其中，D(i)表示嵌套阵列中第i个阵元的位置，D(j)表示嵌套阵列中第j个阵元的位置。
[0056] 步骤6：获得虚拟均匀阵列协方差向量 和虚拟均匀阵列椭圆协方差向量[0057] (6a)根据等效协方差向量rd中所有元素的维数，删除等效协方差向量rd中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量
[0058] (6b)根据等效椭圆协方差向量rs中所有元素的维数，删除等效椭圆协方差向量rs中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量
[0059] 步骤7：根据虚拟均匀阵列协方差向量 和椭圆协方差向量 构造一个稀疏重构向量r：
[0060]
[0061] 步骤8：构造超完备基Ads(θ)并定义第一空域稀疏向量u与第二空域稀疏向量δ，波达方向角范围θ以及常数矩阵T：
[0062] (8a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,将其定义为波达方向角的范围θ＝[θ1,θ2,…,θq,…,θQ]，θq表示第q个刻度值，1≤q≤Q，Q＞＞M；
[0063] (8b)构造一个信号稀疏化后对应的Cd×Q维的差阵列超完备基Ad(θ)和一个Cs×Q维的和阵列超完备基As(θ)：
[0064] Ad(θ)＝[ad(θ1),…,ad(θq),…,ad(θQ)]，
[0065] As(θ)＝[as(θ1),…,as(θq),…,as(θQ)]，
[0066] 其中，αd(θq)表示第q个刻度值对应的差阵列阵列流型向量，αs(θq)表示第q个刻度值对应的和阵列阵列流型向量，其表示公式如下：
[0067] αd(θq)＝(exp(j2πd·(Cd-1)/2·sinθq/λ),.…,1,…exp[-j2πd·(Cd-1)/2·sinθq/λ])T，
[0068] αs(θq)＝(exp(-j2πd·2·sinθq/λ),……,exp[-j2πd·Cs·sinθq/λ])T，1≤q≤Q；
[0069] (8c)根据差阵列超完备基Ad(θ)与和阵列超完备基As(θ)得到超完备基Ads(θ)：
[0070]
[0071] (8d)定义一个2Q×1维的空域稀疏向量：u＝[u1,u2,…,uq,…,uQ,u1,u2,…,uq,…,uQ]T，定义一个Q×1维的空域稀疏向量：δ＝[δ1,δ2,…,δq,…,δQ]T，其中uq和δq分别表示第q个第一未知变量和第q个第二未知变量，1≤q≤Q；
[0072] (8e)定义一个单位矩阵TQ为Q×Q的对角阵，即对角元素全为1，其它元素全为0，并定义一个常数矩阵T＝[TQ,TQ]。
[0073] 步骤9：通过稀疏重构方法获得第一空域稀疏向量δ的最优估计
[0074] (9a)根据稀疏表示的思想，将波达方向角估计问题转化为求解约束优化方程式：
[0075] min||δ||1
[0076] s.t.δ(i)≥||[u(i),u(N+i)]||2
[0077]
[0078] 其中，||·||1和||·||2分别表示求向量的1-范数和向量的2-范数，s.t.表示约束关系，β表示误差的允许值，β的取值为现实环境中的噪声标准差,根据求得u，u(i)表示u中第i个元素，θi∈θ；
[0079] (9b)利用凸优化方法求解上述优化方程，得到第一空域稀疏向量的最优解[0080] 步骤10：以波达方向角范围θ＝[θ1,θ2,…,θq,…,θQ]的值为x轴坐标，以第一空域稀疏向量的最优解 的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。
[0081] 以上描述仅是本发明的一个具体实例，不构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不违背本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。
