CN113376569B - 基于最大似然的嵌套阵稀疏表示波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于最大似然的嵌套阵稀疏表示波达方向估计方法，首先根据目标到嵌套阵的波达方向以及嵌套阵的布阵结构，计算接收阵列的回波信号、进而计算其协方差矩阵；根据稀疏表示理论，在角度维划分整个空域，得到角度集合；矢量化协方差矩阵，在角度集合上稀疏展开，得到嵌套阵波达方向估计稀疏模型；构造块对角矩阵，剔除稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型；再结合协方差矩阵计算噪声白化矩阵，结合去噪后嵌套阵波达方向估计稀疏模型，计算噪声白化后的数据模型，进而计算栅格最大似然估计，最后得到波达方向最大似然估计值。本发明以实现嵌套阵在低信噪比和少快拍下的波达方向估计性能。

Description

基于最大似然的嵌套阵稀疏表示波达方向估计方法

技术领域

本发明属于雷达技术领域，具体涉及一种基于最大似然的嵌套阵稀疏表示波达方向估计方法。

背景技术

波达方向估计技术是实现目标定位不可或缺的技术手段，通过波达方向估计结果可得目标在空域的方位信息。传统典型的方法有多重信号分类方法，旋转不变子空间方法，这两种方法突破瑞利限，可实现目标的超分辨，但在低信噪比或者少快拍环境下性能严重下降。

最大似然估计假设信号源为一具有已知分布的随机过程，利用已知的样本结果信息，反推最具有可能导致这些样本结果出现的待估计参数值。在低信噪比或者相干信号源情况下仍然可以获得较好的估计性能，但是求解比较困难，需要在多维参数空间上寻优，这将导致算法的运算量很大，并且算法的性能依赖于初始值的选取。

嵌套阵是由多个均匀线阵嵌套而成，通过非均匀布阵减少空间采样的冗余性，可提高线性阵列的自由度，并且阵元位置的闭式解和嵌套阵所能提高的自由度均可通过阵元总数计算。

稀疏表示理论的基本思想是用冗余函数集合组成的原子字典代替基函数集合，进而可以将信号表示为原子字典中少数原子列向量的线性组合。当空域分布少量点目标时，对整个空域角度而言目标具有稀疏性，因此可以将稀疏表示理论应用到波达方向估计中。但需假设目标正好落在划分的角度栅格上，这将不可避免的存在模型失配问题。

针对稀疏表示波达方向估计方法中模型失配问题，现有方法主要分为两大类，一类是多次重构法，算法的计算量随着栅格个数增加按指数级增长，另一类是贝叶斯估计法，虽然降低了字典矩阵栅格失配问题，但仍然存在失配量。

发明内容

本发明的目的是提供一种基于最大似然的嵌套阵稀疏表示波达方向估计方法，以实现嵌套阵在低信噪比和少快拍下的波达方向估计性能。

本发明所采用的技术方案是，基于最大似然的嵌套阵稀疏表示波达方向估计方法，具体按照以下步骤实施：

步骤1、根据目标到嵌套阵雷达系统接收阵列的波达方向以及嵌套阵的布阵结构，计算接收阵列的回波信号；

步骤2、根据步骤1得到的接收阵列的回波信号，计算嵌套阵接收数据的协方差矩阵；

步骤3、根据稀疏表示理论，在角度维划分整个空域，得到角度集合；矢量化步骤2得到的嵌套阵接收数据协方差矩阵，并在角度集合上稀疏展开，得到嵌套阵波达方向估计稀疏模型；

步骤4、构造块对角矩阵，剔除步骤3稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型；再结合步骤2得到的协方差矩阵计算噪声白化矩阵，结合得到的去噪后嵌套阵波达方向估计稀疏模型，计算噪声白化后的稀疏模型；

步骤5、根据步骤4建立的噪声白化后的稀疏模型，计算目标波达方向。

本发明的特点还在于，

步骤1中接收阵列的回波信号y(t)计算如下：

y(t)＝A(θ)s(t)+n(t)，

其中，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T表示信号矢量，[·]^T为转置运算，K表示目标个数，n(t)表示通道噪声矢量，假设服从复高斯分布，即 表示均值为μ，协方差矩阵为Σ的复高斯分布，σ²表示噪声功率，I_M×M表示维数为M×M的单位矩阵，M表示阵元个数，A(θ)为阵列流形矩阵，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_k),…,a(θ_K)]，a(θ_k)表示阵列导向矢量，θ_k表示第k个目标的来波方向，k＝1,2,…,K，/>(·)_m表示矢量的第m个元素，D_m表示嵌套阵雷达系统第m个阵元相对参考阵元的位置信息，m＝1,2,…,M，λ表示电磁波的波长，t表示采样率归一化的时间，t＝1,2,…,L，L为总快拍数。

步骤2协方差矩阵的计算具体如下：

其中，(·)^H为共轭转置运算。

步骤3具体如下：

步骤3.1、根据稀疏表示理论，在角度维划分整个空域，得到角度集合Θ：

Θ＝{θ₁,θ₂,…,θ_n,…,θ_N}，

其中，N表示空域角度划分个数，θ_n表示第n个划分的角度，n＝1,2,…,N；

步骤3.2、矢量化步骤2得到的嵌套阵接收数据协方差矩阵

其中，vec(·)表示矢量化运算；

步骤3.3、将步骤3.2得到的矢量化后的y_v在步骤3.1得到的角度集合Θ上稀疏展开：

y_v＝Ψ(Θ)p_Θ+σ²1+Δy_v，

其中， 表示Kronecker积，(·)*表示共轭运算，p_Θ为稀疏向量，非零位置对应的角度信息即为目标的来波方向，/>e_m表示除第m个元素为1外，其它元素均为0的单位矢量，Δy_v表示嵌套阵理论协方差矩阵和实际协方差矩阵差值矢量化后的向量，嵌套阵波达方向估计稀疏模型即为y_v＝Ψ(Θ)p_Θ+σ²1+Δy_v。

步骤4具体如下：

步骤4.1、构造块对角矩阵J：

其中，J_m＝[e₁,…,e_m-1,e_m+1,…,e_M]，m＝2,…,M-1，J₁＝[e₂,…,e_M]，J_M＝[e₁,…,e_M-1]；

步骤4.2、根据步骤2得到的协方差矩阵和步骤4.1得到块对角矩阵J，计算噪声白化矩阵W：

步骤4.3、建立去噪后嵌套阵波达方向估计稀疏模型：

利用步骤4.1得到的块对角矩阵J，剔除步骤3稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型：

y_J＝Jy_v＝JΨ(Θ)p_Θ+JΔy_v；

步骤4.4、根据步骤4.2得到的噪声白化矩阵W及步骤4.3得到的稀疏模型，建立噪声白化后的波达方向估计稀疏模型y_w：

y_w＝W^-1/2y_J＝W^-1/2JΨ(Θ)p_Θ+ε＝Φ(Θ)p_Θ+ε，

其中，Φ(Θ)＝W^-1/2JΨ(Θ)，服从复高斯白噪声分布。

步骤5具体如下：

步骤5.1、假设稀疏向量p_Θ服从复高斯分布，即

其中，P(·|·)表示条件概率，Γ＝diag(γ₁，γ₂，…，γ_N)，diag(·)表示对角操作运算；

步骤5.2、根据步骤4.4中Φ(Θ)、y_w以及步骤5.1中Γ，分别计算稀疏向量p_Θ的均值和协方差矩阵/>

步骤5.3、根据步骤5.2得到的稀疏向量p_Θ的均值和协方差矩阵/>计算栅格最大似然估计/>

其中，ζ是一个极小正数，(·)^q表示第q次迭代，(·)_n,n表示矩阵第n行第n列的元素；

步骤5.4、根据步骤5.3得到的所有的栅格最大似然估计值形成空间谱，由谱峰位置计算协方差矩阵Σ_-k：

Σ_-k＝Φ(Θ_-k)diag(γ_-k)Φ(Θ_-k)^H+I_{M(M-1)×M(M-1)}，

其中，表示从集合Θ中删除空间谱谱峰所对应的第k个信号源对应的栅格/> 表示从步骤5.3得到的栅格最大似然估计结果中删除第k个信号源对应的/>的估计值；

步骤5.5、计算参数

根据步骤5.4得到的协方差矩阵Σ-_k和步骤5.5得到的参数计算目标波达方向最大似然估计θ_k：

其中，Re{·}表示取实部运算，/>表示第k个信号源对应的栅格/>左右领域的角度集合，argmax[·]表示取函数最大值的变量值。

本发明的有益效果是，一种基于最大似然的嵌套阵稀疏表示波达方向估计方法，相比于均匀阵列，本发明由于采用矢量化嵌套阵协方差矩阵后的数据建立波达方向稀疏重构模型，扩大了阵列的虚拟孔径，故可有效的提高算法的估计精度和分辨性能；相比于现有方法，本发明由于采用了最大似然方法，建立了波达方向估计的最大似然模型，可通过一维搜索获得目标的来波方向，避免多次重构带来计算量骤增的问题及泰勒级数近似模型失配问题。

附图说明

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

图2是两层嵌套阵示意图；

图3是用本发明方法对空域9个不相关信号空间谱的仿真结果图；

图4是用现有方法和用本发明方法对空中角度间隔较小的两目标的成功分辨概率随信噪比的变化曲线的仿真结果图；

图5是用现有方法和用本发明方法对空中角度间隔较大的两目标的均方根误差随快拍数的变化曲线的仿真结果图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明基于最大似然的嵌套阵稀疏表示波达方向估计方法，流程图如图1所示，具体按照以下步骤实施：

步骤1、根据目标到嵌套阵雷达系统接收阵列的波达方向以及嵌套阵的布阵结构，计算接收阵列的回波信号。

步骤1中接收阵列的回波信号y(t)计算如下：

y(t)＝A(θ)s(t)+n(t)，

其中，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T表示信号矢量，[·]^T为转置运算，K表示目标个数，n(t)表示通道噪声矢量，假设服从复高斯分布，即 表示均值为μ，协方差矩阵为Σ的复高斯分布，σ²表示噪声功率，I_M×M表示维数为M×M的单位矩阵，M表示阵元个数，A(θ)为阵列流形矩阵，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_k),…a(θ_K)]，a(θ_k)表示阵列导向矢量，θ_k表示第k个目标的来波方向，k＝1,2,…,K，/>(·)_m表示矢量的第m个元素，D_m表示嵌套阵雷达系统第m个阵元相对参考阵元的位置信息，m＝1,2,…,M，λ表示电磁波的波长，t表示采样率归一化的时间，t＝1,2,…,L，L为总快拍数。

步骤2、根据步骤1得到的接收阵列的回波信号，计算嵌套阵接收数据的协方差矩阵

步骤2协方差矩阵的计算具体如下：

其中，(·)^H为共轭转置运算。

步骤3、根据稀疏表示理论，在角度维划分整个空域，得到角度集合；矢量化步骤2得到的嵌套阵接收数据协方差矩阵，并在角度集合上稀疏展开，得到嵌套阵波达方向估计稀疏模型。

步骤3具体如下：

步骤3.1、根据稀疏表示理论，在角度维划分整个空域，得到角度集合Θ：

Θ＝{θ₁,θ₂,…,θ_n,…,θ_N}，

其中，N表示空域角度划分个数，θ_n表示第n个划分的角度，n＝1,2,…,N；

步骤3.2、矢量化步骤2得到的嵌套阵接收数据协方差矩阵

其中，vec(·)表示矢量化运算；

步骤3.3、将步骤3.2得到的矢量化后的y_v在步骤3.1得到的角度集合Θ上稀疏展开：

y_v＝Ψ(Θ)p_Θ+σ²1+Δy_v，

其中， 表示Kronecker积，(·)*表示共轭运算，p_Θ为稀疏向量，非零位置对应的角度信息即为目标的来波方向，/>e_m表示除第m个元素为1外，其它元素均为0的单位矢量，Δy_v表示嵌套阵理论协方差矩阵和实际协方差矩阵差值矢量化后的向量，嵌套阵波达方向估计稀疏模型即为y_v＝Ψ(Θ)p_Θ+σ²1+Δy_v。

步骤4、构造块对角矩阵，剔除步骤3稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型；再结合步骤2得到的协方差矩阵计算噪声白化矩阵，结合得到的去噪后嵌套阵波达方向估计稀疏模型，计算噪声白化后的稀疏模型。

步骤4具体如下：

步骤4.1、构造块对角矩阵J：

其中，J_m＝[e₁,…,e_m-1,e_m+1,…,e_M]，m＝2,…,M-1，J₁＝[e₂,…,e_M]，J_M＝[e₁,…,e_M-1]；

步骤4.2、根据步骤2得到的协方差矩阵和步骤4.1得到块对角矩阵J，计算噪声白化矩阵W：

步骤4.3、建立去噪后嵌套阵波达方向估计稀疏模型。

利用步骤4.1得到的块对角矩阵J，剔除步骤3稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型：

y_J＝Jy_v＝JΨ(Θ)p_Θ+JΔy_v；

步骤4.4、根据步骤4.2得到的噪声白化矩阵W及步骤4.3得到的稀疏模型，建立噪声白化后的波达方向估计稀疏模型y_w：

y_w＝W^-1/2y_J＝W^-1/2JΨ(Θ)p_Θ+ε＝Φ(Θ)p_Θ+ε，

其中，Φ(Θ)＝W^-1/2JΨ(Θ)，服从复高斯白噪声分布。

步骤5、根据步骤4建立的噪声白化后的稀疏模型，计算目标波达方向。

步骤5具体如下：

步骤5.1、假设稀疏向量p_Θ服从复高斯分布，即

其中，P(·|·)表示条件概率，Γ＝diag(γ₁，γ₂，…，γ_N)，diag(·)表示对角操作运算；

步骤5.2、根据步骤4.4中Φ(Θ)、y_w以及步骤5.1中Γ，分别计算稀疏向量p_Θ的均值和协方差矩阵/>

步骤5.3、根据步骤5.2得到的稀疏向量p_Θ的均值和协方差矩阵/>计算栅格最大似然估计/>

其中，ζ是一个极小正数，(·)^q表示第q次迭代，(·)_n,n表示矩阵第n行第n列的元素；

步骤5.4、根据步骤5.3得到的所有的栅格最大似然估计值形成空间谱，由谱峰位置计算协方差矩阵Σ_-k：

Σ_-k＝Φ(Θ_-k)diag(γ_-k)Φ(Θ_-k)^H+I_{M(M-1)×M(M-1)}，

其中，表示从集合Θ中删除空间谱谱峰所对应的第k个信号源对应的栅格/> 表示从步骤5.3得到的栅格最大似然估计结果中删除第k个信号源对应的/>的估计值；

步骤5.5、计算参数

步骤5.6、根据步骤5.4得到的协方差矩阵Σ_-k和步骤5.5得到的参数

计算目标波达方向最大似然估计θ_k：

其中，Re{·}表示取实部运算，/>表示第k个信号源对应的栅格/>左右领域的角度集合，argmax[·]表示取函数最大值的变量值。

本发明基于最大似然的嵌套阵稀疏表示波达方向估计方法，可用于嵌套阵雷达系统，在目标满足空域稀疏性的条件下，通过最大似然方法，建立稀疏栅格迭代过程及角度细估计的表达式，提高嵌套阵雷达系统的波达方向估计性能。

本发明对目标角度信息的估计性能可通过以下仿真进一步验证。

1.实验场景：

采用如图2所示的两层嵌套阵，阵元总数M＝6，阵元位置集合为{0,d,2d,3d,7d,11d}，d等于半波长，空域角度划分间隔为1°，划分后的角度集合为{-90°:1°:90°}。的初始值采用最小二乘估计得到，即/>p0＝(Φ(Θ))⁺y，(·)⁺表示广义逆。迭代终止条件为，达到最大迭代次数或两次迭代更新满足||γ^q+1-γ^q||₂/||γ^q||₂≤ι，其中迭代最大次数为2000次，γ＝[γ₁，γ₂，…，γ_N]^T，ι＝10-⁴。

2.实验内容与分析

实验一：用本发明方法对空间远场9个不相关的等功率信号源进行波达方向估计，得到角度估计的空间谱，如图3所示，其中，“o”表示真实目标的波达方向。

由图3可以看出，本发明方法能成功分辨来自9个不同方向的等功率信号，适用于目标个数大于阵元个数的场景。

实验二：改变信噪比，对每一个信噪比，均利用现有方法和本发明方法进行200次Monte Carlo仿真实验，分别统计出现有方法和本发明方法对空中角度间隔较小的两目标的成功分辨概率随信噪比的变化曲线，如图4所示。

由图4可以看出，随着信噪比增大，现有方法和本发明方法对空中角度间隔较小的两目标的成功分辨概率逐渐增大到100％，但本发明方法在每个信噪比处成功分辨的概率均大于等于现有方法，说明本发明方法的分辨性能优于现有方法。

实验三：改变快拍数，对每一个快拍数，均利用现有方法和本发明方法进行200次Monte Carlo仿真实验，分别统计出现有方法和本发明方法对空中角度间隔较大的两目标的均方根误差随快拍数的变化曲线，如图5所示。

由图5可以看出，随着快拍数的增大，现有方法和本发明方法对空中角度间隔较大的两目标的均方根误差逐渐减小，但本发明方法在每个快拍数处的均方根误差均小于等于现有方法，尤其在快拍数小于等于40时，本发明方法的均方根误差明显低于现有方法，说明本发明方法的估计性能优于现有方法。

综上所述，本发明方法能够对空域目标进行有效的波达方向估计，提高了目标的分辨和估计性能。

Claims (6)

1.基于最大似然的嵌套阵稀疏表示波达方向估计方法，其特征在于，具体按照以下步骤实施：

步骤1、根据目标到嵌套阵雷达系统接收阵列的波达方向以及嵌套阵的布阵结构，计算接收阵列的回波信号；

步骤2、根据步骤1得到的接收阵列的回波信号，计算嵌套阵接收数据的协方差矩阵；

步骤3、根据稀疏表示理论，在角度维划分整个空域，得到角度集合；矢量化步骤2得到的嵌套阵接收数据协方差矩阵，并在角度集合上稀疏展开，得到嵌套阵波达方向估计稀疏模型；

步骤4、构造块对角矩阵，剔除步骤3稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型；再结合步骤2得到的协方差矩阵计算噪声白化矩阵，结合得到的去噪后嵌套阵波达方向估计稀疏模型，计算噪声白化后的稀疏模型；

步骤5、根据步骤4建立的噪声白化后的稀疏模型，计算目标波达方向。

2.根据权利要求1所述的基于最大似然的嵌套阵稀疏表示波达方向估计方法，其特征在于，所述步骤1中接收阵列的回波信号y(t)计算如下：

y(t)＝A(θ)s(t)+n(t)，

其中，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T表示信号矢量，[·]^T为转置运算，K表示目标个数，n(t)表示通道噪声矢量，假设服从复高斯分布，即 表示均值为μ，协方差矩阵为Σ的复高斯分布，σ²表示噪声功率，I_M×M表示维数为M×M的单位矩阵，M表示阵元个数，A(θ)为阵列流形矩阵，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_k),…,a(θ_K)]，a(θ_k)表示阵列导向矢量，θ_k表示第k个目标的来波方向，k＝1,2,…,K，/>(·)_m表示矢量的第m个元素，D_m表示嵌套阵雷达系统第m个阵元相对参考阵元的位置信息，m＝1,2,…,M，λ表示电磁波的波长，t表示采样率归一化的时间，t＝1,2,…,L，L为总快拍数。

3.根据权利要求2所述的基于最大似然的嵌套阵稀疏表示波达方向估计方法，其特征在于，所述步骤2协方差矩阵的计算具体如下：

其中，(·)^H为共轭转置运算。

4.根据权利要求3所述的基于最大似然的嵌套阵稀疏表示波达方向估计方法，其特征在于，所述步骤3具体如下：

步骤3.1、根据稀疏表示理论，在角度维划分整个空域，得到角度集合Θ：

Θ＝{θ₁,θ₂,…,θ_n,…,θ_N}，

其中，N表示空域角度划分个数，θ_n表示第n个划分的角度，n＝1,2,…,N；

步骤3.2、矢量化步骤2得到的嵌套阵接收数据协方差矩阵

其中，vec(·)表示矢量化运算；

步骤3.3、将步骤3.2得到的矢量化后的y_v在步骤3.1得到的角度集合Θ上稀疏展开：

y_v＝Ψ(Θ)p_Θ+σ²1+Δy_v，

其中， 表示Kronecker积，(·)*表示共轭运算，p_Θ为稀疏向量，非零位置对应的角度信息即为目标的来波方向，/>e_m表示除第m个元素为1外，其它元素均为0的单位矢量，Δy_v表示嵌套阵理论协方差矩阵和实际协方差矩阵差值矢量化后的向量，嵌套阵波达方向估计稀疏模型即为y_v＝Ψ(Θ)p_Θ+σ²1+Δy_v。

5.根据权利要求4所述的基于最大似然的嵌套阵稀疏表示波达方向估计方法，其特征在于，所述步骤4具体如下：

步骤4.1、构造块对角矩阵J：

其中，J_m＝[e₁,…,e_m-1,e_m+1,…,e_M]，m＝2,…,M-1，J₁＝[e₂,…,e_M]，J_M＝[e₁,…,e_M-1]；

步骤4.2、根据步骤2得到的协方差矩阵和步骤4.1得到块对角矩阵J，计算噪声白化矩阵W：

步骤4.3、建立去噪后嵌套阵波达方向估计稀疏模型；

利用步骤4.1得到的块对角矩阵J，剔除步骤3稀疏模型中的噪声项，得到去噪后嵌套阵波达方向估计稀疏模型：

y_J＝Jy_v＝JΨ(Θ)p_Θ+JΔy_v；

步骤4.4、根据步骤4.2得到的噪声白化矩阵W及步骤4.3得到的稀疏模型，建立噪声白化后的波达方向估计稀疏模型y_w：

y_w＝W^-1/2y_J＝W^-1/2JΨ(Θ)p_Θ+ε＝Φ(Θ)p_Θ+ε，

其中，Φ(Θ)＝W^-1/2JΨ(Θ)，服从复高斯白噪声分布。

6.根据权利要求5所述的基于最大似然的嵌套阵稀疏表示波达方向估计方法，其特征在于，所述步骤5具体如下：

步骤5.1、假设稀疏向量p_Θ服从复高斯分布，即

其中，P(·|·)表示条件概率，Γ＝diag(γ₁，γ₂，…，γ_N)，diag(·)表示对角操作运算；

步骤5.2、根据步骤4.4中Φ(Θ)、y_w以及步骤5.1中Γ，分别计算稀疏向量p_Θ的均值和协方差矩阵/>

步骤5.3、根据步骤5.2得到的稀疏向量p_Θ的均值和协方差矩阵/>计算栅格最大似然估计/>

其中，ζ是一个极小正数，(·)^q表示第q次迭代，(·)_n,n表示矩阵第n行第n列的元素；

步骤5.4、根据步骤5.3得到的所有的栅格最大似然估计值形成空间谱，由谱峰位置计算协方差矩阵Σ_-k：

Σ_-k＝Φ(Θ_-k)diag(γ_-k)Φ(Θ_-k)^H+I_{M(M-1)×M(M-1)}，

其中，表示从集合Θ中删除空间谱谱峰所对应的第k个信号源对应的栅格/> 表示从步骤5.3得到的栅格最大似然估计结果中删除第k个信号源对应的/>的估计值；

步骤5.5、计算参数

步骤5.6、根据步骤5.4得到的协方差矩阵Σ_-k和步骤5.5得到的参数计算目标波达方向最大似然估计θ_k：

其中，Re{·}表示取实部运算，/>表示第k个信号源对应的栅格/>左右领域的角度集合，argmax[·]表示取函数最大值的变量值。