基于稀疏分解理论的声矢量阵宽带测向方法

摘要

本发明公开了基于稀疏分解理论的声矢量阵宽带测向方法，它涉及水声矢量传感器阵探测领域中针对目标宽带连续谱噪声的一种测向技术及实现方法。本发明将稀疏化思想引入声矢量阵测向系统，根据宽带连续谱噪声信号带宽内各频点分量建立一种联合稀疏约束，最终得到统一的空域稀疏分解形式，从而实现对宽带连续谱噪声的测向。本方法能够形成较为尖锐的谱峰和幅度较低的噪底；在阵列快拍数较低时性能稳定，且能分辨相干信号源；在目标来波方向接近阵列轴向时，也不存在常规的波束形成(CBF)法和最小方差无失真响应(MVDR)法在低信噪比下出现的噪底起伏现象。

说明书

技术领域

本发明属于水声信号处理领域，特别涉及一种基于稀疏分解理论的声矢量阵宽带连续谱噪声测向方法，可用于被动水声矢量传感器阵列探测目标辐射的宽带连续谱噪声信号。

背景技术

对水面水下各类目标自身辐射的噪声包括窄带线谱噪声和宽带连续谱噪声。在被动声纳窄带线谱探测方面，传统方法主要通过空域的常规波束合成(CBF)处理获得空间增益，从而完成测向、检测等任务。而对于宽带连续谱噪声信号，CBF方法不能直接使用，通常的做法是将宽带信号拆分成窄带信号，再利用CBF进行测向。这种处理思路没有充分利用宽带信号的信息，它的方位分辨能力受瑞利限限制，不能有效分辨位于主瓣内的两个目标。最小方差无失真响应(MVDR)法在阵列快拍数较低时的目标方位估计偏差较大。上述两种方法在低信噪比下会出现噪底起伏现象。

基于稀疏分解理论的空间谱估计是空间谱估计的另一个方向。信号的稀疏分解理论于上世纪90年代提出，广泛应用于信号压缩、识别、提取等领域，其特有的信号稀疏性理论适用于信号在空间内的分布情况，因此，信号的稀疏分解理论及其算法为声矢量阵列的空间谱估计提供了一个新的解决思路和方向。

发明内容

本发明所要解决的技术问题在于提供一种基于稀疏分解理论的声矢量阵宽带测向方法。本发明将稀疏化思想引入声矢量阵测向系统中，并基于多重分片(Multiple Measurement Slice)的正交匹配追踪(MMS-OMP)算法，将宽带信号带宽内各频点分量建立一种联合稀疏约束，最终得到统一的空域稀疏分解形式，从而实现对宽带连续谱噪声信号的测向。

本发明所要解决的技术问题是由以下技术方案实现的：

基于稀疏分解理论的声矢量阵宽带测向方法，采用分片稀疏表示方法，实现对声矢量阵接收的数据分片稀疏表示，通过联合稀疏约束求解该稀疏表示模型，对求解得到的稀疏矩阵进行空间谱计算，从而实现目标信号的测向。具体包括以下步骤：

(1)将声矢量阵在预设时间段内的接收信号进行频域建模，得到接收信号在多个预定频率处的频域快拍；

(2)计算每个频域快拍对应的协方差矩阵；

(3)根据空域扫描范围和宽带连续谱噪声的带宽建立过完备原子库，并对协方差矩阵进行稀疏表示；

(4)将稀疏表示后的所有协方差矩阵的同一列按频率前后顺序堆叠形成三维张量协方差矩阵，同时，过完备原子库也按照频率前后顺序堆叠形成三维张量过完备原子库后，将三维张量协方差矩阵用三维张量过完备原子库和三维系数矩阵来表示；所述的三维系数矩阵为未知量；

(5)分别初始化信号支撑集合、已选原子集合和迭代次数；

(6)计算过完备原子库中每个原子与当前残差信号的内积，利用准则公式对内积进行处理得到过完备原子库中每个原子与当前残差信号的近似程度值；所述的残差信号的初始值为三维张量协方差矩阵；

(7)选择最大的近似程度值，根据最大的近似程度值找出对应原子的索引值，并更新信号支撑集合和已选原子集合；

(8)利用最小二乘法计算更新后的已选原子集合在当前残差信号方向上的正交投影后，根据正交投影和当前残差信号计算下一次残差信号；

(9)判定更新后的信号支撑集合个数是否大于迭代次数或者当前残差信号的2-范数是否小于预设阈值，如果更新后的信号支撑集合个数大于迭代次数或者当前残差信号的二范数小于预设阈值，则执行步骤(10)；否则，将下一次残差信号更新为当前残差信号，跳转到步骤(6)；

(10)根据最终的正交投影和最终的已选原子集合计算得到三维系数矩阵；根据最终的信号支撑集合和宽带连续谱噪声的来波方向的对应关系，计算三维系数矩阵的F-范数，并根据F-范数得到宽带连续谱噪声的空间谱估计；

(11)对空间谱估计进行谱峰搜索得到前K个空间角度值；所述的K为迭代次数。

其中，步骤(3)所述的根据空域扫描范围和宽带连续谱噪声的带宽建立过完备原子库，具体包括步骤：

(301)将协方差矩阵中的每一列表示成

$\begin{array}{c}{r}_n\left(f_m\right)=A\left(f_m\right)R_x\left(f_m\right)A^H\left(f_m\right)\\ =A\left(f_m\right)s\left(f_m\right)\end{array};$

其中，R_x(f_m)为协方差矩阵；f_m为频率，且f_m∈[f₁,f_M]，M为宽带连续谱噪声带宽内的频点数；A(f_m)为导向矢量矩阵，其每一列对应一个来波方向；s(f_m)为宽带连续谱噪声；A^H(f_m)为导向矢量矩阵的共轭转置；r_n(f_m)为协方差矩阵中的某一列，n＝1,2,…,N，N为声矢量阵的阵元数；

(302)将-180度～180度的空间进行均匀划分，构造过完备化的阵列流型矩阵，满足

r_n(f_m)＝G(f_m)f_n(f_m)；

其中，f_n(f_m)为系数矩阵，其非零行对应真实的信号来波方向，其它非真实来波方向的数据均为0；G(f_m)为过完备化的阵列流型矩阵即过完备原子库。

其中，所述的步骤(4)具体包括以下步骤：

(401)将稀疏表示后的不同频率的协方差矩阵，按频率前后顺序堆叠成一个三维张量协方差矩阵；

(402)将协方差矩阵中的每一列放入稀疏表示后的协方差矩阵对应的位置，同时，将过完备原子库和系数矩阵也按频率前后顺序堆叠存放形成三维张量过完备原子库与三维系数矩阵，并且满足

其中，为三维张量协方差矩阵；为三维张量过完备原子库；为三维系数矩阵，它的每一水平切片都具有相同的稀疏结构。

其中，步骤(6)所述的利用准则公式对内积进行处理得到过完备原子库中每个原子与当前残差信号的近似程度值；具体为：

第k次迭代时，原子g_i与当前残差信号的近似程度值为：

${\Delta }_{i,k}=\underset{n=1}{\overset{3N}{\Sigma }}\underset{m=1}{\overset{M}{\Sigma }}\left|\right|g_i^H\left(f_m\right)R_n^{(k-1)}\left(f_m\right)|{|}_2;$

其中，Δ_i,k为近似程度值；N为声矢量阵的阵元数；M为宽带连续谱噪声带宽内的频点数；||·||₂为2-范数；g_i(f_m)为三维张量过完备原子库中对应角度θ_i的垂直分片中对应频率f_m的列向量；为当前残差信号中第n个垂直切分片中对应频率f_m的列向量；i为三维张量过完备原子库中原子的索引值，i∈I，I为过完备原子库的原子的索引值个数。

其中，所述的步骤(7)具体为：挑选符合条件的原子g补充已选原子集合，选择条件为其中，Q^(k)为已选原子集合；

更新后的信号支撑集合Ω^(k)＝Ω^(k-1)∪{i_k}，更新后的已选原子集合其中，i_k为第k次迭代时，三维张量过完备原子库中原子的索引值；为第k次迭代时，三维张量过完备原子库的第i_k列数据。

其中，所述的步骤(8)具体包括以下步骤：

(801)利用最小二乘法计算更新后的已选原子集合在当前残差信号方向上的正交投影矩阵，计算公式为：

$P_{Q\left(f\right)}=Q_f{Q}_f^{+}$

其中，P_Q(f)为正交投影矩阵；Q_f为更新后的已选原子集合，Q_f＝[g_i1>i2>

(802)根据正交投影矩阵和当前残差信号，计算下一次残差信号，计算方法为：

$R_x^{\left(k\right)}\left(f_m\right)=R_x^{(k-1)}\left(f_m\right)(I-P_{Q\left(f\right)});$

其中，为当前残差信号中第x个水平切分片中对应频率f_m的列向量，m＝1,2,…,M；I为单位矩阵。

其中，步骤(10)所述的根据最终的信号支撑集合和宽带连续谱噪声来波方向的对应关系，计算三维系数矩阵的F-范数，并根据F-范数得到宽带连续谱噪声的空间谱估计；具体包括以下步骤：

(111)将三维系数矩阵按三维张量过完备原子库中原子的索引值划分成列片结构；

(112)根据列片结构划分，计算三维系数矩阵的功率谱估计；计算公式为：

$P_{mms-omp}\left({\theta }_i\right)=\left|\right|F_{\theta }\left(i\right)|{|}_F^2$

其中，P_mms-omp(θ_i)为三维系数矩阵的功率谱估计；i为三维张量过完备原子库中原子的索引值；θ_i为索引值i对应的来波方向，当I个原子组成的过完备原子库对应的来波方向在整个空域均匀分布时，有||·||_F表示F-范数；F_θ(i)为每一列片结构的数据矩阵。

本发明与现有技术相比具有如下优点：

1、本发明与传统方法相比，具有尖锐的谱峰，具有更低的旁瓣水平，更有利于信号的分辨和检测。

2、本发明单次快拍即可完成对目标辐射的宽带连续谱信号进行快速测向，测向精度高。

3、本发明能分辨相干信号源，在目标来波方向接近阵列轴向时，也不存在传统常规方法的噪底起伏现象。

附图说明

图1是本发明的声矢量阵列接收信号原理图；

图2是本发明的声矢量阵宽带测向方法流程图；

图3是本发明的三维张量协方差矩阵三维系数矩阵和三维的构造方式示意图；

图4是本发明三维张量过完备原子库的列片划分示意图；

图5是本发明三维张量协方差矩阵的稀疏分解示意图；

图6是本发明三维张量协方差矩阵的列片分解示意图；

图7是本发明三维系数矩阵的稀疏分片结构示意图；

图8是使用本发明方法声矢量阵宽带处理结果；

图9是本发明方法声矢量阵宽带处理与其它方法的处理结果对比；

图10是本发明方法估计的角度方向的RMSE(均方根误差)与其它方法的对比。

具体实施方式

下面结合附图，对本发明的技术方案和效果作进一步详细说明。

由于舰艇辐射噪声的宽带特性，在将稀疏分解理论应用到矢量水听器阵列测向中，必然会产生宽带信号的稀疏分解问题。宽带信号频率能量分布的不均匀性，将导致频带内各频点的空域稀疏特性不相同，若对各频点进行独立的窄带空域稀疏分解，再将分解结果简单相加，不仅损失了信噪比，也不能处理相干信号。为了得到带宽内统一的空域稀疏分解，需要一种方式能够将各频点的稀疏结构进行综合。多重分片处理是能够对多个平行子分解进行联合统一约束的一种处理方式。基于MMS-OMP算法将窄带信号的稀疏分解扩展为多个窄带信号的联合统一稀疏分解。通过MMS-OMP算法，多个窄带信号可以得到一个统一的稀疏表示形式，从而得到宽带信号的稀疏分解。在基于MMS-OMP的矢量阵宽带测向算法中，根据宽带信号带宽内各频点分量建立带宽内的联合约束，最终得到统一的空域稀疏分解形式，从而确定宽带信号的来波方向。

图1是本发明的声矢量阵列接收信号原理图。

图中声矢量阵列共有N个阵元，按照等间距排列，阵元间距为d，N个阵元对应有N个接收通道。每个接收通道接收到水声宽带连续谱噪声信号以后经过滤波放大调理和AD变换之后，输出一路数字信号以供后续处理。

本发明的基于稀疏分解理论的声矢量阵宽带测向方法，具体流程如图2所示，实现步骤如下：

步骤1：参考图1，对声矢量阵在预设时间段内的接收信号进行频域模型建模。

(101)考虑由K个目标发出或反射的宽带连续谱噪声信号以波速c经水声信道传播，在目标距离声矢量阵足够远的条件下，宽带连续谱噪声信号近似为远场平面波，并以来波方向θ_i(i＝1,2,...,K)入射至声矢量阵。若上述宽带连续谱噪声信号分布在频率f₁～f_M之间，则声矢量阵在时间段T内的接收信号在频域的数学模型可表示为：

X＝[x(f₁)>2) ...>M)]

其中，x(f_m)(m＝1,2,...,M)为声矢量阵在频率f_m处的频域快拍。

(102)在基于扩展法的声矢量阵处理中，将x(f_m)表示成如下形式：

$\left(\begin{array}{c}x\left(f_m\right)={\left(\begin{array}{cccc}{x}_1\left(f_m\right)& x_2\left(f_m\right)& \mathrm{...}& x_{3N}\left(f_m\right)\end{array}\right)}^T\\ =\underset{i=1}{\overset{K}{\Sigma }}a({\theta }_i,f_m)s_i\left(f_m\right)\otimes u\left({\theta }_i\right)+N\left(f_m\right)\end{array}\right)$

其中，s(f_m)为声矢量阵接收的目标声压信号在频率f_m处的分量；a(θ,f_m)为目标在频率f_m处对应的导向矢量，对于第i个目标，有其中ω_m＝2πf_m，τ_i＝dcosθ_i/c；u(θ)＝[cosθ>T为声矢量传感器方向向量；符号为Kronecker积；N(f_m)＝[n₁(f_m)>2(f_m) ...>3N(f_m)]^T为空间加性噪声矢量。

(103)将上式改写为矩阵乘积的形式为：

x(f_m)＝A(f_m)S(f_m)+N(f_m)

其中，为导向矢量矩阵；S(f_m)＝[s₁(f_m)>2(f_m) ...>K(f_m)]^T为信号在频率f_m处的分量。

步骤2：计算频率f_m对应的频域快拍的协方差矩阵。

x(f_m)的协方差矩阵R_x(f_m)为：

R_x(f_m)＝E[x(f_m)x^H(f_m)]

＝A(f_m)E[S(f_m)S^H(f_m)]A^H(f_m)+R_n(f_m)

＝A(f_m)R_s(f_m)A^H(f_m)+R_n(f_m)

由前述数学模型可知，矩阵A(f_m)中含有信号在频率f_m处分量的来波方向信息，对宽带信号的来波方向估计即是对宽带信号所在的各频点处分量来波方向的综合估计。

步骤3：根据空域扫描范围和宽带连续谱噪声的带宽建立过完备原子库G(f_m)，并对协方差矩阵进行稀疏表示。

建立过完备原子库G(f_m)的过程：

(301)首先，在忽略声矢量阵列接收噪声的基础上，将R_x(f_m)中的每一列表示成形式：

r_n(f_m)＝A(f_m)R_x(f_m)A^H(f_m)

＝A(f_m)s(f_m)

其中，A(f_m)为导向矢量矩阵，它中的每一列对应一个来波方向。

(302)将{θ_0°～θ_180°}角度空间完备化，构造一个完备化的阵列流型矩阵G(f_m)(即过完备原子库)，满足

r_n(f_m)＝G(f_m)f_n(f_m)n＝1,2,…,N

其中，f_n(f_m)中非零行对应真实的信号来波方向，其它非真实来波方向的数据均为0。

在尽可能充分利用信号带宽的前提下，对宽带信号的稀疏表示需在频带[f₁f_M]内进行。

步骤4：对稀疏表示后的协方差矩阵R_x(f_m)和过完备原子库分别按频率前后顺序堆叠形成三维张量协方差矩阵和三维张量过完备原子库后，将三维张量协方差矩阵用三维张量过完备原子库和三维系数矩阵来表示。

(401)将稀疏表示后的不同频率的协方差矩阵R_x(f_m)，f_m∈[f₁f_M]按频率顺序堆叠成一个三维张量协方差矩阵

(402)将r_n(f_m)放入R_x(f_m)对应的位置时，过完备原子库G(f_m)和系数矩阵f(f_m)也按频率顺序堆叠存放形成三维张量过完备原子库与三维系数矩阵并且满足

其中，三维系数矩阵为未知量，它的每一水平切片F(f_m)都具有相同的稀疏结构。

参考图3、4完成数据构造，参考图5、6对三维张量协方差矩阵进行稀疏表示，建立分片稀疏表示模型。

步骤5：初始化参数设置，包括信号支撑集合、已选原子集合和迭代次数。

初始化信号支撑集合已选原子集合迭代次数K等参数。

步骤6：寻找与当前残差信号最为匹配的原子：计算过完备原子库G(f_m)中每个原子与当前残差信号的内积并利用准则公式得到过完备原子库G(f_m)中每个原子与当前残差信号的近似程度值Δ_i,k。

设三维张量过完备原子库为第k次迭代时，定义原子g_i与当前残差信号的近似程度值为：

${\Delta }_{i,k}=\underset{n=1}{\overset{3N}{\Sigma }}\underset{m=1}{\overset{M}{\Sigma }}\left|\right|g_i^H\left(f_m\right)R_n^{(k-1)}\left(f_m\right)|{|}_2$

其中，||·||₂为2-范数；N为阵元数；M为带宽内频点数；为当前残差信号中第n个垂直切分片中对应频率f_m的列向量；g_i(f_m)为三维张量过完备原子库中对应角度θ_i的垂直分片中对应频率f_m的列向量。

步骤7：选择最大Δ_i,k值对应原子的索引值，根据索引值更新信号支撑集合和已选原子集合。

挑选符合条件的原子g补充原子集，选择条件为：

$g={\arg }_{g_i}max\left({\Delta }_{i,k}\right)$

更新信号支撑集合Ω^(k)＝Ω^(k-1)∪{i_k}，更新已选原子集合其中，i_k为第k次迭代时，三维张量过完备原子库中原子的索引值；为第k次迭代时，三维张量过完备原子库的第i_k列数据。

步骤8：通过最小二乘法计算更新残差信号：利用最小二乘法计算更新后的已选原子集合在当前残差信号方向上的正交投影P_Q(f)，进而重新组合得到残差矩阵完成更新残差矩阵解算。

(801)利用最小二乘法计算更新后的已选原子集合在当前残差信号方向上的正交投影矩阵P_Q(f)，计算公式为

$P_{Q\left(f\right)}=Q_f{Q}_f^{+}$

其中，Q_f为已选原子集合Q_f＝[g_i1>i2>

(802)根据正交投影矩阵和当前残差信号，计算下一次残差信号计算方法为：

$R_x^{\left(k\right)}\left(f_m\right)=R_x^{(k-1)}\left(f_m\right)(I-P_{Q\left(f\right)}),m=1,2,...,M$

其中，为当前残差信号中第x个水平切分片中对应频率f_m的列向量，m＝1,2,…,M；I为单位矩阵。

步骤9：判定是否终止迭代：根据更新后的信号支撑集合个数是否满足设定的迭代次数K或者通过判断当前残差信号的2-范数是否小于预设阈值，来判定是否终止步骤(6)～(8)的迭代。

步骤10：根据最终的正交投影和最终的已选原子集合计算得到三维系数矩阵；根据最终的信号支撑集合和宽带连续谱噪声的来波方向的对应关系，计算三维系数矩阵的Frobenius范数(简称F-范数)，并根据F-范数得到宽带连续谱噪声的空间谱估计；

(111)由于三维系数矩阵的每一水平分片F(f_m)都具有相同的行稀疏结构，因此将三维系数矩阵按三维张量过完备原子库中原子索引值i划分成列片结构，每一列片数据矩阵表示为F_θ(i)；

(112)根据列片划分，计算三维系数矩阵的功率谱估计，计算公式为：

$P_{mms-omp}\left({\theta }_i\right)=\left|\right|F_{\theta }\left(i\right)|{|}_F^2$

其中，||·||_F表示F-范数；θ_i为索引值i对应的来波方向。当I个原子组成的过完备原子库对应的来波方向在整个空域均匀分布时，有

参考图7，计算阵列接收数据的空间谱。经过K次分解后得到三维系数矩阵利用原子集索引值与噪声信号来波方向具有对应关系，通过计算相应的原子集索引对应系数矩阵的F-范数得到空间谱估计P_mms-omp(θ_i)。

步骤11：对空间谱估计进行谱峰搜索得到前K个空间角度值；所述的K为迭代次数。

本发明的效果可以通过以下仿真说明：

1.仿真条件与方法

设有一16元声矢量阵接收带宽分别为150Hz～170Hz和180Hz～200Hz的两个目标信号，以阵列轴向为0°，其来波方向分别为78°和-146°，接收机带宽为140Hz～210Hz，带内信噪比为0dB。

2.仿真内容与结果

图8是将迭代次数限制为2次的稀疏分解结果，迭代次数等价于已知宽带信号包含的目标数，当目标数未知时，可采用残差能量阈值来迭代次数，但是有可能会增加不必要的计算量。

在不改变前述仿真参数的条件下，将本文提出的基于MMS-OMP的矢量阵宽带测向算法与常用的宽带非相干MVDR、非相干CBF和STMV算法进行比较，可得如图9所示的结果。本发明的方法能够形成较为尖锐的谱峰和幅度较低的噪底，这也是基于稀疏分解的阵列测向算法的一个特点，在目标来波方向接近阵列轴向时，也不存在上述三种算法的噪底起伏现象。

不改变其余仿真参数，对本文提出的算法在不同带内信噪比下的测向误差进行仿真，可得如图10所示的结果。在不同带内信噪比下，本发明的矢量阵宽带测向方法的测向均方根误差与STMV算法相近，其相对较窄的目标谱峰是与常规空间谱估计方法相比的优势所在。综上所述，本发明的声矢量阵宽带信号测向方法能够形成较为尖锐的谱峰和幅度较低的噪底。相比传统方法增强了目标信号的检测和分辨能力，在阵列快拍数较低时性能较为稳定。在目标来波方向接近阵列轴向时，有效解决了上述传统方法的噪底起伏现象。