CN113109760A - 一种基于组稀疏的多线谱联合doa估计和聚类方法及系统 - Google Patents

Abstract

本发明公开了一种基于组稀疏的多线谱联合DOA估计和聚类方法，包括：步骤1、获取均匀线阵水听器的观测阵列信号；步骤2、基于常规波束形成算法粗略估计目标方位

步骤3、对

方向做跟踪波束，检测目标信号的L个信噪比最大的线谱；步骤4、计算L根线谱对应频域值的协方差矩阵R(f_l)；步骤5、将其转化为压缩感知问题，通过在目标函数中添加l₁范数和l₂范数，联合求解L根线谱的空间谱P(θ,f_l)的稀疏解；步骤6、对L个线谱的空间谱进行聚类,将L根线谱分为来自目标方位

以及来自非目标方位

两类。本发明的一种基于组稀疏的多线谱联合DOA估计和聚类方法实现了在多目标强干扰环境下多线谱的提取、DOA估计以及分类，应用简单直接，经济代价低且效果明显，运算量较小，分类效果好。

Description

一种基于组稀疏的多线谱联合DOA估计和聚类方法及系统

技术领域

本发明属于声呐信号处理技术领域，具体涉及一种在多目标强干扰环境中基于组稀疏的多线谱联合估计DOA方法及系统。

背景技术

水声目标通常处于复杂的海洋环境中，多目标强干扰、环境噪声干扰等因素给基于目标特征的被动声呐识目标识别带来了很大的困难。水下目标辐射噪声特征分析是实现目标探测、识别和分类的基础，线谱是水下目标检测与识别的相当重要的特征。

水下目标的辐射噪声的功率谱通常由连续谱和线谱构成。线谱通常具有较好的相位稳定性，比连续谱拥有更高的信噪比，是被动声呐目标检测的重要方式，也是全世界各国水声信号处理科研人员的研究热点。

现有的线谱检测方法主要分为以下几类：1.离散傅里叶基线谱检测方法，例如参数模型类谱估计方法，包括自回归模型(Auto regressive model，AR)、滑动平均模型(moving average model,MA)以及自回归滑动平均模型(Auto regressive movingaverage model，ARMA)模型法。其功率谱估计结果较平滑，频率分辨率较高，但是其模型的阶次不易选择。子空间分解类谱估计典型方法有多信号分类方法(Multiple SignalClassification，MUSIC)、基于旋转不变技术的信号参数估计方法(Estimating SignalParameter via Rotational Invariance Techniques，ESPRIT)。子空间分解类方法可以获得较高的频率分辨率，但是需要信号子空间维数的先验信息，且在低信噪比、色噪声等背景下性能退化严重。2.自适应线谱增强器方法，利用线谱信号的时间相关半径大于噪声信号的时间相关半径的特征，对接收信号中的宽带信号进行解相干延时，保持周期信号相干，进而对二者进行LMS自适应抵消处理，分离宽带信号，增强周期信号。3.高阶谱方法，高阶谱能够更全面地描述舰船的辐射噪声，目前被更多地用于舰船的目标识别。但目前高阶谱的物理意义不明确，有碍直观地对其进行分析和理解，高阶谱用于线谱检测还需要更深入的研究。4.稀疏重构类方法，正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法为一种经典的稀疏重构类方法。但是该方法需要事先知道目标线谱个数，在实际条件下很难满足要求。基于迭代最小化的稀疏学习(Sparse Learning via Iterative Minimization,SLIM)方法采用迭代的方式对线谱信号频率、幅度及背景噪声功率进行估计，无需目标线谱个数信息，参数设置相对简单。5.非线性混沌振子类方法，这是一种非线性处理方法，仿照非线性动力学系统，构建一个处于混沌状态的系统模型。该方法对噪声信号不敏感，实现了对微弱线谱信号的有效监测，但在系统模型的状态定量判断、算法的复杂度等方面还有待于改进。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供了一种基于组稀疏的多线谱联合DOA估计和聚类方法，该方法能够同时估计检测得到的多根线谱的来波方向，并将多根线谱分为来自目标方向和非目标方向两类，实现了在多目标强干扰水声环境下对线谱的提取、DOA估计以及分类。

技术方案：一种基于组稀疏的多线谱联合DOA估计和聚类方法，包括：

步骤1，获取均匀线阵(ULA)水听器的阵元接收到的时域信号x_i(t)，i＝1,2,...,M，i为均匀线阵中阵元序号，M为均匀线阵中阵元数目；

步骤2，基于常规波束形成算法(CBF)粗略估计目标方位θ，θ为目标信号波束能量最大时的引导角；步骤3，对步骤(2)中估计得到的目标方位θ做跟踪波束，得到波束形成后的阵列信号，并对阵列信号作FFT得到阵列信号的频谱图，在频谱图中搜索检测出L根信噪比最大的线谱；

步骤4，计算L根线谱对应的频域值的协方差矩阵R(f_l)，计算L根线谱对应的导向矢量矩阵D(f_l)，计算对应的噪声功率σ_n ²；

步骤5，构造目标函数

通过添加l₁范数和l₂范数约束，将其转化为压缩感知的问题，联合求解L根线谱的空间谱P(θ,f_l)的稀疏解，

为协方差矩阵的估计值，I为单位矩阵；

步骤6，使用经典层次聚类算法对联合估计得到的L根线谱的空间谱进行聚类，将其分为来自目标方位

以及非目标方位

两类。

进一步地，所述步骤2包括：

步骤2.1，计算理想均匀线阵(ULA)在引导角θ_j下相邻阵元时延τ_j：

其中j＝1,2,...,J+1，j为引导角序号，J+1为总引导角个数，引导角范围为：0～180°，d为相邻阵元间的距离，v是声音在水中的传播速度；

步骤2.2，对各阵元接收的阵列信号进行延时相加，得到目标信号波束能量图B[b(1),…,b(J+1)]，其中b(j)为：

其中，M表示均匀线阵水听器的观测阵列个数，τ_j表示阵元的时延；

步骤2.3，通过对波束能量图的检测搜索找到波束能量最大值位置时的引导角为目标方位的粗估计θ。

进一步地，所述步骤3包括：

步骤3.1，根据粗估计的目标方位θ计算各阵元的时延估计

其中v是声音在水中的传播速度，d为相邻阵元间的距离；

步骤3.2，将各阵元数据按时延估计

与参考阵元对齐，将对齐后的阵元数据相干相加获得目标跟踪波束g(t)：

步骤3.3，对g(t)进行傅里叶变换获得目标信号频谱G(f)，同时利用滑动窗平滑技术估计目标信号连续谱G_c(f)，在目标信号频谱G(f)中删除连续谱G_c(f)的影响，获得目标信号的线谱G_line(f)：

G_line(f)＝G(f)-G_c(f)

获得目标的线谱G_line(f)之后，通过计算对应频率f处强线谱的信噪比，通过排序算法得到L根信噪比最大的强线谱，强线谱对应频率为f_l,l＝1,2,3,...,L。

进一步地，所述步骤4包括：

步骤4.1，对阵列接收的信号进行分帧处理，每一帧信号的采样点数为S，信号重叠部分长度为S/2，共分为T帧数据；

步骤4.2，对每一帧信号中单个阵元接收信号x_m(t)作傅里叶变换得到x_fft,m(t)，得到L根线谱频率f_l,l＝1,2,...,L处的值x_fft,m(t,f_k)，包含了相位和幅度信息，可以得到每一帧中的线谱阵列向量：

x_fft(t,f_k)＝[x_fft,1(t,f_k),...x_fft,M(t,f_k)]^T；

步骤4.3，利用T帧数据，对得到的L根信噪比最大的线谱，求每一根线谱f_l其频域值的协方差矩阵为：

步骤4.4，计算L根线谱对应的导向矢量矩阵D(f_l)：

上式中矩阵的导向矢量具体表达式为:

步骤4.5，计算噪声功率σ_n ²，用每一根线谱f_l其频域值的协方差矩阵R(f_l)的最小特征值作为噪声功率的估计值

进一步地，所述步骤5包括：

步骤5.1，在阵列信号处理中，信源的数目被认为少于传感器阵列的数目，因此L根线谱对应的空间谱分布P(θ,f_l)认为是稀疏的；

步骤5.2，实际的水声信号处理过程中，精确的干扰加噪声协方差矩阵R_i+n不易获取，我们通常采用样本的协方差矩阵代替

当采样点数S足够多时，样本的协方差矩阵近似等于干扰加噪声协方差矩阵，即

用L根线谱对应的频域值的协方差矩阵

近似等于精确的干扰加噪声协方差矩阵R_i+n；

步骤5.3，基于上述步骤的说明，构造了多线谱空间谱的联合估计稀疏约束优化问题：

其中，矩阵P(f_l)是频率为f_l强线谱的空间谱分布p(f_l)的对角矩阵，即P(f_l)＝diag{p(f_l)}，

代表频率为f_l的噪声的误差的方差，||·||_F代表矩阵的Frobenius范数，||·||₀代表向量的l₀范数，在这里表示强线谱空间谱中的非0元素个数，系数γ控制着强线谱空间谱的稀疏度和剩余范数之间的权衡；

步骤5.4，步骤5.2中的稀疏约束优化目标函数中由于l₀范数的存在，难以求解，本发明中l₁范数代替l₀范数作为该优化目标函数的近似解，即||p(f_l)||₀＝||p(f_l)||₁；

步骤5.5，对于线谱检测估计得到的L根信噪比最大的线谱，对于来自目标信号方向上的线谱，其空间谱的分布p(f_l)可以认为是近似相等的，即引导角范围0～180°上的值近似相等，引入了组稀疏(Sparse-group)的概念，将0～180°分为181组，第i组数据中包含L根强线谱在引导角θ_i上的功率{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}，对于每一组中的强线谱的功率，本发明中认为该组中的功率要么是非零值，要么是零值，即181组L根强线谱的功率是组稀疏的；

步骤5.6，通过在第i组数据p⁽ⁱ⁾＝{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}中添加l₂范数的形式，实现组稀疏(Sparse-group)；

步骤5.7，对于第i组数据中包含L根强线谱在引导角θ_i上的功率{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}，对于来自不同方向的强线谱不一定能保证该组中的数据都是非零值或是零值，因此对于每一组中的数据是稀疏的，本发明中引入l₁范数代替l₀范数保证每一组中数据的稀疏性，从而保证联合估计L根强线谱的空间谱更加接近真实的空间谱分布；

步骤5.8，引入组稀疏的概念，本发明改进了步骤(5.3)中优化目标函数，重新构造L根强线谱的组稀疏回归约束(Sparse-group lasso)优化目标函数：

上式中p_i为第i组数据的长度，J+1为总引导角个数，引导角范围为：0～180°，即L根强线谱空间谱分组的组数，α是一个权重系数，α∈[0,1]，α用来控制组稀疏和组内稀疏性的权衡，经过验证，上述优化目标函数为凸函数，可以求解L根强线谱的空间谱分布。

进一步地，所述步骤6包括：

步骤6.1，计算L根强线谱空间谱p(f_l)两两之间的Frechet距离，构成矩阵

其中P_ij表示第i个强线谱空间谱分布与第j个强线谱空间谱分布之间的Frechet距离；

步骤6.2，以Frechet距离矩阵P作为距离矩阵，进行层次聚类，对L个强线谱的空间谱分布进行分类得到分类簇

其中c_θ包含来自目标方向的强线谱的空间谱分布，

包含来自非目标方向的强线谱的空间谱分布。

一种基于组稀疏的多线谱联合估计DOA系统，包括依次连接的观测阵列信号获取模块、目标方位粗估计模块、线谱检测模块、强线谱相关参数计算模块、强线谱的空间谱联合估计模块、强线谱的空间谱聚类模块；

观测阵列信号获取模块，用于获取观测阵列信号x_i(t)，i＝1,2,...,M，M为均匀线阵中(ULA)阵元数目；

目标方位粗估计模块，用于粗略估计目标方位

为目标信号波束能量最大时的引导角；

线谱检测模块，利用估计得到的目标方位

对阵列信号做跟踪波束想干相加后得到g(t)，对g(t)做傅里叶变换后，利用线谱检测技术检测出L根信噪比最大的线谱；

强线谱相关参数计算模块，对阵列接收的信号进行分帧处理，对每一帧信号中单个阵元接收信号x_m(t)作傅里叶变换得到x_fft,m(t)，求每一根线谱f_l其频域值的协方差矩阵R(f_l)。计算L根线谱对应的导向矢量矩阵D(f_l)，计算噪声功率σ_n ²；

强线谱的空间谱联合估计模块，引入组稀疏的概念，构造了L根强线谱的组稀疏回归约束(Sparse-group lasso)优化目标函数，求解L根强线谱的空间谱分布；

强线谱的空间谱聚类模块，计算L根强线谱空间谱p(f_l)两两之间的Frechet距离，将L根强线谱分为来自目标方向

的强线谱，以及来自非目标方向

的强线谱两类。

有益效果：与现有技术相比，本发明公开的多线谱联合估计DOA方法具有以下优点：波束形成直接从接收到的阵元数据提取信噪比大的强线谱，构建强线谱的组稀疏回归约束优化目标函数，联合求解多线谱的空间谱分布，实现了在多目标强干扰环境下多线谱的提取、DOA估计以及分类，应用简单直接，经济代价低且效果明显，运算量较小，分类效果好。

附图说明

图1为实施例一种基于组稀疏的多线谱联合DOA估计和聚类方法的流程图；

图2为基于均匀线阵(ULA)的波束能量图；

图3为强线谱的挑选结果图；

图4为实施例一中基于组稀疏的多线谱联合估计DOA方法的多线谱的空间谱估计图；

图5为实施例一基于传统Capon谱的多线谱的空间谱估计图；

图6为实施例一中多线谱的分类结果图，类别一：来自目标方向，类别二：来自非目标方向；

图7为基于组稀疏的多线谱联合估计DOA系统的组成图；

图8为实施例二中不同信噪比下多线谱联合估计方法与传统Capon谱估计方法的对比图。

具体实施方式

下面结合附图对本发明作进一步说明：

实施例一：

本发明公开了一种基于组稀疏的多线谱联合DOA估计和聚类方法，该方法可以同时估计检测得到的多根线谱的来波方向，并将多根线谱分为来自目标方向和非目标方向两类。实现了在多目标强干扰水声环境下对线谱的提取、DOA估计以及分类。为了验证此方法的有效性，本实施例以均匀线阵(ULA)为例，包括64个阵元，即M＝64，阵元间距有微小差别，本实施例中忽略间距差别，认为阵元间距近似相等，为：d＝1.5米。将左侧第一个阵元作为参考阵元，以参考阵元的位置为原点建立坐标系。本实施例中仿真了舰船的辐射噪声频谱包括连续谱和线谱。本实施例中，添加了八根仿真线谱，其中来自目标方向60°的线谱频率为:59Hz、97Hz、163Hz、198Hz、232Hz。来自非目标方向63°的线谱频率为：125Hz、138Hz以及280Hz。

本实施例公开的基于组稀疏的多线谱联合估计DOA方法具体步骤如图1所示，包括：

步骤1，获取均匀线阵(ULA)水听器的阵元接收到的时域信号x_i(t)，i＝1,2,...,M，M为均匀线阵中(ULA)阵元数目；

步骤2，基于常规波束形成算法(CBF)粗略估计目标方位

为目标信号波束能量最大时的引导角；步骤2.1，计算理想均匀线阵(ULA)在引导角θ_j下相邻阵元时延τ_j：

其中j＝1,2,...,J+1，j为引导角序号，J+1为总引导角个数，引导角范围为：0～180°，d为相邻阵元间的距离，v是声音在水中的传播速度；

步骤2.2，对各阵元接收的阵列信号进行延时相加，得到目标信号波束能量图B[b(1),…,b(J+1)]，其中b(j)为：

其中，M表示均匀线阵水听器的观测阵列个数，τ_j表示阵元的时延；

步骤2.3，通过对波束能量图的检测搜索找到波束能量最大值位置时的引导角为目标方位的粗估计θ。基于理想阵形的波束形成如图2所示，通过能量检测找到波束能量最大值的位置，得到目标方位的粗估计θ＝60°。

步骤3、对步骤2中估计得到的目标方位

做跟踪波束，得到波束形成后的阵列信号，并对阵列信号作FFT得到阵列信号的频谱图，在频谱图中搜索检测出L根信噪比最大的线谱；L为估计线谱的数目，本实施例中取值为12；强线谱检测的结果如图3所示。

步骤3.1，根据粗估计的目标方位θ计算各阵元的时延估计

其中v是声音在水中的传播速度，d为相邻阵元间的距离；

步骤3.2，将各阵元数据按时延估计

与参考阵元对齐，将对齐后的阵元数据相干相加获得目标跟踪波束g(t)：

步骤3.3，对g(t)进行傅里叶变换获得目标信号频谱G(f)，同时利用滑动窗平滑技术估计目标信号连续谱G_c(f)，在目标信号频谱G(f)中删除连续谱G_c(f)的影响，获得目标信号的线谱G_line(f)：

G_line(f)＝G(f)-G_c(f)

获得目标的线谱G_line(f)之后，通过计算对应频率f处强线谱的信噪比，通过排序算法得到L根信噪比最大的强线谱，强线谱对应频率为f_l,l＝1,2,3,...,L。

步骤4，计算L根线谱对应的频域值的协方差矩阵R(f_l)，计算L根线谱对应的导向矢量矩阵D(f_l)，计算对应的噪声功率σ_n ²；

步骤5、构造目标函数

通过添加l₁范数和l₂范数约束，将其转化为压缩感知的问题，联合求解L根线谱的空间谱P(θ,f_l)的稀疏解，

为协方差矩阵的估计值，I为单位矩阵；

步骤5.1，在阵列信号处理中，信源的数目被认为少于传感器阵列的数目，因此L根线谱对应的空间谱分布P(θ,f_l)认为是稀疏的；

步骤5.2，实际的水声信号处理过程中，精确的干扰加噪声协方差矩阵R_i+n不易获取，我们通常采用样本的协方差矩阵代替

当采样点数S足够多时，样本的协方差矩阵近似等于干扰加噪声协方差矩阵，即

用L根线谱对应的频域值的协方差矩阵

近似等于精确的干扰加噪声协方差矩阵R_i+n；

步骤5.3，基于上述步骤的说明，构造了多线谱空间谱的联合估计稀疏约束优化问题：

其中，矩阵P(f_l)是频率为f_l强线谱的空间谱分布p(f_l)的对角矩阵，即P(f_l)＝diag{p(f_l)}，

代表频率为f_l的噪声的误差的方差，||·||_F代表矩阵的Frobenius范数，||·||₀代表向量的l₀范数，在这里表示强线谱空间谱中的非0元素个数，系数γ控制着强线谱空间谱的稀疏度和剩余范数之间的权衡；

步骤5.4，步骤5.2中的稀疏约束优化目标函数中由于l₀范数的存在，难以求解，本发明中l₁范数代替l₀范数作为该优化目标函数的近似解，即||p(f_l)||₀＝||p(f_l)||₁；

步骤5.5，对于线谱检测估计得到的L根信噪比最大的线谱，对于来自目标信号方向上的线谱，其空间谱的分布p(f_l)可以认为是近似相等的，即引导角范围0～180°上的值近似相等，因此本发明中，引入了组稀疏(Sparse-group)的概念，将0～180°分为181组，第i组数据中包含L根强线谱在引导角θ_i上的功率{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}，对于每一组中的强线谱的功率，本发明中认为该组中的功率要么是非零值，要么是零值，即181组L根强线谱的功率是组稀疏的；

步骤5.6，本发明中通过在第i组数据p⁽ⁱ⁾＝{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}中添加l₂范数的形式，实现组稀疏(Sparse-group)；

步骤5.7，对于第i组数据中包含L根强线谱在引导角θ_i上的功率{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}，对于来自不同方向的强线谱不一定能保证该组中的数据都是非零值或是零值，因此对于每一组中的数据是稀疏的，本发明中引入l₁范数代替l₀范数保证每一组中数据的稀疏性，从而保证联合估计L根强线谱的空间谱更加接近真实的空间谱分布；

步骤5.8，引入组稀疏的概念，本发明改进了步骤5.3中优化目标函数，重新构造L根强线谱的组稀疏回归约束(Sparse-group lasso)优化目标函数：

上式中p_i为第i组数据的长度，J+1为总引导角个数，引导角范围为：0～180°，即L根强线谱空间谱分组的组数，α是一个权重系数，α∈[0,1]，α用来控制组稀疏和组内稀疏性的权衡，经过验证，上述优化目标函数为凸函数，可以求解L根强线谱的空间谱分布。

图4给出了使用基于组稀疏的强线谱联合估计DOA方法，求解仿真的8根强线谱的空间谱分布。从图中可以看出本发明公开的方法可以有效的同时估计出多根强线谱的空间谱分布，并且精度较高，对于目标方向60°和非目标方向63°这种入射角度间距较小的两个目标，本发明也能很好的将来自目标方向和非目标方向的线谱区分开来。

图5给出了利用常规DOA估计方法——Capon谱估计方法。本发明中公开的方法对目标的方位分辨能力明显优于常规的Capon谱估计方法。说明该方法对目标的方位分辨能力非常好，主要表现为在强线谱的空间谱估计中有尖锐的峰，并且没有其余杂散的峰。

步骤6、使用经典层次聚类算法对联合估计得到的L根线谱的空间谱进行聚类，将其分为来自目标方位θ以及非目标方位θ_i两类；

步骤6.1，计算L根强线谱空间谱p(f_l)两两之间的Frechet距离，构成矩阵

其中P_ij表示第i个强线谱空间谱分布与第j个强线谱空间谱分布之间的Frechet距离；

步骤6.2，以Frechet距离矩阵P作为距离矩阵，进行层次聚类，对L个强线谱的空间谱分布进行分类得到分类簇

其中c_θ包含来自目标方向的强线谱的空间谱分布，

包含来自非目标方向的强线谱的空间谱分布。

本发明将压缩感知的思想用于水声信号多线谱的联合估计DOA，通过聚类算法，可以将来自目标方位的线谱以及来自于非目标方位的线谱有效分类，选取我们需要的线谱。该发明不同于传统的空间谱估计方法，需要将多根不同频率的线谱分别估计DOA。图6给出了多根强线谱的分类结果。本发明中的方法可以有效地将来自目标方向60°的强线谱和非目标方向63°的强线谱有效地区分开来。

本实施例还公开了实现上述基于组稀疏的多线谱联合估计DOA系统，如图7所示，包括依次连接的观测阵列信号获取模块、目标方位粗估计模块、线谱检测模块、强线谱相关参数计算模块、强线谱的空间谱联合估计模块、强线谱的空间谱聚类模块；

观测阵列信号获取模块，用于获取观测阵列信号x_i(t)，i＝1,2,...,M，M为均匀线阵中(ULA)阵元数目；

目标方位粗估计模块，用于粗略估计目标方位

为目标信号波束能量最大时的引导角；

线谱检测模块，利用估计得到的目标方位

对阵列信号做跟踪波束想干相加后得到g(t)，对g(t)做傅里叶变换后，利用线谱检测技术检测出L根信噪比最大的线谱。

强线谱相关参数计算模块，对阵列接收的信号进行分帧处理，对每一帧信号中单个阵元接收信号x_m(t)作傅里叶变换得到x_fft,m(t)，求每一根线谱f_l其频域值的协方差矩阵R(f_l)。计算L根线谱对应的导向矢量矩阵D(f_l)，计算噪声功率σ_n ²。

强线谱的空间谱联合估计模块，引入组稀疏的概念，构造了L根强线谱的组稀疏回归约束(Sparse-group lasso)优化目标函数，求解L根强线谱的空间谱分布。

强线谱的空间谱聚类模块，计算L根强线谱空间谱p(f_l)两两之间的Frechet距离，将L根强线谱分为来自目标方向θ的强线谱，以及来自非目标方向θ_i的强线谱两类。

实施例二：

为验证本发明公开的方法在不同信噪比情况下的效果，本实施例在实施例一的基础上，改变舰船的辐射噪声频谱包括连续谱和线谱中的信噪比，使信噪比为s，-30dB＜s＜0dB，观察在不同信噪比下本发明方法的强线谱DOA估计结果。实验结果如图8所示。

随着信噪比的提升，本发明的强线谱联合估计DOA方法的估计误差随着信噪比的提高，其强线谱联合估计的空间谱的误差与真实的空间谱误差明显减少。并且在低信噪比情况下，也能达到较好的DOA估计精度，误差较小。在低信噪比情况下，随着信噪比变化，该方法的估计误差变化不明显，说明本发明的方法的鲁棒性较好，在信噪比较低情况下，也能完成对多强线谱的空间谱估计任务。与传统的Capon谱估计算法比较，在低信噪比情况下，本发明的方法的估计误差要明显小于传统的Capon谱估计算法。

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

Claims (7)

1.一种基于组稀疏的多线谱联合DOA估计和聚类方法，其特征在于，包括：

步骤1，获取均匀线阵(ULA)水听器的阵元接收到的时域信号x_i(t)，i＝1,2,...,M，i为均匀线阵中阵元序号，M为均匀线阵中阵元数目；

步骤2，基于常规波束形成算法(CBF)粗略估计目标方位θ，θ为目标信号波束能量最大时的引导角；步骤3，对步骤2中估计得到的目标方位θ做跟踪波束，得到波束形成后的阵列信号，并对阵列信号作FFT得到阵列信号的频谱图，在频谱图中搜索检测出L根信噪比最大的线谱；

步骤4，计算L根线谱对应的频域值的协方差矩阵R(f_l)，计算L根线谱对应的导向矢量矩阵D(f_l)，计算对应的噪声功率σ_n ²；

步骤5，构造目标函数

通过添加l₁范数和l₂范数约束，将其转化为压缩感知的问题，联合求解L根线谱的空间谱P(θ,f_l)的稀疏解，

为协方差矩阵的估计值，I为单位矩阵；

步骤6，使用经典层次聚类算法对联合估计得到的L根线谱的空间谱进行聚类，将其分为来自目标方位

以及非目标方位

两类。

2.根据权利要求1基于稀疏重构的多线谱联合DOA估计和聚类方法，其特征在于，所述步骤2包括：

步骤2.1，计算理想均匀线阵(ULA)在引导角θ_j下相邻阵元时延τ_j：

其中j＝1,2,...,J+1，j为引导角序号，J+1为总引导角个数，引导角范围为：0～180°，d为相邻阵元间的距离，v是声音在水中的传播速度；

步骤2.2，对各阵元接收的阵列信号进行延时相加，得到目标信号波束能量图B[b(1),…,b(J+1)]，其中b(j)为：

其中，M表示均匀线阵水听器的观测阵列个数，τ_j表示阵元的时延；

步骤2.3，通过对波束能量图的检测搜索找到波束能量最大值位置时的引导角为目标方位的粗估计θ。

3.根据权利要求1所述的一种基于组稀疏的多线谱联合DOA估计和聚类方法，其特征在于，所述步骤3包括：

步骤3.1，根据粗估计的目标方位θ计算各阵元的时延估计

其中v是声音在水中的传播速度，d为相邻阵元间的距离；

步骤3.2，将各阵元数据按时延估计

与参考阵元对齐，将对齐后的阵元数据相干相加获得目标跟踪波束g(t)：

步骤3.3，对g(t)进行傅里叶变换获得目标信号频谱G(f)，同时利用滑动窗平滑技术估计目标信号连续谱G_c(f)，在目标信号频谱G(f)中删除连续谱G_c(f)的影响，获得目标信号的线谱G_line(f)：

G_line(f)＝G(f)-G_c(f)

获得目标的线谱G_line(f)之后，通过计算对应频率f处强线谱的信噪比，通过排序算法得到L根信噪比最大的强线谱，强线谱对应频率为f_l,l＝1,2,3,...,L。

4.根据权利要求1所述的一种基于组稀疏的多线谱联合DOA估计和聚类方法，其特征在于，所述步骤4包括：

步骤4.1，对阵列接收的信号进行分帧处理，每一帧信号的采样点数为S，信号重叠部分长度为S/2，共分为T帧数据；

步骤4.2，对每一帧信号中单个阵元接收信号x_m(t)作傅里叶变换得到x_fft,m(t)，得到L根线谱频率f_l,l＝1,2,...,L处的值x_fft,m(t,f_k)，包含了相位和幅度信息，可以得到每一帧中的线谱阵列向量：

x_fft(t,f_k)＝[x_fft,1(t,f_k),...x_fft,M(t,f_k)]^T；

步骤4.3，利用T帧数据，对得到的L根信噪比最大的线谱，求每一根线谱f_l其频域值的协方差矩阵为：

步骤4.4，计算L根线谱对应的导向矢量矩阵D(f_l)：

上式中矩阵的导向矢量具体表达式为:

步骤4.5，计算噪声功率σ_n ²，用每一根线谱f_l其频域值的协方差矩阵R(f_l)的最小特征值作为噪声功率的估计值

5.根据权利要求1所述的一种基于组稀疏的多线谱联合DOA估计和聚类方法，其特征在于，所述步骤5包括：

步骤5.1，在阵列信号处理中，信源的数目被认为少于传感器阵列的数目，因此L根线谱对应的空间谱分布P(θ,f_l)认为是稀疏的；

步骤5.2，实际的水声信号处理过程中，精确的干扰加噪声协方差矩阵R_i+n不易获取，采用样本的协方差矩阵代替

当采样点数S足够多时，样本的协方差矩阵近似等于干扰加噪声协方差矩阵，即

用L根线谱对应的频域值的协方差矩阵

近似等于精确的干扰加噪声协方差矩阵R_i+n；

步骤5.3，基于上述步骤的说明，构造了多线谱空间谱的联合估计稀疏约束优化问题：

其中，矩阵P(f_l)是频率为f_l强线谱的空间谱分布p(f_l)的对角矩阵，即P(f_l)＝diag{p(f_l)}，

代表频率为f_l的噪声的误差的方差，||·||_F代表矩阵的Frobenius范数，||·||₀代表向量的l₀范数，在这里表示强线谱空间谱中的非0元素个数，系数γ控制着强线谱空间谱的稀疏度和剩余范数之间的权衡；

步骤5.4，步骤5.2中的稀疏约束优化目标函数中由于l₀范数的存在，难以求解，用l₁范数代替l₀范数作为该优化目标函数的近似解，即||p(f_l)||₀＝||p(f_l)||₁；

步骤5.5，对于线谱检测估计得到的L根信噪比最大的线谱，对于来自目标信号方向上的线谱，其空间谱的分布p(f_l)可以认为是近似相等的，即引导角范围0～180°上的值近似相等，因此，引入了组稀疏(Sparse-group)的概念，将0～180°分为181组，第i组数据中包含L根强线谱在引导角θ_i上的功率{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}，对于每一组中的强线谱的功率，本发明中认为该组中的功率要么是非零值，要么是零值，即181组L根强线谱的功率是组稀疏的；

步骤5.6，通过在第i组数据p⁽ⁱ⁾＝{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}中添加l₂范数的形式，实现组稀疏(Sparse-group)；

步骤5.7，对于第i组数据中包含L根强线谱在引导角θ_i上的功率{p(θ_i,f₁),p(θ_i,f₂),...,p(θ_i,f_L)}，对于来自不同方向的强线谱不一定能保证该组中的数据都是非零值或是零值，因此对于每一组中的数据是稀疏的，引入l₁范数代替l₀范数保证每一组中数据的稀疏性，联合估计L根强线谱的空间谱更加接近真实的空间谱分布；

步骤5.8，改进了步骤5.3中优化目标函数，重新构造L根强线谱的组稀疏回归约束(Sparse-group lasso)优化目标函数：

上式中p_i为第i组数据的长度，J+1为总引导角个数，引导角范围为：0～180°，即L根强线谱空间谱分组的组数，α是一个权重系数，α∈[0,1]，α用来控制组稀疏和组内稀疏性的权衡，经过验证，上述优化目标函数为凸函数，可以求解L根强线谱的空间谱分布。

6.根据权利要求1所述的一种基于组稀疏的多线谱联合DOA估计和聚类方法，其特征在于，所述步骤6包括：

步骤6.1，计算L根强线谱空间谱p(f_l)两两之间的Frechet距离，构成矩阵

其中P_ij表示第i个强线谱空间谱分布与第j个强线谱空间谱分布之间的Frechet距离；

步骤6.2，以Frechet距离矩阵P作为距离矩阵，进行层次聚类，对L个强线谱的空间谱分布进行分类得到分类簇

其中c_θ包含来自目标方向的强线谱的空间谱分布，

包含来自非目标方向的强线谱的空间谱分布。

7.一种基于组稀疏的多线谱联合估计DOA系统，其特征在于，包括依次连接的观测阵列信号获取模块、目标方位粗估计模块、线谱检测模块、强线谱相关参数计算模块、强线谱的空间谱联合估计模块、强线谱的空间谱聚类模块；

观测阵列信号获取模块，用于获取观测阵列信号x_i(t)，i＝1,2,...,M，M为均匀线阵中(ULA)阵元数目；

目标方位粗估计模块，用于粗略估计目标方位

为目标信号波束能量最大时的引导角；

线谱检测模块，利用估计得到的目标方位

对阵列信号做跟踪波束想干相加后得到g(t)，对g(t)做傅里叶变换后，利用线谱检测技术检测出L根信噪比最大的线谱；

强线谱相关参数计算模块，对阵列接收的信号进行分帧处理，对每一帧信号中单个阵元接收信号x_m(t)作傅里叶变换得到x_fft,m(t)，求每一根线谱f_l其频域值的协方差矩阵R(f_l)，计算L根线谱对应的导向矢量矩阵D(f_l)，计算噪声功率σ_n ²；

强线谱的空间谱联合估计模块，引入组稀疏的概念，构造了L根强线谱的组稀疏回归约束(Sparse-group lasso)优化目标函数，求解L根强线谱的空间谱分布；

强线谱的空间谱聚类模块，计算L根强线谱空间谱p(f_l)两两之间的Frechet距离，将L根强线谱分为来自目标方向

的强线谱，以及来自非目标方向

的强线谱两类。

Images