CN109116337B - 一种基于矩阵滤波的稀疏近似最小方差方位估计方法 - Google Patents

Abstract

本发明涉及一种基于矩阵滤波的稀疏近似最小方差方位估计方法，使用矩阵滤波器作为预处理器处理阵列接收信号，可抑制强干扰信号，避免其掩蔽弱目标或者影响后续DOA估计算法对弱目标的定位精度；使用SAMV算法进行DOA估计，可在小快拍和低信噪比的情况下保持高分辨性能，并且可以处理相干信号的方位估计问题。同时，整个算法仅需提供迭代停止的门限值η，避免了正则参数的选择，增强了算法的实用性。除此之外，当干扰的阵列流形变化破坏了协方差矩阵的结构时，SAMV算法仍可实现弱目标信号的DOA估计。

Description

一种基于矩阵滤波的稀疏近似最小方差方位估计方法

技术领域

本发明信号处理等领域，涉及一种基于矩阵滤波的稀疏近似最小方差方位估计方法，适用于强干扰环境下弱目标信号的方位估计，涉及信号处理等领域。

背景技术

被动声纳是进行水下目标方位(Direction of Arrival,DOA)估计的有效工具。不同于主动声纳自主发射信号并通过接收反射回波进行目标探测，被动声纳是通过接收舰船辐射噪声来进行目标探测的，因而具有更好的隐蔽性。但当目标信号周围存在强干扰声源如拖船噪声时，强干扰声源将影响DOA估计算法对弱目标信号的定位精度甚至掩蔽弱目标信号。

矩阵滤波器(Matrix Filter,MF)作为一种强干扰抑制的有效手段，被广泛应用于被动声纳系统中。它将空间划分为通带区域和阻带区域，使通带内的信号尽可能无失真地通过，同时抑制阻带内的信号。在现存的基于矩阵滤波器的DOA估计算法中，一般采用多重信号分类(Multiple Signal Classification,MUSIC)算法完成后续DOA估计工作。MUSIC算法通过将协方差矩阵划分为信号子空间和噪声子空间，利用信号子空间与噪声子空间的正交性实现DOA估计。但考虑到在强干扰环境下弱目标信号容易被强干扰所掩蔽，弱目标信号个数的先验信息一般难以获取。在缺少该先验信息的情况下，信号子空间和噪声子空间的划分难以实现，不正确的子空间划分严重影响了MUSIC算法的性能。此外，MUSIC算法也无法处理相干信号的DOA估计问题，在实际使用时受到了较大的限制。

国内学者杨龙和杨益新等(Long Yang,Yixin Yang,and Jiannan Zhu.SourceLocalization based on Sparse Spectral Fitting and Spatial Filtering[C].Oceans2016MTS/IEEE Monterey,2016:1-4)将MF与稀疏谱拟合算法(Sparse Spectrum Fitting,SpSF)相结合，提出SpSF-MF算法。SpSF算法是一种稀疏信号处理类的DOA估计算法，它利用阵列输出数据的二阶统计量信息，通过扫描网格信号功率的l₁范数和协方差矩阵拟合误差的l₂范数联合最小化实现DOA估计。与MUSIC算法相比，该算法可适用于低信噪比和小快拍情况下，同时还可以处理相干信号的DOA估计问题。因此，SpSF-MF算法的性能要远远优于基于矩阵滤波器的MUSIC算法。但SpSF-MF算法在求解时需要提供一个正则参数，由于该参数的选取通常较为困难，使得SpSF-MF在实际信号处理中难以应用。除此之外，干扰通过矩阵滤波器后阵列流形会发生改变，当矩阵滤波器输出的干扰功率仍强于弱目标信号时，输出信号的协方差矩阵的结构会由于干扰阵列流形的改变而被破坏，这将导致SpSF-MF算法的失效。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于矩阵滤波的稀疏近似最小方差方位估计方法，在小快拍、低信噪比的情况下实现强干扰环境下弱目标信号的DOA估计，同时避免正则参数的选取以及减小干扰通过MF后阵列流形变化对DOA估计的影响。本发明给出一种基于MF的稀疏近似最小方差(Sparse Asymptotic Minimum Variance,SAMV)算法(SAMV-MF)，在无需提供正则参数的条件下实现弱目标信号的DOA估计，并且减小了干扰阵列流形变化对DOA估计的影响。

技术方案

一种基于矩阵滤波的稀疏近似最小方差方位估计方法，其特征在于步骤如下：

步骤1：采用阵元间距为d的M元均匀线列阵作为接收阵列，均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(t),0≤t≤N_T,i＝1,...,M；

将接收的时域信号通过希尔伯特变换转化为解析信号，并均匀划分为N段，在每段中进行傅里叶变换将信号划分至子带中，第l个子带内的阵列接收信号可表示为

x_l(n)＝A(f_l,Θ)s_l(n)+e_l(n),n＝1,...,N

其中：

和

分别为各阵元接收信号、信号源以及各阵元接收的环境噪声的傅里叶变换系数组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，对于均匀线列阵来说，

f_l为第l个子带的中心频率，c为声速；Θ＝[θ₁,θ₂,...,θ_Q]是空间[-90°,90°]内Q个离散网格点所代表的方向角组成的向量，所述90°为端射方向；

步骤2：将在第l个子带上设计的通带范围为[θ_pL,θ_pR]、阻带范围为[-90°,θ_sL]∪[θ_sR,90°]的矩阵滤波器转化为在中心频率为f'的子带上设计、的通带范围为：

[arcsin(f_l sin(θ_pL)/f'),arcsin(f_l sin(θ_pR)/f')]

阻带范围为[-90°,arcsin(f_lsin(θ_sL)/f')]∪[arcsin(f_lsin(θ_sR)/f'),90°]的矩阵滤波器，其中θ_pL和θ_pR以及θ_sL和θ_sR分别为在第l个子带上所设计的矩阵滤波器的通带边界和阻带边界，θ_sL＜θ_pL，θ_sR＞θ_pR；

第l个子带上的矩阵滤波器：

min||G^H(f_l)A(f',Θ'_P)-A(f',Θ'_P)||_F

s.t.||G^H(f_l)a(f',θ′_s)||₂≤ε,θ′_s∈Θ'_S

其中：

表示在第l个子带上的矩阵滤波器；ε表示阻带衰减量；||·||₂和||·||_F分别表示l₂范数和Frobenius范数；A(f',Θ'_P)＝[a(f',θ'_pL),...,a(f',θ'_pR)]，θ'_pL和θ'_pR分别表示Θ'_P的左边界和右边界；Θ'_P和Θ'_S表示第l个子带的矩阵滤波器在中心频率为f'的子带上所对应的通带范围和阻带范围；

步骤3、计算第l个子带内矩阵滤波输出信号的采样协方差矩阵：

其中：y_l(n)＝G^H(f_l)x_l(n)为矩阵滤波的输出信号；A(f_l,Θ)被简记为A；E[·]表示期望算子；

为信号协方差矩阵，diag(·)表示对角矩阵，p_q,q＝1,...,Q表示信号功率；

E[e_l(n)e_l ^H(n)]＝σ²I为噪声协方差矩阵，σ²代表噪声功率；令C＝G^H(f_l)A(f_l,Θ)并带入上式中，得

c_q表示矩阵C的第q列，b_i为矩阵G^H(f_l)的第i列；

步骤4：根据SAMV算法，采用迭代的方式计算出第l个子带上的信号功率：

其中：

为矩阵滤波器输出信号的采样协方差矩阵；

和

分别为第i次迭代时第q个信号的功率和噪声功率；

采用迭代的方式计算第l个子带上的噪声功率：

将所求的M个噪声功率进行算数平均，作为该次迭代所计算的噪声功率值

迭代的初始值为：

当前后两次迭代得到的信号与噪声功率组成的向量差值的l₂范数小于所选的迭代停止门限值η时，迭代终止；

步骤5：将所有子带内的信号功率进行叠加，即可得到宽带信号的功率谱。在通带范围内，功率谱的峰值所对应的角度即为弱目标信号的方位角。

有益效果

本发明提出的一种基于矩阵滤波的稀疏近似最小方差方位估计方法，使用矩阵滤波器作为预处理器处理阵列接收信号，可抑制强干扰信号，避免其掩蔽弱目标或者影响后续DOA估计算法对弱目标的定位精度；使用SAMV算法进行DOA估计，可在小快拍和低信噪比的情况下保持高分辨性能，并且可以处理相干信号的方位估计问题。同时，整个算法仅需提供迭代停止的门限值η，避免了正则参数的选择，增强了算法的实用性。除此之外，当干扰的阵列流形变化破坏了协方差矩阵的结构时，SAMV算法仍可实现弱目标信号的DOA估计。

附图说明

图1：采用SAMV-MF方法在强干扰环境下对弱目标信号的方位估计总体流程

图2：各子带内矩阵滤波器的幅度响应

图3：SAMV算法方位估计结果

图4：SAMV-MF方位估计结果

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明解决其技术问题所采用的技术方案包括以下步骤：

1)建立接收信号模型

采用阵元间距为d的M元均匀线列阵作为接收阵列，接收水下目标的辐射噪声。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),0≤n≤N,i＝1,...,M。将接收的时域信号通过希尔伯特变换转化为解析信号，并均匀划分为N段。考虑到辐射噪声一般为宽带信号，通常将其转化为窄带信号进行处理。因此，需对每段信号进行傅里叶变换，将宽带信号划分至各个子带，各子带上的信号可视为窄带信号。

将空间[-90°,90°](其中90°为端射方向)划分为Q个网格，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]。在第l个子带中，阵列的接收信号可表示为x_l(n)＝A(f_l,Θ)s_l(n)+e_l(n),n＝1,...,N，其中

和

分别为各阵元接收信号、信号源以及各阵元接收的环境噪声的傅里叶变换系数组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，对于均匀线列阵来说，

f_l为第l个子带的中心频率，c为声速。

2)设计矩阵滤波器

为保证各个子带内所设计的矩阵滤波器的特性基本相同，需将矩阵滤波器设计问题转化至同一子带中进行。对于第l个子带，其阵列流形中第m个元素可表示为

可以看出，从θ方向上来的中心频率为f_l的窄带信号所对应的阵列流形等价于从arcsin(f_l sin(θ)/f')方向上来的中心频率为f'的窄带信号的阵列流形。因此在频率f_l上设计通带范围为[θ_pL,θ_pR]、阻带范围为[-90°,θ_sL]∪[θ_sR,90°]的矩阵滤波器可转化为在频率f'上设计通带范围为[arcsin(f_l sin(θ_pL)/f'),arcsin(f_l sin(θ_pR)/f')]、阻带范围为[-90°,arcsin(f_l sin(θ_sL)/f')]∪[arcsin(f_l sin(θ_sR)/f'),90°]的矩阵滤波器，其中θ_pL和θ_pR以及θ_sL和θ_sR分别为在中心频率为f_l的子带上设计的矩阵滤波器的通带边界和阻带边界，θ_sL＜θ_pL，θ_sR＞θ_pR。

使用Θ_P和Θ_S表示第l个子带上矩阵滤波器的通带范围和阻带范围，并使用Θ'_P和Θ'_S表示将第l个子带转化至中心频率为f'的子带上矩阵滤波器的通带范围和阻带范围。根据矩阵滤波器使通带内信号尽可能无失真地通过并抑制阻带内信号的性质，将第l个子带上所对应的矩阵滤波器设计问题表示为：

其中

表示所设计的矩阵滤波器；ε表示阻带衰减量；||·||₂和||·||_F分别表示l₂范数和Frobenius范数；A(f',Θ'_P)＝[a(f',θ'_pL),...,a(f',θ'_pR)]，θ'_pL和θ'_pR分别表示Θ'_P的左边界和右边界。

考虑到|f_l sin(θ)/f'|≤1，所选择的矩阵滤波器的设计频率f'必须满足f'≥f_lmax[|sin(θ_pL)|,|sin(θ_pR)|,|sin(θ_sL)|,|sin(θ_sR)|]。当f'大于基阵的工作频率f_c时，可能会出现a(f',θ₁)＝a(f',θ₂),θ₁≠θ₂的情况，若θ₁∈Θ'_P，θ₂∈Θ'_S，则矩阵滤波器的设计将出现错误，因此所选择的频率f'还需满足f'≤f_c。

3)基于矩阵滤波器的稀疏近似最小方差算法

第l个子带上矩阵滤波器的输出信号可表示为y_l(n)＝G^H(f_l)x_l(n),n＝1,...,N，上标“H”表示共轭转置。假设环境噪声为高斯白噪声，矩阵滤波器输出信号的协方差矩阵可计算为

其中A(f_l,Θ)被简记为A；E[·]表示期望算子；

为信号协方差矩阵，diag(·)表示对角矩阵，p_q,q＝1,...,Q表示信号功率；

为噪声协方差矩阵，σ²代表噪声功率。令

并带入到式(3)中，得

c_q表示矩阵C的第q列，b_i为矩阵G^H(f_l)的第i列。可以看出，噪声可看作M个功率为σ²的信号入射至阵列上，它们所对应的阵列流形为b_i,i＝1,...,M。定义一个新的矩阵W＝[c₁,c₂,...,c_Q,b₁,...,b_M]，输出信号的协方差矩阵可重新表示为R_MF＝WP'W^H，其中P'＝diag(p₁,...,p_Q,σ²,...,σ²)。

根据SAMV算法，采用迭代的方式计算出信号功率。信号功率的迭代式如下：

其中，

为矩阵滤波器输出信号的采样协方差矩阵；

第q个信号在第i次迭代时的功率。对于噪声功率，对应的迭代式如下：

在式(4)和(5)中，

从

的表达式中可以看出，估计矩阵滤波器输出信号的协方差矩阵时仅需要一个噪声功率值，因此需对所得的M个噪声功率进行算数平均，作为该次迭代所计算的噪声功率值。

迭代的初始值可由

提供。当前后两次迭代得到的信号与噪声功率组成的向量差值的l₂范数小于所选的迭代停止门限值η时，迭代终止。

当完成各子带内的信号功率估计后，需将所有子带内的信号功率进行叠加，即可得到宽带信号的功率谱。在通带范围内，功率谱的峰值所对应的角度即为弱目标信号的方位角。

具体实施例子：

三个宽带信号入射到阵元间距为4m的32元均匀线列阵上,其中一个为强干扰信号，另外两个为弱目标信号，三个信号的频率范围均为[90,180]Hz。强干扰信号的入射角为20°，干噪比为20dB；两个弱目标信噪比均为0dB，分别从-5°和-8°入射至阵列上。信号采样率为5,120Hz，将采集的信号按重叠率50％划分为40组，每组中进行1024点的傅里叶变换。方位扫描网格为-90°～90°(其中90°方向为端射方向)，扫描角度为1°。

选择180Hz为设计频率，在各个子带内设计通带范围为[-40°,0°]，阻带范围为[-90°,-45°]∪[5°,90°]的矩阵滤波器，阻带衰减级为-15dB。计算各个子带上矩阵滤波器幅度响应

结果见图2。从图2中可以看出，在各个子带上设计的矩阵滤波器具有相同的通带区域和阻带区域。通带范围内的幅度响应均为0dB，可使信号无失真地通过；阻带衰减级均控制在为-15dB以下，以抑制阻带内的信号，符合设计要求。

图1为采用SAMV-MF方法在强干扰环境下对弱目标信号进行方位估计的总体流程，具体实施如下：

1)采用阵元间距为d的M元均匀线列阵作为接收阵列，接收水下目标辐射噪声。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),0≤n≤N,i＝1,...,M。将接收信号通过希尔伯特变换转化为解析信号，并划分为N段，在每段中进行傅里叶变换。观察任意一段信号的频谱图，选择合适的频段进行阵列信号处理。将空间[-90°,90°](其中90°为端射方向)划分为Q个网格点，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]。在第l个子带中，阵列的接收信号可表示为x_l(n)＝A(f_l,Θ)s_l(n)+e_l(n),n＝1,...,N，其中

和

分别为各阵元接收信号、信号源以及各阵元接收的环境噪声的傅里叶变换系数组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，对于均匀线列阵来说，

f_l为第l个子带的中心频率，c为声速。

2)采用CBF算法对阵列接收信号进行方位估计，得到强干扰信号的分布信息。根据该分布信息，将空间划分为通带区域Θ_P和阻带区域Θ_S，其中通带区域为感兴趣的区域，而阻带区域需覆盖所有干扰所在方位。为在每个子带中得到特性相同的矩阵滤波器，需将在中心频率为f_l的子带上设计通带范围为[θ_pL,θ_pR]、阻带范围为[-90°,θ_sL]∪[θ_pR,90°]的矩阵滤波器问题转化至在中心频率为f'的子带上设计通带范围为[arcsin(f_l sin(θ_pL)/f'),arcsin(f_l sin(θ_pR)/f')]、阻带范围为[-90°,arcsin(f_l sin(θ_sL)/f')]∪[arcsin(f_l sin(θ_sR)/f'),90°]的矩阵滤波器问题，其中θ_pL和θ_pR以及θ_sL和θ_sR分别为在中心频率为f_l的子带上设计的矩阵滤波器的通带边界和阻带边界，θ_sL＜θ_pL，θ_sR＞θ_pR。所选择的设计频率f'必须满足条件f_l max[|sin(θ_pL)|,|sin(θ_pR)|,|sin(θ_sL)|,|sin(θ_sR)|]≤f'≤f_c，f_c为阵列的工作频率。第l个子带上所对应的矩阵滤波器设计问题表示为：

其中

表示矩阵滤波器；ε表示阻带衰减量；||·||₂和||·||_F分别表示l₂范数和Frobenius范数；Θ'_P和Θ'_S表示将第l个子带转化至中心频率为f'的子带上矩阵滤波器的通带范围和阻带范围；A(f',Θ'_P)＝[a(f',θ'_pL),...,a(f',θ'_pR)]，θ'_pL和θ'_pR分别表示Θ'_P的左边界和右边界；上标“H”表示矩阵共轭转置。

3)在各子带内使用设计好的矩阵滤波器对阵列接收信号进行滤波，抑制强干扰源，第l个子带上矩阵滤波器的输出可表示为y_l(n)＝G^H(f_l)x_l(n),n＝1,...,N。考虑环境噪声为高斯白噪声的情况，对于第l个子带，矩阵滤波器输出信号的协方差矩阵可表示为

其中E[·]表示期望算子；C＝G^H(f_l)A(f_l,Θ)，c_q表示矩阵C的第q列；P＝diag(p₁,...,p_Q)为信号协方差矩阵，diag(·)表示对角矩阵，p_q,q＝1,...,Q表示信号功率；σ²代表噪声功率；b_i为矩阵G^H(f_l)的第i列。可以看出，噪声可看作M个功率为σ²的信号入射至阵列上，它们所对应的阵列流形为b_i,i＝1,...,M。利用SAMV算法采用迭代的方式进行弱目标信号方位估计，迭代式如下：

其中

和

分别为第i次迭代时计算的第q个信号的功率和第j个噪声功率值；

为第i次迭代得到的M个噪声功率的平均值；

为矩阵滤波器输出信号的采样协方差矩阵；

迭代的初始值可由

提供。选择合适的迭代停止门限η，当前后两次迭代得到的信号与噪声功率组成的向量差值的l₂范数小于迭代停止门限值η时，迭代终止。将各个子带上的方位谱进行叠加，得到宽带信号的方位谱，在通带范围内方位谱峰值对应的角度即为感兴趣的弱目标信号的方位角。

分别使用SAMV和SAMV-MF进行弱目标方位估计，结果分别见图3和图4，其中图上蓝色虚线为真实信号方位，红色圆点为估计的信号方位。从图中可以看出，由于受到强干扰的影响，弱目标被强干扰所掩蔽，SAMV算法无法分辨出弱目标，而本发明则可同时估计出弱目标信号和强干扰所在方位。通过仿真数据可以看出，本发明能够实现强干扰环境下弱目标信号的方位估计。

Claims (1)

1.一种基于矩阵滤波的稀疏近似最小方差方位估计方法，其特征在于步骤如下：

步骤1：采用阵元间距为d的M元均匀线列阵作为接收阵列，均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(t),0≤t≤N_T,i＝1,...,M；

将接收的时域信号通过希尔伯特变换转化为解析信号，并均匀划分为N段，在每段中进行傅里叶变换将信号划分至子带中，第l个子带内的阵列接收信号可表示为

x_l(n)＝A(f_l,Θ)s_l(n)+e_l(n),n＝1,...,N

其中：

和

分别为各阵元接收信号、信号源以及各阵元接收的环境噪声的傅里叶变换系数组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，对于均匀线列阵来说，

f_l为第l个子带的中心频率，c为声速；Θ＝[θ₁,θ₂,...,θ_Q]是空间[-90°,90°]内Q个离散网格点所代表的方向角组成的向量，所述90°为端射方向；

步骤2：将在第l个子带上设计的通带范围为[θ_pL,θ_pR]、阻带范围为[-90°,θ_sL]∪[θ_sR,90°]的矩阵滤波器转化为在中心频率为f'的子带上设计的通带范围为：

[arcsin(f_lsin(θ_pL)/f'),arcsin(f_lsin(θ_pR)/f')]

阻带范围为[-90°,arcsin(f_lsin(θ_sL)/f')]∪[arcsin(f_lsin(θ_sR)/f'),90°]的矩阵滤波器，其中θ_pL和θ_pR以及θ_sL和θ_sR分别为在第l个子带上所设计的矩阵滤波器的通带边界和阻带边界，θ_sL<θ_pL，θ_sR>θ_pR；

第l个子带上的矩阵滤波器：

min||G^H(f_l)A(f',Θ'_P)-A(f',Θ'_P)||_F

s.t.||G^H(f_l)a(f',θ′_s)||₂≤ε,θ_s'∈Θ'_S

其中：

表示在第l个子带上的矩阵滤波器；ε表示阻带衰减量；||·||₂和||·||_F分别表示l₂范数和Frobenius范数；A(f',Θ'_P)＝[a(f',θ'_pL),...,a(f',θ'_pR)]，θ'_pL和θ'_pR分别表示Θ'_P的左边界和右边界；Θ'_P和Θ'_S表示第l个子带的矩阵滤波器在中心频率为f'的子带上所对应的通带范围和阻带范围；

步骤3、计算第l个子带内矩阵滤波输出信号的采样协方差矩阵：

其中：y_l(n)＝G^H(f_l)x_l(n)为矩阵滤波的输出信号；A(f_l,Θ)被简记为A；E[·]表示期望算子；

为信号协方差矩阵，diag(·)表示对角矩阵，p_q,q＝1,...,Q表示信号功率；

为噪声协方差矩阵，σ²代表噪声功率；令C＝G^H(f_l)A(f_l,Θ)并带入上式中，得

c_q表示矩阵C的第q列，b_i为矩阵G^H(f_l)的第i列；

步骤4：根据SAMV算法，采用迭代的方式计算出第l个子带上的信号功率：

其中：

为矩阵滤波器输出信号的采样协方差矩阵；

和

分别为第k次迭代时第q个信号的功率和噪声功率；

采用迭代的方式计算第l个子带上的噪声功率：

将所求的M个噪声功率进行算数平均，作为该次迭代所计算的噪声功率值

迭代的初始值为：

当前后两次迭代得到的信号与噪声功率组成的向量差值的l₂范数小于所选的迭代停止门限值η时，迭代终止；

步骤5：将所有子带内的信号功率进行叠加，即可得到宽带信号的功率谱；在通带范围内，功率谱的峰值所对应的角度即为弱目标信号的方位角。

Images