CN106814343B - 一种分步估计的空域信号空间谱估计方法 - Google Patents

Abstract

本发明公开了一种分步估计的空域信号空间谱估计方法，其特征是，根据空域信号具有稀疏性的特点，将空域信号与压缩感知方法结合，对感兴趣的空间划分为较大区域，对空间谱进行一次预估计，通过自适应的方法将以预估计结果为中心的邻域视为感兴趣空间，对此空间进行细分，再次利用压缩感知方法得出空间谱的精确估计。本发明所达到的有益效果：本方法具有在少测量数、低信噪比条件下即可对含有空间上相距较近信号的多目标空间谱进行精确估计，且更加易于硬件存储实现的优点。

Description

一种分步估计的空域信号空间谱估计方法

技术领域

本发明涉及一种分步估计的空域信号空间谱估计方法，属于阵列信号处理技术领域。

背景技术

空域信号的空间谱估计在阵列信号处理中有着重要地位，它涉及雷达、声呐、通信、射电天文等以及医学诊断等多种国民经济和军事等应用领域。因此，对空域信号空间谱进行参数估计受到越来越多的重视。

以多重信号分类(multiple signal classification，MUSIC)算法为代表的现代谱估计方法是空间谱估计向超分辨率发展的一次飞跃，但其需要较大的快拍数，大的快拍导致数据增多，这增加了系统复杂性和计算复杂度。MUSIC算法要求较高的信噪比，并且此方法虽然是一种超分辨率算法，但其对空间上相聚较近的目标分辨效果仍然较差。之后，在此基础上又出现了许多改进算法，如Root-MUSIC,平滑子空间MUSIC算法，使得算法性能显著提高，但它们都有一个共同特点计算量大。

D.Malioutov和M.Cetin等人将稀疏表示的思想引入到DOA估计中。随后，学者们将压缩感知的方法推广到空间谱估计中，提出了基于正交匹配追踪的压缩感知算法，它突破了奈奎斯特定理中信号带宽对于采样频率的限制，具有所需快拍数少，计算简单且在较小信噪比条件下即可对空间谱进行高分辨率估计的优点。但在对含有空间上相距较近信号的多个目标进行谱估计时算法性能较差。

发明内容

为解决现有技术的不足，本发明的目的在于提供一种分步估计的空域信号空间谱估计方法，在阵元数较少条件下，利用少的测量数、在低信噪比时对含有空间上相距较近信号的多个目标空间谱进行精确估计。

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

一种分步估计的空域信号空间谱估计方法，其特征是，包括如下步骤：

1)对感兴趣的空间Θ以步长π/λ₁划分为LN₁份，其满足其中span(Θ)表示空间Θ中最大值和最小值之间的差值，INT表示将一个数值向下取整为最接近的整数的函数，λ₁为预估计步长因子；

这样形成的正交完备稀疏字典构造预估计的稀疏基矩阵并对信号x∈C^N×1进行预估计稀疏表示，Ν为阵源数，表示复数域的N×LN₁维矩阵，其中构造稀疏基目标信号可被稀疏表示为x＝Ψ₁y₁+w₁，其中，为空域信号x的预估计稀疏表示，w₁∈C^N×1为高斯白噪声，j为虚数表示，d为阵元间距，λ为波长，表示空间划分角度，i＝1,2,…,LN₁；

2)利用双结构系统构造预估计测量矩阵其中Φ₃为单位对角矩阵，M₁为预估计测量数；

3)将空域信号x投影到测量矩阵Φ_y上得到观测信号s₁＝Φ_yx＝Φ_y(Ψ₁y₁+w₁)＝T₁y₁+e₁,T₁＝Φ_yΨ₁,e₁＝Φ_yw₁，其中，是预估计观测信号，是预估计恢复矩阵，为预估计观测信号的噪声矢量，M₁为预估计测量数；

4)在得到观测信号后利用OMP算法对步骤3)中的等式进行求解，在OMP算法执行结束后即为预估计结果其中，为预估计所得角度，i＝1、2、…、K，K为稀疏度；

5)缩小搜索范围，生成精确估计正交完备稀疏字典，对信号进行精确估计稀疏表示；在步骤4)得到的预估计值附近进行精确搜索得到空间谱的准确信息，以预估计值为中心，以为半径，自适应生成感兴趣空间其中，为以为中心为半径的中心邻域，依据经验值设定；

6)以为感兴趣空间，以步长π/λ₂划分为LN₂份，其满足其中表示空间中最大值和最小值之间的差值，INT表示将一个数值向下取整为最接近的整数的函数，λ₂为精确估计步长因子，λ₂＞λ₁；构造稀疏基矩阵j为虚数表示，d为阵元间距，λ为波长，表示空间划分角度，i＝1,2,…,LN₂，表示复数域的N×LN₂维矩阵；

7)在新的空间区域用精确估计完备稀疏基将目标信号稀疏表示为x＝Ψ₂y₂+w₂，其中为空域信号x的预估计稀疏表示，w₂∈C^N×1为高斯白噪声；

8)利用双结构系统构造精确测量矩阵其中 Φ₃为单位对角矩阵，M₂为预估计测量数；

9)将空域信号x投影到Φ_j上得到观测信号s₂＝Φ_jx＝Φ_j(Ψ₂y₂+w₂)＝T₂y₂+e₂,T₂＝Φ_jΨ₂,e₂＝Φ_jy₂，其中，是精确估计观测信号，是精确估计恢复矩阵，为精确估计观测信号的噪声矢量；

10)在得到观测信号后利用OMP算法对上述等式进行求解，在OMP算法执行结束后即为精确估计结果AccurateEstimation＝[φ₁ φ₂ … φ_K]，其中，φ_i为精确估计所得角度，i＝1、2、…、K，K为稀疏度。

优选地，所述步骤2)中的双结构系统具体构造步骤如下：

选用Logistic映射，通过映射方程x_n+1＝μx_n(1-x_n)，n＝0,1,2,3...，构造混沌序列{x₀ x₁ … x_n}，上述映射方程中x_n∈(0,1)表示第n次迭代数，n表示混沌序列迭代次数；

舍弃前t(t＜n)个数，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3…，得到序列{z₀ z₁ … z_k}；取其中前(M₁×M₁-1)个值生成一个混沌矩阵M₁是测量数；

将矩阵Γ稀疏化得到其中

构造一个单位对角矩阵将Φ₂、Φ₃组合成为一个新矩阵输出Φ作为预估计测量矩阵Φ_y。

优选地，选取混沌系统参数μ＝4，初值x₀＝0.256。

优选地，所述步骤4)中利用OMP算法的求解步骤为：

41)数据初始化：残差r₀＝s₁，迭代次数inter＝1，T₀为空矩阵；

42)在T₁中选出与残差相关性最大的列：n_inter＝argmax<r_inter-1,t_i>，i＝1,2,…,LN₁，t_i表示T₁的第i列；

43)更新所选列空间：

44)通过对最小二乘问题的求解，保证残差最小，获得在已选列上的最优投影，求解满足的获得估计值；

45)更新残差：

46)更新迭代次数：inter＝inter+1，如果达到最终迭代次数则输出估计值否则返回执行42)。

优选地，所述步骤8)中的双结构系统具体构造步骤如下：

选用Logistic映射，通过映射方程x_n+1＝μx_n(1-x_n)，n＝0,1,2,3...，构造混沌序列{x₀ x₁ … x_n}，上述映射方程中x_n∈(0,1)表示第n次迭代数，n表示混沌序列迭代次数；

舍弃前t(t＜n)个数，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3…，得到序列{z₀ z₁ … z_k}；取其中前(M₂×M₂-1)个值生成一个混沌矩阵M₂是测量数；

将矩阵Γ稀疏化得到其中

构造一个单位对角矩阵将Φ₂、Φ₃组合成为一个新矩阵输出Φ作为预估计测量矩阵Φ_j。

优选地，所述步骤10)中利用OMP算法的求解步骤如下：

101)数据初始化：残差r₀＝s₂，迭代次数inter＝1，T₀为空矩阵；

102)在T₂中选出与残差相关性最大的列：n_inter＝argmax<r_inter-1,t_i>，i＝1,2,…,LN₂，t_i表示T₂的第i列；

103)更新所选列空间：

104)通过对最小二乘问题的求解，保证残差最小，获得在已选列上的最优投影，求解满足的获得估计值；

105)更新残差：

106)更新迭代次数：inter＝inter+1，如果达到最终迭代次数则输出估计值否则返回执行102)。

本发明所达到的有益效果：本方法具有在少测量数、低信噪比条件下即可对含有空间上相距较近信号的多目标空间谱进行精确估计，且更加易于硬件存储实现的优点。

附图说明

图1是Logistic混沌二叉图；

图2是算法框架流程图；

图3是测量矩阵构造流程图；

图4(a)(b)分别是步骤4)和步骤10)中OMP算法流程图；

图5(a)(b)分别是在预估计测量数不同的情况下，不同算法的RMSE与精确估计测量数关系示意图；

图6(a)(b)分别是不同信噪比条件下的精确估计信噪比与RMSE关系示意图；

图7(a)(b)(c)是不同稀疏度条件下精确测量数与RMSE关系示意图；

图8是稀疏度与RMSE关系示意图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

压缩感知中测量矩阵的研究，一直是压缩感知的热点问题，学者们提出了多种测量矩阵，其中高斯矩阵以其较强的随机性、各列之间低的相关性特点普遍应用于压缩感知的各类问题，但其产生的矩阵都是随机的，这对于硬件实现造成了极大的困难。一些学者提出使用混沌序列构造测量矩阵，利用Logistic映射生成测量矩阵已被证明满足约束等距性(RIP)，性能较之高斯矩阵更为优异，然而它属于随机性测量矩阵，具有与高斯矩阵具有相似的性质。

基于此本发明中提出一种新的测量矩阵的构造方法，并在此构造方法的基础上形成一种分步估计的空域信号空间谱估计方法，可简单描述为：第一步对感兴趣的空间划分为较大区域，利用压缩感知方法对空间谱进行一次预估计，第二步通过自适应的方法将以预估计结果为中心的邻域视为感兴趣空间，对此空间进行细分，再次利用压缩感知方法得出空间谱的精确估计。之所以选择这种分步估计的方法在于，在压缩感知稀疏基字典中如果维数过大，将会导致稀疏基矩阵中原子相关性增大，这不利于稀疏信号的恢复，因此在第一步中，对感兴趣的空间范围进行粗略的划分以确保稀疏基维度较小。然而，由于这种对感兴趣空间的粗略划分也导致了空域信号空间谱估计具有较大偏差。因而结合第二步对预估计结果进行局部寻优以增强空间谱估计的准确性，从而得到空域信号空间谱的精确估计。

下面首先对本方法的步骤进行介绍：

步骤1)对感兴趣的空间Θ以步长π/λ₁划分为LN₁份，其满足其中span(Θ)表示空间Θ中最大值和最小值之间的差值，INT表示将一个数值向下取整为最接近的整数的函数，λ₁为预估计步长因子；

这样形成的正交完备稀疏字典构造预估计的稀疏基矩阵并对信号x∈C^N×1进行预估计稀疏表示，Ν为阵源数，表示复数域的N×LN₁维矩阵，其中构造稀疏基目标信号可被稀疏表示为x＝Ψ₁y₁+w₁，其中，为空域信号x的预估计稀疏表示，w₁∈C^N×1为高斯白噪声，j为虚数表示，d为阵元间距，λ为波长，表示空间划分角度，i＝1,2,…,LN₁。

步骤2)利用双结构系统构造预估计测量矩阵其中Φ₃为单位对角矩阵，M₁为预估计测量数，具体构造的过程如下：

选用Logistic映射，通过映射方程x_n+1＝μx_n(1-x_n)，n＝0,1,2,3...，构造混沌序列{x₀ x₁ … x_n}，上述映射方程中x_n∈(0,1)表示第n次迭代数，n表示混沌序列迭代次数；

舍弃前t(t＜n)个数，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3…，得到序列{z₀ z₁ … z_k}；取其中前(M₁×M₁-1)个值生成一个混沌矩阵M₁是测量数；

将矩阵Γ稀疏化得到其中

构造一个单位对角矩阵将Φ₂、Φ₃组合成为一个新矩阵输出Φ作为预估计测量矩阵Φ_y。

本实施例中，优选地选取混沌系统参数μ＝4，初值x₀＝0.256。

步骤3)将空域信号x投影到测量矩阵Φ_y上得到观测信号s₁＝Φ_yx＝Φ_y(Ψ₁y₁+w₁)＝T₁y₁+e₁,T₁＝Φ_yΨ₁，其中，是预估计观测信号，是预估计恢复矩阵，为预估计观测信号的噪声矢量，M₁为预估计测量数。

步骤4)在得到观测信号后利用OMP算法对步骤3)中的等式进行求解，在OMP算法执行结束后即为预估计结果其中，为预估计所得角度，i＝1、2、...、K，K为稀疏度。具体步骤为：

41)数据初始化：残差r₀＝s₁，迭代次数inter＝1，T₀为空矩阵；

42)在T₁中选出与残差相关性最大的列：n_inter＝argmax<r_inter-1,t_i>，i＝1,2,…,LN₁，t_i表示T₁的第i列；

43)更新所选列空间：

44)通过对最小二乘问题的求解，保证残差最小，获得在已选列上的最优投影，求解满足的获得估计值；

45)更新残差：

46)更新迭代次数：inter＝inter+1，如果达到最终迭代次数则输出估计值否则返回执行42)。

步骤5)缩小搜索范围，生成精确估计正交完备稀疏字典，对信号进行精确估计稀疏表示；在步骤4)得到的预估计值附近进行精确搜索得到空间谱的准确信息，以预估计值为中心，以为半径，自适应生成感兴趣空间其中，为以为中心为半径的中心邻域，依据经验值设定；

步骤6)以为感兴趣空间，以步长π/λ₂划分为LN2份，其满足其中表示空间中最大值和最小值之间的差值，INT表示将一个数值向下取整为最接近的整数的函数，λ₂为精确估计步长因子，λ₂＞λ₁；构造稀疏基矩阵j为虚数表示，d为阵元间距，λ为波长，表示空间划分角度，i＝1,2,…,LN₂，表示复数域的N×LN₂维矩阵。对于本步骤的内容，进行补充说明，本步骤的内容与步骤1)基本相似，步骤1)时预估计，而本步骤是精确估计，内容上是有先后关联的，但是具体的计算方法相同，在这里变量上选用了相同的标记，但是不代表取值是与步骤1)相同的，具体的取值不同，能够根据文字的区别来体现，而且，同一个变量，步骤1)所涉及的变量的取值范围在本步骤以及本步骤之后均不会出现，所以不会存在取值范围不清楚的问题。

步骤7)在新的空间区域用精确估计完备稀疏基将目标信号稀疏表示为x＝Ψ₂y₂+w₂，其中为空域信号x的预估计稀疏表示，w₂∈C^N×1为高斯白噪声。

步骤8)利用双结构系统构造精确计测量矩阵其中Φ₃为单位对角矩阵，M₂为预估计测量数。本步骤中的具体内容与步骤2)中相同，因此所涉及的一些中间量均用相同的字母表示，但是对于具体地实施例来说其对应的计算结果不一样。

具体构造步骤如下：

选用Logistic映射，通过映射方程x_n+1＝μx_n(1-x_n)，n＝0,1,2,3...，构造混沌序列{x₀ x₁ … x_n}，上述映射方程中x_n∈(0,1)表示第n次迭代数，n表示混沌序列迭代次数；优选地选取混沌系统参数μ＝4，初值x₀＝0.256。

舍弃前t(t＜n)个数，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3…，得到序列{z₀ z₁ … z_k}；取其中前(M₂×M₂-1)个值生成一个混沌矩阵M₂是测量数；

将矩阵Γ稀疏化得到其中

构造一个单位对角矩阵将Φ₂、Φ₃组合成为一个新矩阵输出Φ作为预估计测量矩阵Φ_j。

步骤9)将空域信号x投影到Φ_j上得到观测信号s₂＝Φ_jx＝Φ_j(Ψ₂y₂+w₂)＝T₂y₂+e₂,T₂＝Φ_jΨ₂，其中，是精确估计观测信号，是精确估计恢复矩阵，为精确估计观测信号的噪声矢量；

步骤10)在得到观测信号后利用OMP算法对上述等式进行求解，在OMP算法执行结束后即为精确估计结果AccurateEstimation＝[φ₁ φ₂ … φ_K]，其中，φ_i为精确估计所得角度，i＝1、2、…、K，K为稀疏度。利用OMP算法的求解步骤如下：

101)数据初始化：残差r₀＝s₂，迭代次数inter＝1，T₀为空矩阵；

102)在T₂中选出与残差相关性最大的列：n_inter＝argmax<r_inter-1,t_i>，i＝1,2,…,LN₂，t_i表示T₂的第i列；

103)更新所选列空间：

104)通过对最小二乘问题的求解，保证残差最小，获得在已选列上的最优投影，求解满足的获得估计值；

105)更新残差：

106)更新迭代次数：inter＝inter+1，如果达到最终迭代次数则输出估计值否则返回执行102)。

对于上面的方法步骤，进行如下的补充说明：

1：本方法在步骤3)和8)中选用研究较为成熟的Logistic映射生成混沌序列，从图1中可以看出，当系统参数μ＝4时，x_n取值会遍历0到1的整个区域，系统进入混沌状态，各点分布具有伪随机性，因此本方法中混沌系统参数优选μ＝4，初值x₀＝0.256。

2：本方法以Logistic混沌序列构造混沌矩阵，在步骤3)和8)中为增强序列随机性，舍弃前t(t＜n)个数，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3...，得到序列{z₀ z₁ … z_k}，该序列是一个以0.5为均值，0.5为对称的伪随机数列，取其中前M×M-1个值生成一个混沌矩阵M为M₁或M₂，利用混沌序列构造的矩阵已被证明满足RIP性质，所以以序列{z₀ z₁ … z_M×M-1}构造的矩阵Γ∈C^M×M也是满足RIP性质的。

在一些文献中混沌矩阵Γ∈C^M×M被直接用作测量矩阵，并且取得了较好的效果，然而这是一个稠密矩阵，其伪随机性使得它与高斯矩阵具有相似的性质，不仅难以硬件实现而且处理大规模数据时速度较慢。考虑到稀疏矩阵在矩阵运算中的优良性能，本方法对混沌矩阵Γ∈C^M×M稀疏处理以优化混沌矩阵性能。

3：在本发明所涉及的构建测量矩阵时引入了稀疏矩阵的概念，稀疏化以后的矩阵占用内存小且混沌系统具有伪随机的特点，这就便于硬件存储和实现。下表反应不同测量数条件下高斯矩阵与本文所提矩阵占用内存比较

测量率	所提矩阵(bytes)	高斯矩阵(bytes)
			0.25	1128	3200
0.5	3000	6400
			0.75	4520	9600
1	6600	12800

表1所提矩阵与高斯矩阵内存比较

由表中数据可看出本文所用矩阵占用内存远小于同等规模高斯矩阵，由于本文测量矩阵是稀疏的，因此所提矩阵的采样率在同等情况下远小于高斯矩阵，这就为处理高维信号问题提供了便利，且节省内存便于硬件存储和实现。

本发明以已知稀疏度的空域信号为模型验证所提方法的可行性与准确性，分别在不同测量数、不同稀疏度、不同信噪比条件下对所提方法进行测试，为直观比较该方法性能优劣，引入高斯随机矩阵，利用高斯随机矩阵单次估计以及分步估计的空间谱估计结果与本文所提方法进行对照比较。

为准确评价算法性能采用蒙特卡罗方法对算法进行仿真，利用均方根误差来描述空间谱估计精度，空间谱估计的均方根误差(RMSE)定义为其中K为稀疏度，CNT为蒙特卡罗循环次数，φ_k,cnt为第k个角度在第cnt次蒙特卡罗实验中所得估计值，θ_k为第k个角度的实际位置。

实施例一：天线阵元数N＝40，稀疏度K＝4，信噪比SNR＝15，信号实际角度信息θ＝[1° 2° -20° 35°]，预估计步长因子λ₁＝30，精确估计步长因子λ₂＝100，预估计测量数分别为M₁＝21、M₁＝25，蒙特卡罗次数为100，在不同精确估计测量数条件下，观察不同算法RMSE随着测量数的变化情况。

在图5中，(a)的预估计数M₁＝21，(b)的预估计数M₁＝25，从图中可看出，空域信号空间谱估计RMSE随着测量数增大而减小，对比(a)、(b)两图，在同样的估计方法下所提矩阵估计精度要好于高斯矩阵，验证了本文所提测量矩阵的可靠性。高斯分步估计曲线与高斯矩阵的单次估计曲线相比，分步估计法所得RMSE在相同测量数条件下小于高斯矩阵单次估计所得结果。从最后的估计结果看本文所提方法的估计误差可以控制在0.5°左右，且此误差主要来源于1°和2°的估计过程。对比(a)、(b)中的高斯矩阵分步估计与本文所提算法，(b)中RMSE均较(a)中有所降低，这是因为随着预估计测量数增大，得到了更加精确的空间谱预估计值，从而增强了精确估计的估计精度。因此可以得出：对于有两个物体较为靠近的苛刻条件下的空域信号，利用一个阵元较少的线阵系统，在观测数较少的条件下本文所提方法即可对空域信号空间谱进行较为精确的估计。

实施例二：天线阵元数N＝40，稀疏度K＝4，预估计和精确估计测量数均为21，信号实际角度信息θ＝[1° 2° -20° 35°]，预估计步长因子λ₁＝30，精确估计步长因子λ₂＝100，蒙特卡罗次数为100，在不同信噪比条件下，观察不同算法RMSE随着SNR的变化情况。

在图6中，(a)中预估计SNR＝10，(b)中预估计SNR＝15，从图中可看出，空域信号空间谱RMSE随着SNR增大而减小，对比高斯分步估计曲线与高斯矩阵的单次估计曲线，分步估计的方法得到的空域信号空间谱估计RMSE更小，因此分步估计的方法在本文条件下性能优于单次估计。比较高斯矩阵分步估计与本文所提算法，可发现本文所提算法在相同条件下空域信号空间谱估计RMSE要小于高斯方法，这说明本文提出的测量矩阵对于小信噪比条件下性能优于高斯矩阵。对比(a)、(b)中的高斯分步估计与本文所提方法，发现(b)中高斯方法RMSE小于同等条件下的(a)中结果，说明空域信号空间谱预估计SNR较大时精确估计结果更加准确，但本文所提方法在(a)、(b)中变化不大，这是因为SNR＝10时预估计的空域信号空间谱已经较为精确，与SNR＝15所得预估计结果相差不大。

实施例三：天线阵元数N＝40，预估计步长因子λ₁＝30，精确估计步长因子λ₂＝100，蒙特卡罗次数为100，信噪比SNR＝15。

图(7)表示不同稀疏度条件下精确测量数与RMSE的关系，(a)中K＝2，信号实际角度信息θ＝[1° 2°]，(b)中K＝4，信号实际角度信息θ＝[1° 2° -20° 35°]，(c)中K＝6，信号实际角度信息θ＝[1° 2° -20° 35° -40° 70°]，在(a)、(b)、(c)中预估计测量数均为21。

图(8)表示不同稀疏度的空域信号在固定测量数下的RMSE，其中各稀疏信号均包括来波方向为1°,2°的信号，其余为中(除1°和2°)方向随机的信号，预估计测量数为21，精确估计测量数为35。

观察图(7)，在(a)中本文所提方法与高斯矩阵分步估计方法性能相差不大，但观察(b)、(c)可发现随着测量数增大本文所提方法的空域信号空间谱估计RMSE均小于高斯矩阵分步估计及高斯矩阵单次估计，从实验数据来看估计误差主要来自于对1°和2°的估计。

在图(8)中可直观看到，空域信号空间谱估计在固定测量数条件下增大信号稀疏度会增大谱估计的误差，在相同条件下，所提算法空域信号空间谱估计性能好于高斯分步估计，图中，高斯矩阵分步估计在稀疏度大于4后RMSE迅速增大，而本文所提算法在稀疏度为6时精确度虽然有所变差但RMSE远小于同等条件下的高斯矩阵分步估计结果。因此，所提方法在目标数较多的条件下同样具有较好的性能。

以上实验结果表明，本方法在少的测量数、小的信噪比条件下即可对空间谱进行精确估计，且在信号条件较为苛刻即含有非常靠近的多个目标条件下也能取得很好的效果，验证了此方法的可靠性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (6)

1.一种分步估计的空域信号空间谱估计方法，其特征是，包括如下步骤：

1)对感兴趣的空间Θ以步长π/λ₁划分为LN₁份，其满足其中span(Θ)表示空间Θ中最大值和最小值之间的差值，INT表示将一个数值向下取整为最接近的整数的函数，λ₁为预估计步长因子；

这样形成的正交完备稀疏字典构造预估计的稀疏基矩阵并对信号x∈C^N×1进行预估计稀疏表示，Ν为阵源数，表示复数域的N×LN₁维矩阵，其中构造稀疏基目标信号可被稀疏表示为x＝Ψ₁y₁+w₁，其中，为空域信号x的预估计稀疏表示，w₁∈C^N×1为高斯白噪声，j为虚数表示，d为阵元间距，λ为波长，表示空间划分角度，i＝1,2,…,LN₁；

2)利用双结构系统构造预估计测量矩阵其中Φ₁₃由M₁×M₁单位对角矩阵中取N-M₁列构成，N-M₁≤M₁，M₁为预估计测量数；

3)将空域信号x投影到测量矩阵Φ_y上得到观测信号s₁＝Φ_yx＝Φ_y(Ψ₁y₁+w₁)＝T₁y₁+e₁,T₁＝Φ_yΨ₁,e₁＝Φ_yw₁，其中，是预估计观测信号，是预估计恢复矩阵，为预估计观测信号的噪声矢量，

4)在得到观测信号后利用OMP算法对步骤3)中的等式进行求解，在OMP算法执行结束后即为预估计结果其中，为预估计所得角度，i＝1、2、…、K，K为稀疏度；

5)缩小搜索范围，生成精确估计正交完备稀疏字典，对信号进行精确估计稀疏表示；在步骤4)得到的预估计值附近进行精确搜索得到空间谱的准确信息，以预估计值为中心，以为半径，自适应生成感兴趣空间其中，是以为中心为半径的中心邻域，依据经验值设定；

6)以为感兴趣空间，以步长π/λ₂划分为LN₂份，其满足其中表示空间中最大值和最小值之间的差值，INT表示将一个数值向下取整为最接近的整数的函数，λ₂为精确估计步长因子，λ₂＞λ₁；构造稀疏基矩阵j为虚数表示，d为阵元间距，λ为波长，表示空间划分角度，i＝1,2,…,LN₂，表示复数域的N×LN₂维矩阵；

7)在新的空间区域用精确估计完备稀疏基将目标信号稀疏表示为x＝Ψ₂y₂+w₂，其中为空域信号x的精确估计稀疏表示，w₂∈C^N×1为高斯白噪声；

8)利用双结构系统构造精确测量矩阵其中 Φ₂₃由M₂×M₂单位对角矩阵中取N-M₂列构成，N-M₂≤M₂，M₂为精确测量数；

9)将空域信号x投影到Φ_j上得到观测信号s₂＝Φ_jx＝Φ_j(Ψ₂y₂+w₂)＝T₂y₂+e₂,T₂＝Φ_jΨ₂，其中，是精确估计观测信号，是精确估计恢复矩阵，为精确估计观测信号的噪声矢量；

10)在得到观测信号后利用OMP算法对步骤9)中的等式进行求解，在OMP算法执行结束后即为精确估计结果AccurateEstimation＝[φ₁ φ₂ … φ_K]，其中，φ_i为精确估计所得角度，i＝1、2、…、K，K为稀疏度。

2.根据权利要求1所述的一种分步估计的空域信号空间谱估计方法，其特征是，所述步骤2)中的双结构系统具体构造步骤如下：

选用Logistic映射，通过映射方程x_n+1＝μx_n(1-x_n)，n＝0,1,2,3...，构造混沌序列{x₀x₁ … x_n}，上述映射方程中x_n∈(0,1)表示第n次迭代数，n表示混沌序列迭代次数，μ表示混沌系统参数；

舍弃前t个数，t＜n，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3…，得到序列{z₀ z₁ … z_k}；取其中前(M₁×M₁-1)个值生成一个混沌矩阵M₁是预估计测量数；

将矩阵Γ稀疏化得到其中

构造一个单位对角矩阵将Φ₁₂、Φ₁₃组合成为一个新矩阵输出Φ作为预估计测量矩阵Φ_y。

3.根据权利要求1所述的一种分步估计的空域信号空间谱估计方法，其特征是，所述步骤4)中利用OMP算法的求解步骤为：

41)数据初始化：残差r₀＝s₁，迭代次数inter＝1，T₀为空矩阵；

42)在T₁中选出与残差相关性最大的列：n_inter＝argmax<r_inter-1,t_i>，i＝1,2,…,LN₁，t_i表示T₁的第i列；

43)更新所选列空间：

44)通过对最小二乘问题的求解，保证残差最小，获得在已选列上的最优投影，求解满足的获得估计值；

45)更新残差：

46)更新迭代次数：inter＝inter+1，如果达到最终迭代次数则输出估计值否则返回执行42)。

4.根据权利要求1所述的一种分步估计的空域信号空间谱估计方法，其特征是，所述步骤8)中的双结构系统具体构造步骤如下：

选用Logistic映射，通过映射方程x_n+1＝μx_n(1-x_n)，n＝0,1,2,3...，构造混沌序列{x₀x₁ … x_n}，上述映射方程中x_n∈(0,1)表示第n次迭代数，n表示混沌序列迭代次数，μ表示混沌系统参数；

舍弃前t个数，t＜n，生成新的序列{x_t x_t+1 … x_n}，对此混沌序列以间隔d进行等间隔抽样得到z_k＝x_t+kd，k＝0,1,2,3…，得到序列{z₀ z₁ … z_k}；取其中前(M₂×M₂-1)个值生成一个混沌矩阵M₂是精确估计测量数；

将矩阵Γ稀疏化得到其中

构造一个单位对角矩阵将Φ₂₂、Φ₂₃组合成为一个新矩阵输出Φ作为精确估计测量矩阵Φ_j。

5.根据权利要求1所述的一种分步估计的空域信号空间谱估计方法，其特征是，所述步骤10)中利用OMP算法的求解步骤如下：

101)数据初始化：残差r₀＝s₂，迭代次数inter＝1，T₀为空矩阵；

102)在T₂中选出与残差相关性最大的列：n_inter＝argmax<r_inter-1,t_i>，i＝1,2,…,LN₂，t_i表示T₂的第i列；

103)更新所选列空间：

104)通过对最小二乘问题的求解，保证残差最小，获得在已选列上的最优投影，求解满足的获得估计值；

105)更新残差：

106)更新迭代次数：inter＝inter+1，如果达到最终迭代次数则输出估计值否则返回执行102)。

6.根据权利要求2或4所述的一种分步估计的空域信号空间谱估计方法，其特征是，选取混沌系统参数μ＝4，初值x₀＝0.256。