CN110244272A - 基于秩一去噪模型的波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于秩一去噪模型的波达方向估计方法，该方法的具体步骤为：建立雷达的接收信号模型，确定接收信号的测量矩阵；根据接收信号的测量矩阵，构造基于秩一去噪模型的信号协方差矩阵；对基于秩一去噪模型的信号协方差矩阵进行稀疏重构，得到稀疏信号矢量；根据稀疏重构信号矢量，采用交替网格优化算法估计目标信源的波达方向。本发明在未知非均匀噪声条件下，基于秩一去噪模型，通过重新设计降维信号的协方差矩阵消除非均匀噪声，然后利用矢量化操作来获得等效的源矢量，最后采用确定信号的波达方向，实现有限快拍数的高分辨率波达方向估计和抗非均匀噪声的高鲁棒性。

Description

基于秩一去噪模型的波达方向估计方法

技术领域

本发明属于雷达信号处理技术领域，尤其涉及基于秩一去噪模型的波达方向估计方法。

背景技术

通过传感器阵列确定信号的波达方向(DOA)是阵列信号处理领域中的基本问题，涉及领域包括雷达，声纳，通信和射电天文学。

现有研究表明，用于波达方向估计的传统信号子空间方法，例如多信号分类(MUSIC)算法和旋转不变子空间(ESPRIT)算法，在具有足够快拍的高斯白噪声的条件下，可以实现相对高精度和超高分辨率的波达方向估计。而匹配追踪(MP)算法和正交匹配追踪(OMP)算法，则利用传感器阵列获得信号的稀疏特性，通过稀疏重构算法极大地改善了波达方向估计的性能。

基于压缩感知理论的L1奇异值分解(L1-SVD)算法，该算法利用1范数约束和SVD可以在低信噪比环境下提高DOA估计性能。然而，上述稀疏重构算法的计算复杂度大，对于快拍数的要求也相当高。一种基于具有反馈的零空间调整算法(NST+HT+FB)的交替网格优化(AGO)算法解决了该问题。但是，这种方法只是简单地考虑接收信号的第一列互相关信息，同时也会导致阵列孔径减小，从而造成在DOA估计分辨率方面有某种程度上的性能损失。

上述基于子空间和稀疏重建的算法仅在满足均匀白噪声的假设时才能很好地执行，在均匀白噪声中，噪声协方差矩阵是缩放的单位矩阵。但是在某些应用中，例如稀疏阵列，由传感器阵列获得的噪声协方差矩阵彼此不相同，也就是噪声在空间上是白色的，而方差彼此不相同，就会产生非均匀高斯白噪声。此外，随机非均匀ML-DOA估计方法，类似于确定性ML估计器中使用的迭代方式。但是基于ML的非均匀DOA估计器由于包含迭代过程和高度非线性优化问题而耗费大量时间，而且该方法的性能与初始值密切相关，计算复杂度大，稳定性差。

一种用于确定噪声协方差矩阵的DOA估计方法可以改善非均匀噪声环境中的DOA估计性能。该方法将接收信号协方差矩阵划分为多个子矩阵，然后根据子矩阵之间的关系得到噪声协方差矩阵；以迭代方式确定信号/噪声子空间的DOA估计方法的原理是通过求解对数似然函数或最小二乘问题来估计噪声协方差矩阵和信号/噪声子空间，然后可以从估计的信号/噪声子空间或噪声协方差矩阵完成DOA估计。但是该方法由于包含迭代过程导致耗时较大，且移除操作可能导致信号协方差矩阵中包含的部分信息丢失。

上述传统的稀疏重构算法对于密集信源和多目标的DOA估计性能不是很理想，且计算复杂度大。

发明内容

针对上述问题，本发明的目的是提出一种基于秩一去噪模型的波达方向估计方法。本发明在未知非均匀噪声条件下可估计波达方向，该算法基于秩一去噪模型，通过重新设计降维信号的协方差矩阵消除非均匀噪声，然后利用矢量化操作来获得等效的源矢量，最后采用交替网格优化算法(AGO)确定信号的波达方向，实现有限快拍数的高分辨率波达方向估计和抗非均匀噪声的高鲁棒性。

为了达到上述目的，本发明采用以下技术方案予以解决。

基于秩一去噪模型的波达方向估计方法，包括以下步骤：

步骤1，建立雷达的接收信号模型x(t)，确定接收信号的测量矩阵X。

(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列，则t时刻接收到的目标信源的回波信号x(t)，其表达式为：

x(t)＝As(t)+n(t)；

其中，A为导向矩阵，s(t)是信号波形矢量，n(t)＝[n₁(t),n₂(t),…,n_M(t)]^T为零均值加性非均匀复数高斯白噪声向量，且n(t)～CN(0,Q)，Q是与n(t)相关的噪声协方差矩阵。

(1.2)导向矩阵A的表达式为：

A＝[a(θ₁),a(θ₂),a(θ_l),…,a(θ_L)]；

其中，L是目标信源数，l＝1,2,...,L；a(θ_l)是M×1的导向矢量，M代表雷达阵元数，θ_l代表第l个入射信号角度，θ_l∈Θ，Θ表示角度搜索范围，其表达式为：

a(θ_l)＝[1,e^-jα,…,e^-j(M-1)α]^T；

其中，α表示对应第l个DOA的阵列几何形状，α＝2πdsin(θ_l)/λ，d代表阵元间距，λ代表波长；

s(t)是信号波形矢量，其表达式为：

s(t)＝[s₁(t),s₂(t),…,s_L(t)]^T；

(1.3)根据接收信号模型x(t)，得到接收信号的测量矩阵X：

其中，J是快拍数，t_z表示第z个快拍时刻，x(t_z)表示t_z快拍时刻天线阵列所接收的M×1维回波信号数据，x₁(t_z)表示t_z快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据，*表示共轭操作。

步骤2，根据接收信号的测量矩阵X，构造基于秩一去噪模型的信号协方差矩阵

(2.1)根据接收信号的测量矩阵X，构建基于秩一去噪模型的信号协方差矩阵的第m列其表达式为：

其中，E_m是统计期望；是一个秩一相关矢量；X((1:M)_m,t_z)表示除去第m个阵元的测量矩阵；是第m个阵元的测量矩阵的共轭转置；t_z表示第z个快拍时刻。

(2.2)由M个秩一相关矢量得到基于秩一去噪模型的信号协方差矩阵的组合表达式：

其中，A((1:M)₁,:)表示导向矩阵中第1到M行中的第1行的所有列；表示信号功率矢量；P_l表示第l个信号功率。

(2.3)根据基于秩一去噪模型的信号协方差矩阵的组合表达式，当矩阵中的行序号与列序号不相等时，即x≠y，基于秩一去噪模型的信号协方差矩阵的第x行第y列元素的表达式为：

(2.4)根据基于秩一去噪模型的信号协方差矩阵的第x行第y列元素的表达式，得到基于秩一去噪模型的信号协方差矩阵为：

其中，表示L个信号功率的总和， j表示虚数单位，M代表雷达阵元数，θ_l代表第l个入射信号角度，d代表阵元间距，λ代表波长。

步骤3，对基于秩一去噪模型的信号协方差矩阵进行稀疏重构，得到稀疏信号矢量Υ。

(3.1)根据基于秩一去噪模型的信号协方差矩阵信号构造一个维数为2(M-1)×1的信号矢量Υ'，信号矢量Υ'的第项的表达式为：

其中，表示的第行第列的元素。

(3.2)将信号矢量r'进行线性表示为：

Υ'＝B(θ)P；

其中，P表示信号功率矢量，P＝[P₁,P₂,…,P_L]^T；B(θ)表示虚拟流形矩阵，B(θ)＝[b(θ₁),…,b(θ_l),…,b(θ_L)]，第l列导向矢量b(θ_l)为j表示虚数单位，M代表雷达阵元数，θ_l代表第l个入射信号角度，d代表阵元间距，λ代表波长。

(3.3)根据实际回波信号在空域中的稀疏性，采用重构算法对信号矢量Υ'进行稀疏重构，得到稀疏重构信号矢量Υ：

其中，表示划分的一组字典角度，一般 是具有2(M-1)＜＜N_θ属性的过完备字典，且满足受限等距属性(RIP)，表示字典角度为时的稀疏导向矢量，即第1列稀疏导向矢量，N_θ是字典数；表示稀疏信号功率矢量，其仅有L列非零系数，N_θ＞＞L，表示第1个稀疏信号功率。

步骤4，根据稀疏重构信号矢量Υ，采用交替网格优化算法估计目标信源的波达方向。

与现有技术相比，本发明的有益效果为：

(1)本发明基于秩一去噪模型，通过构建降维的信号协方差矩阵，并对该矩阵进行矢量化，获得了包含波达方向和功率参数的全部信息。

(2)本发明通过构建的降维的信号协方差矩阵，获得了扩展的虚拟孔径阵列，该阵列可提供有限快拍数的高分辨率波达方向估计和抗非均匀噪声的高鲁棒性。

附图说明

下面结合附图和具体实施例对本发明做进一步详细说明。

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

图2是本发明实施例中分别采用交替网格优化算法(AGO)、基于压缩感知理论的L1奇异值分解算法(L1-SVD)和He算法，在信噪比为0时的目标信源波达方向的估计结果图；

图3是本发明实施例中采用本发明方法在信噪比为0时的目标信源波达方向的估计结果图；

图4是本发明实施例中对网格优化算法(AGO)、基于压缩感知理论的L1奇异值分解算法(L1-SVD)、He算法和本发明方法，在不同信噪比条件下，目标信源波达方向的估计均方根误差与信噪比的关系图。

具体实施方式

下面结合附图对本发明的实施例及效果作进一步详细描述。

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

步骤1，建立雷达的接收信号模型x(t)，确定接收信号的测量矩阵X。

(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列，则t时刻接收到的目标信源的回波信号x(t)，其表达式为：

x(t)＝As(t)+n(t)

其中，A为导向矩阵，s(t)是信号波形矢量，n(t)＝[n₁(t),n₂(t),…,n_M(t)]^T为零均值加性非均匀复数高斯白噪声向量，且n(t)～CN(0,Q)，Q是与n(t)相关的噪声协方差矩阵。

(1.2)导向矩阵A的表达式为：

A＝[a(θ₁),a(θ₂),a(θ_l),…,a(θ_L)]；

其中，L是目标信源数，l＝1,2,...,L；a(θ_l)是M×1的导向矢量，M代表雷达阵元数，θ_l代表第l个入射信号角度，θ_l∈Θ，Θ表示角度搜索范围，其表达式为：

a(θ_l)＝[1,e^-jα,…,e^-j(M-1)α]^T；

其中，α＝2πdsin(θ_l)/λ对应第l个DOA的阵列几何形状，d代表阵元间距，λ代表波长；

信号波形矢量s(t)的表达式为：

s(t)＝[s₁(t),s₂(t),…,s_L(t)]^T；

(1.3)根据接收信号模型x(t)，得到接收信号的测量矩阵X：

其中，J是快拍数，t_z表示第z个快拍时刻，x(t_z)表示t_z快拍时刻天线阵列所接收的M×1维回波信号数据，x₁(t_z)表示t_z快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据，*表示共轭操作。

步骤2，根据接收信号的测量矩阵X，构造基于秩一去噪模型的信号协方差矩阵模型

(2.1)根据接收信号的测量矩阵X，构建基于秩一去噪模型的信号协方差矩阵的第m列其表达式为：

其中，E_m是统计期望；是一个秩一相关矢量；X((1:M)_m,t_z)表示除去第m个阵元的测量矩阵；是第m个阵元的测量矩阵的共轭转置；t_z表示第z个快拍时刻。

(2.2)由M个秩一相关矢量得到基于秩一去噪模型的信号协方差矩阵的组合表达式：

其中，A((1:M)₁,:)表示导向矩阵中第1到M行中的第1行的所有列；表示信号功率矢量；P_l表示第l个信号功率。

(2.3)根据基于秩一去噪模型的信号协方差矩阵的组合表达式，当矩阵中的行序号与列序号不相等时，即x≠y，基于秩一去噪模型的信号协方差矩阵的第x行第y列元素的表达式为：

(2.4)根据基于秩一去噪模型的信号协方差矩阵的第x行第y列元素的表达式，得到基于秩一去噪模型的信号协方差矩阵为：

其中，表示L个信号功率的总和， d代表阵元间距，λ代表波长。

步骤3，对基于秩一去噪模型的信号协方差矩阵模型进行稀疏重构，得到稀疏信号矢量Υ。其具体实施步骤为：

(3.1)根据基于秩一去噪模型的信号协方差矩阵模型构造一个维数为2(M-1)×1的信号矢量Υ'，信号矢量Υ'的第项的表达式为：

其中，表示的第行第列的元素。

(3.2)将信号矢量Υ'进行线性表示为：

Υ'＝B(θ)P；

其中，P表示信号功率矢量，P＝[P₁,P₂,…,P_L]^T；B(θ)表示虚拟流形矩阵，B(θ)＝[b(θ₁),…,b(θ_l),…,b(θ_L)]，第l列导向矢量b(θ_l)为d代表阵元间距，λ代表波长。

(3.3)根据实际回波信号在空域中的稀疏性，采用重构算法对信号矢量Υ'进行稀疏重构，得到稀疏重构信号矢量Υ：

其中， 表示划分的一组字典角度，一般 是具有2(M-1)＜＜N_θ属性的过完备字典，且满足受限等距属性(RIP)，N_θ是字典数；表示稀疏信号功率矢量，其仅有L列非零系数，N_θ＞＞L。

步骤4，根据稀疏重构信号矢量Υ，采用AGO算法估计目标信源的波达方向。

其具体实施步骤为：

(4.1)按照下式，计算导向矩阵A的右逆矩阵：

Φ_a＝A^H(A·A^H)^-1；

其中，Φ_a表示导向矩阵A的右逆矩阵，A表示导向矩阵，H表示共轭转置操作，上标-1表示求逆操作；

(4.2)按照下式，计算导向矩阵A的正交投影矩阵：

Φ＝I-Φ_aA；

其中，Φ表示导向矩阵A的正交投影矩阵，I表示单位矩阵，Φ_a表示导向矩阵A的右逆矩阵，A表示导向矩阵；

(4.3)按照下式，计算初始恢复矢量：

其中，表示初始恢复矢量，Φ_a表示导向矩阵A的右逆矩阵，Υ表示稀疏重构信号；

(4.4)估计信源数：

(4.4a)将动态信源数初始化为1；

(4.4b)按照下式，对第k次内循环中的恢复矢量进行排序操作：

其中，表示恢复矢量取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(·|,'descend')表示降序排列操作；

(4.4c)按照下式，计算第k+1次内循环中的恢复矢量：

其中，表示第k+1次内循环中的恢复矢量，表示第k次内循环中的恢复矢量，u^k表示第k次内循环中的中间辅助矢量，Φ表示导向矩阵A的正交投影矩阵；

(4.4d)按照下式，计算内循环相对误差值：

其中，H₂表示内循环相对误差值，和分别表示内循环次数为k+1和k时的恢复矢量，||·||₂表示取2范数操作；

(4.4e)判断内循环相对误差值H₂是否大于10-³，若是，则执行步骤(4.4b)，否则，执行步骤(4.4f)；

(4.4f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数；

(4.4g)按照下式，计算失配相对误差：

其中，γ^l+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示导向矩阵，u^k表示第k次内循环中的中间辅助矢量，X表示测量矩阵，||·||₂表示取2范数操作；

(4.4h)按照下式，计算外循环相对误差值：

H₁＝|γ^l+1-γ^l|；

其中，H₁表示外循环相对误差值，γ^l+1和γ^l分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作；

(4.4i)判断外循环相对误差值H₁是否大于0.05，若是，则执行步骤(4.4b)，否则，执行步骤(4.4j)；

(4.4j)将外循环结束时的动态信源数值作为信源数的估计值；

(4.5)第一次估计目标到达角：

(4.5a)按照下式，寻找峰值位置矢量：

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.5b)将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量；

(4.5c)取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量中；

(4.6)按照下式，寻找峰值位置矢量：

其中，pV表示恢复矢量取模值后按降序重排的峰值矢量，pI表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.7)按照下式，计算当前代价函数：

其中，F^(M)表示当前代价函数，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，∑(·)表示求和，表示恢复矢量中的第i个元素， 表示不属于符号，表示峰值位置矢量p_I的前个元素，表示信源数的估计值，|·|表示取模值操作；

(4.8)按照下式，计算当前峰值：

其中，表示当前峰值，表示当前峰值对应的位置下标，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，p_I(s)表示峰值位置矢量p_I中的第s个元素；

(4.9)计算左移代价函数：

(4.9a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围；

(4.9b)按照下式，计算左移导向矩阵：

A^(L)＝[a(θ₁),a(θ₂),…,a(θ_i),…]

其中，A^(L)表示左移导向矩阵，M表示天线阵列中的阵元个数，θ_i∈Θ^(L)，∈表示属于符号，Θ(^L)表示左移角度搜索范围，i表示搜索角度θ_i在左移角度搜索范围Θ(^L)中的序号；

(4.9c)按照下式，计算左移导向矩阵的正交投影矩阵：

Φ^(L)＝I-(A^(L))^H(A^(L)(A^(L))^H)^-1A^(L)

其中，Φ^(L)表示左移导向矩阵A^(L)的正交投影矩阵，I表示单位矩阵，A^(L)表示左移导向矩阵，H表示共轭转置操作，-1表示求逆操作；

(4.9d)按照下式，计算左移恢复矢量：

其中，表示左移恢复矢量，u^(L)表示左移辅助矢量，Φ^(L)表示左移导向矩阵A(^L)的正交投影矩阵；

(4.9e)按照下式，寻找左移峰值位置矢量：

其中，表示左移恢复矢量取模值后按降序重排的左移峰值矢量，表示左移恢复矢量元素取模值降序重排后，左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.9f)按照下式，计算左移代价函数：

其中，F^(L)表示左移代价函数，表示左移峰值矢量中的第s个元素，s表示动态信源数，∑(·)表示求和，表示左移恢复矢量中的第i个元素， 表示不属于符号，表示左移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(4.10)计算右移代价函数：

(4.10a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围；

(4.10b)按照下式，计算右移导向矩阵：

A^(R)＝[a(θ₁),a(θ₂),…,a(θ_i),…]

其中，A^(R)为右移导向矩阵，M表示天线阵列中的阵元个数，θ_i∈Θ^(R)，∈表示属于符号，Θ^(R)表示右移角度搜索范围，i表示搜索角度θ_i在右移角度搜索范围Θ^(R)中的序号；

(4.10c)按照下式，计算右移导向矩阵的正交投影矩阵：

Φ^(R)＝I-(A^(R))^H(A^(R)(A^(R))^H)^-1A^(R)

其中，Φ^(R)表示右移导向矩阵A(^R)的正交投影矩阵，I表示单位矩阵，A(^R)表示右移导向矩阵，H表示共轭转置操作，-1表示求逆操作；

(4.10d)按照下式，计算右移恢复矢量：

其中，表示右移恢复矢量，u^(R)表示右移辅助矢量，Φ^(R)表示右移导向矩阵A(^R)的正交投影矩阵；

(4.10e)按照下式，寻找右移峰值位置矢量：

其中，表示右移恢复矢量取模值后按降序重排的右移峰值矢量，表示右移峰值矢量元素取模值降序重排后，右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.10f)按照下式，计算右移代价函数：

其中，F^(R)表示右移代价函数，表示右移峰值矢量中的第s个元素，s表示动态信源数，∑(·)表示求和，表示右移恢复矢量中的第i个元素， 表示不属于符号，表示右移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(4.11)更新当前栅格参数：

(4.11a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(4.11b)，否则，执行步骤(4.11e)；

(4.11b)将步骤(4.7)的当前代价函数的值更新为左移代价函数的值；

(4.11c)将恢复矢量的元素值更新为左移恢复矢量的元素值；

(4.11d)按照下式，更新目标到达角矢量：

其中，表示更新后的目标到达角矢量，Θ^(L)表示左移角度搜索范围，表示左移峰值位置矢量的第s个元素，s表示动态信源数；

(4.11e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(4.11f)，否则，执行步骤(4.12)；

(4.11f)将当前代价函数的值更新为右移代价函数的值；

(4.11g)将恢复矢量的元素值更新为右移恢复矢量的元素值；

(4.11h)按照下式，更新目标到达角估计矢量：

其中，表示更新后的目标到达角估计矢量，Θ^(R)表示右移角度搜索范围，表示右移峰值位置矢量的第s个元素，s表示动态信源数；

(4.12)判断可调栅格步长△是否大于栅格优化精度ξ_θ，若是，则执行步骤(4.9)，否则，执行步骤(4.13)；

(4.13)按照下式，更新目标到达角估计矢量：

其中，表示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，表示更新后的目标到达角估计矢量，表示当前峰值对应的位置下标；

(4.14)判断动态信源数是否小于信源数，若是，则执行步骤(4.6)，否则，执行步骤(4.15)；

(4.15)获得第二次目标到达角估计矢量：

将步骤(4.13)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。

仿真实验

本发明的效果可通过以下仿真实验进一步说明。

(1)仿真参数：

本发明采用M＝24个阵元，且阵元间距为半波长(d＝λ2)的均匀线性阵列。采用在-90°至90°范围内具有1°间隔的一般网格。阈值τ和误差常数ξ分别为5M和5×10^-4。假定非均匀噪声协方差矩阵为Q＝diag{Q_n,3*Q_n,1.5*Q_n}，其中Q_n＝[2.0,10,2.5,5.0,0.5,1.5,3.0,5.0]。采用三个不相关的信号，其中DOA为{-5°,1°,4°}，信噪比SNR分别为-6dB、-4dB、-2dB、0dB、2dB、4dB、6dB。

参数设置如表1：

表1系统仿真参数

(2)仿真内容：

仿真1，在上述仿真参数下，分别采用AGO算法、L1-SVD算法和He算法在信噪比为0的条件下，对目标的回波信号进行波达方向估计，结果如图2所示。

由图2可以看出，由于存在非均匀噪声，L1-SVD算法无法正常工作；AGO算法和He算法可以确定三个不相关信号。

仿真2，在上述仿真参数下，采用本发明方法在信噪比为0的条件下，对回波信号进行波达方向估计，结果如图3所示。

由图3可以看出，本发明方法能够确定三个不相关信号；比较图2和图3可以看出，本发明方法得到的目标波达方向估计结果图的目标处的峰较尖锐，说明本发明方法的估计结果精度最高。且在旁瓣较低的情况下，本发明方法稍微优于AGO算法。

在上述仿真参数下，分别采用AGO算法、L1-SVD算法、He算法和本发明方法在信噪比为分别为-6dB、-4dB、-2dB、0dB、2dB、4dB、6dB的条件下，对回波信号进行波达方向估计，并采用均方根误差RMSE对估计结果进行评价：

其中，K是蒙特卡洛实验次数，表示第k次蒙特卡洛实验中第l个目标的估计角度，θ_l表示第l个目标的真实角度。

以上试验结果如图4所示，由图4可以看出，随着信噪比SNR的增大，四种方法的估计均方根误差RMSE均呈现下降趋势，且本发明方法与He算法的估计性能相当，且均优于AGO算法和L1-SVD算法，但本发明方法在信噪比为2时的性能优于AGO算法，说明本发明方法在信噪比较高时，具有更优的目标波达方向估计结果。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

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

Claims (7)

1.基于秩一去噪模型的波达方向估计方法，其特征在于，包括以下步骤：

步骤1，建立雷达的接收信号模型x(t)，确定接收信号的测量矩阵X；

步骤2，根据接收信号的测量矩阵X，构造基于秩一去噪模型的信号协方差矩阵

步骤3，对基于秩一去噪模型的信号协方差矩阵进行稀疏重构，得到稀疏信号矢量Υ；

步骤4，根据稀疏重构信号矢量Υ，采用交替网格优化算法，估计目标信源的波达方向。

2.根据权利要求1所述的基于秩一去噪模型的波达方向估计方法，其特征在于，步骤1的具体步骤为：

(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列，则t时刻接收到的目标信源的回波信号，即t时刻的接收信号x(t)的表达式为：

x(t)＝As(t)+n(t)；

其中，A为导向矩阵，s(t)是信号波形矢量，n(t)＝[n₁(t)，n₂(t)，…，n_M(t)]^T为零均值加性非均匀复数高斯白噪声向量，且n(t)～CN(0，Q)，Q是与n(t)相关的噪声协方差矩阵；

(1.2)根据接收信号模型x(t)，得到接收信号的测量矩阵X：

其中，J是快拍数，t_z表示第z个快拍时刻，x(t_z)表示t_z快拍时刻天线阵列所接收的M×1维回波信号数据，x₁(t_z)表示t_z快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据，*表示共轭操作。

3.根据权利要求2所述的基于秩一去噪模型的波达方向估计方法，其特征在于，所述导向矩阵A的表达式为：

A＝[a(θ₁)，a(θ₂)，a(θ_l)，…，a(θ_L)]；

其中，L是目标信源数，l＝1，2，...，L；M代表雷达阵元数，θ_l代表第l个入射信号角度，θ_l∈Θ，Θ表示角度搜索范围；a(θ_l)是M×1的导向矢量，其表达式为：

a(θ_l)＝[1，e^-jα，…，e^-j(M-1)α]^T；

其中，α表示对应第l个DOA的阵列几何形状，α＝2πdsin(θ_l)/λ，d代表阵元间距，λ代表波长；

所述信号波形矢量s(t)的表达式为：

s(t)＝[s₁(t)，s₂(t)，…，s_L(t)]^T；

其中，s_L(t)表示第L个信号波形。

4.根据权利要求1所述的基于秩一去噪模型的波达方向估计方法，其特征在于，步骤2的具体步骤为：

(2.1)根据接收信号的测量矩阵X，构建基于秩一去噪模型的信号协方差矩阵的第m列其表达式为：

其中，E_m是统计期望；是一个秩一相关矢量；X((1：M)_m，t_z)表示除去第m个阵元的测量矩阵；是第m个阵元的测量矩阵的共轭转置；t_z表示第z个快拍时刻，M表示雷达的阵元总数；

(2.2)由M个秩一相关矢量得到基于秩一去噪模型的信号协方差矩阵的组合表达式：

其中，A((1∶M)₁，：)表示导向矩阵中第1到M行中的第1行的所有列；表示信号功率矢量；P_l表示第l个信号功率；

(2.3)根据基于秩一去噪模型的信号协方差矩阵的组合表达式，当矩阵中的行序号与列序号不相等时，即x≠y，基于秩一去噪模型的信号协方差矩阵的第x行第y列元素的表达式为：

其中，A表示导向矩阵，上标H表示共轭转置，θ_l代表第l个入射信号角度，θ_l∈Θ，Θ表示角度搜索范围；d代表阵元间距，λ代表波长；

(2.4)根据基于秩一去噪模型的信号协方差矩阵的第x行第y列元素的表达式，得到基于秩一去噪模型的信号协方差矩阵为：

其中，表示L个信号功率的总和，j表示虚数单位。

5.根据权利要求1所述的基于秩一去噪模型的波达方向估计方法，其特征在于，步骤3的具体步骤为：

(3.1)根据基于秩一去噪模型的信号协方差矩阵构造一个2(M-1)×1的信号矢量Υ′；

其中，M表示雷达的阵元总数；

(3.2)将信号矢量Υ′进行线性表示，其线性表达式为：

Υ′＝B(θ)P；

其中，P表示信号功率矢量，P＝[P₁，P₂，…，P_L]^T；B(θ)表示虚拟流形矩阵，B(θ)＝[b(θ₁)，…，b(θ_l)，…，b(θ_L)]，b(θ_l)表示第l列导向矢量，且θ_l代表第l个入射信号角度，θ_l∈Θ，Θ表示角度搜索范围；d代表阵元间距，λ代表波长；j表示虚数单位；

(3.3)根据信号矢量Υ′的线性表达式，对信号矢量Υ′的进行稀疏重构，得到稀疏重构信号矢量Υ：

其中， 表示划分的一组字典角度；是具有2(M-1)＜＜N_θ属性的过完备字典，且满足受限等距属性；表示字典角度为时的稀疏导向矢量，即第1列稀疏导向矢量，N_θ是字典数；表示稀疏信号功率矢量，且 表示第1个稀疏信号功率。

6.根据权利要求5所述的基于秩一去噪模型的波达方向估计方法，其特征在于，所述信号矢量Υ′的第项的表达式为：

其中，表示的第行第列的元素。

7.根据权利要求5所述的基于秩一去噪模型的波达方向估计方法，其特征在于，步骤4的具体步骤为：

(4.1)按照下式，计算导向矩阵A的右逆矩阵：

Φ_a＝A^H(A·A^H)^-1；

其中，Φ_a表示导向矩阵A的右逆矩阵，A表示导向矩阵，H表示共轭转置操作，上标-1表示求逆操作；

(4.2)按照下式，计算导向矩阵A的正交投影矩阵：

Φ＝I-Φ_aA；

其中，Φ表示导向矩阵A的正交投影矩阵，I表示单位矩阵，Φ_a表示导向矩阵A的右逆矩阵，A表示导向矩阵；

(4.3)根据稀疏重构信号矢量Υ，按照下式，计算初始恢复矢量：

其中，表示初始恢复矢量，Φ_a表示导向矩阵A的右逆矩阵；

(4.4)估计信源数：

(4.4a)将动态信源数初始化为1；

(4.4b)按照下式，对第k次内循环中的恢复矢量进行排序操作：

其中，表示恢复矢量取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(|·|，′descend′)表示降序排列操作；

(4.4c)按照下式，计算第k+1次内循环中的恢复矢量：

其中，表示第k+1次内循环中的恢复矢量，表示第k次内循环中的恢复矢量，u^k表示第k次内循环中的中间辅助矢量，Φ表示导向矩阵A的正交投影矩阵；

(4.4d)按照下式，计算内循环相对误差值：

其中，H₂表示内循环相对误差值，和分别表示内循环次数为k+1和k时的恢复矢量，||·||₂表示取2范数操作；

(4.4e)判断内循环相对误差值H₂是否大于10^-3，若是，则执行步骤(4.4b)，否则，执行步骤(4.4f)；

(4.4f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数；

(4.4g)按照下式，计算失配相对误差：

其中，γ^l+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示导向矩阵，u^k表示第k次内循环中的中间辅助矢量，X表示测量矩阵，||·||₂表示取2范数操作；

(4.4h)按照下式，计算外循环相对误差值：

H₁＝|γ^l+1-γ^l|；

其中，H₁表示外循环相对误差值，γ^l+1和γ^l分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作；

(4.4i)判断外循环相对误差值H₁是否大于0.05，若是，则执行步骤(4.4b)，否则，执行步骤(4.4j)；

(4.4j)将外循环结束时的动态信源数值作为信源数的估计值；

(4.5)第一次估计目标到达角：

(4.5a)按照下式，寻找峰值位置矢量：

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.5b)将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量；

(4.5c)取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量中；

(4.6)按照下式，寻找峰值位置矢量：

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.7)按照下式，计算当前代价函数：

其中，F^(M)表示当前代价函数，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，∑(·)表示求和，表示恢复矢量中的第i个元素， 表示不属于符号，表示峰值位置矢量p_I的前个元素，表示信源数的估计值，|·|表示取模值操作；

(4.8)按照下式，计算当前峰值：

其中，表示当前峰值，表示当前峰值对应的位置下标，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，p_I(s)表示峰值位置矢量p_I中的第s个元素；

(4.9)计算左移代价函数：

(4.9a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围；

(4.9b)按照下式，计算左移导向矩阵：

A^(L)＝[a(θ₁)，a(θ₂)，…，a(θ_i)，…]

其中，A^(L)表示左移导向矩阵，M表示天线阵列中的阵元个数，θ_i∈Θ^(L)，∈表示属于符号，Θ^(L)表示左移角度搜索范围，i表示搜索角度θ_i在左移角度搜索范围Θ^(L)中的序号；

(4.9c)按照下式，计算左移导向矩阵的正交投影矩阵：

Φ^(L)＝I-(A^(L))^H(A^(L)(A^(L))^H)^-1A^(L)

其中，Φ^(L)表示左移导向矩阵A^(L)的正交投影矩阵，I表示单位矩阵，A^(L)表示左移导向矩阵，H表示共轭转置操作，-1表示求逆操作；

(4.9d)按照下式，计算左移恢复矢量：

其中，表示左移恢复矢量，u^(L)表示左移辅助矢量，Φ^(L)表示左移导向矩阵A^(L)的正交投影矩阵；

(4.9e)按照下式，寻找左移峰值位置矢量：

其中，表示左移恢复矢量取模值后按降序重排的左移峰值矢量，表示左移恢复矢量元素取模值降序重排后，左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.9f)按照下式，计算左移代价函数：

其中，F^(L)表示左移代价函数，表示左移峰值矢量中的第s个元素，s表示动态信源数，∑(·)表示求和，表示左移恢复矢量中的第i个元素， 表示不属于符号，表示左移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(4.10)计算右移代价函数：

(4.10a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围；

(4.10b)按照下式，计算右移导向矩阵：

A^(R)＝[a(θ₁)，a(θ₂)，…，a(θ_i)，…]；

其中，A^(R)为右移导向矩阵，M表示天线阵列中的阵元个数，θ_i∈Θ^(R)，∈表示属于符号，Θ^(R)表示右移角度搜索范围，i表示搜索角度θ_i在右移角度搜索范围Θ^(R)中的序号；

(4.10c)按照下式，计算右移导向矩阵的正交投影矩阵：

Φ^(R)＝I-(A^(R))^H(A^(R)(A^(R))^H)^-1A^(R)

其中，Φ^(R)表示右移导向矩阵A^(R)的正交投影矩阵，I表示单位矩阵，A^(R)表示右移导向矩阵，H表示共轭转置操作，-1表示求逆操作；

(4.10d)按照下式，计算右移恢复矢量：

其中，表示右移恢复矢量，u^(R)表示右移辅助矢量，Φ^(R)表示右移导向矩阵A^(R)的正交投影矩阵；

(4.10e)按照下式，寻找右移峰值位置矢量：

其中，表示右移恢复矢量取模值后按降序重排的右移峰值矢量，表示右移峰值矢量元素取模值降序重排后，右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(4.10f)按照下式，计算右移代价函数：

其中，F^(R)表示右移代价函数，表示右移峰值矢量中的第s个元素，s表示动态信源数，∑(·)表示求和，表示右移恢复矢量中的第i个元素， 表示不属于符号，表示右移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(4.11)更新当前栅格参数：

(4.11a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(4.11b)，否则，执行步骤(4.11e)；

(4.11b)将步骤(4.7)的当前代价函数的值更新为左移代价函数的值；

(4.11c)将恢复矢量的元素值更新为左移恢复矢量的元素值；

(4.11d)按照下式，更新目标到达角矢量：

其中，表示更新后的目标到达角矢量，Θ^(L)表示左移角度搜索范围，表示左移峰值位置矢量的第s个元素，s表示动态信源数；

(4.11e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(4.11f)，否则，执行步骤(4.12)；

(4.11f)将当前代价函数的值更新为右移代价函数的值；

(4.11g)将恢复矢量的元素值更新为右移恢复矢量的元素值；

(4.11h)按照下式，更新目标到达角估计矢量：

其中，表示更新后的目标到达角估计矢量，Θ^(R)表示右移角度搜索范围，表示右移峰值位置矢量的第s个元素，s表示动态信源数；

(4.12)判断可调栅格步长Δ是否大于栅格优化精度ξ_θ，若是，则执行步骤(4.9)，否则，执行步骤(4.13)；

(4.13)按照下式，更新目标到达角估计矢量：

其中，表示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，表示更新后的目标到达角估计矢量，表示当前峰值对应的位置下标；

(4.14)判断动态信源数是否小于信源数，若是，则执行步骤(4.6)，否则，执行步骤(4.15)；

(4.15)获得第二次目标到达角估计矢量：

将步骤(4.13)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。