基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法

摘要

本发明涉及一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，首先将多帧阵元域回波做波束形成、匹配滤波处理，得到角度‑距离二维声呐图像序列；然后划分子块，通过计算连续几帧图像间的离散系数，得到各子块对应的特征值，并利用特征值重新分配各子块的正则化参数；接着建立BS‑RPCA模型，利用低运算量的非精确拉格朗日乘子法求解目标函数，得到分离后的稀疏矩阵，实现混响抑制。本方法可有效抑制强旁瓣干扰，提高输出信混比。通过分块正则化参数设置提高了低秩稀疏矩阵分解的精度，进一步改善了混响抑制的效果。通过联合当前帧与其之前连续多帧进行处理，可满足声纳系统实时处理的需要。

说明书

技术领域

本发明属于水声探测领域，特别涉及一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法。

背景技术

水下入侵小目标的探测是岛礁、海岸线、港口、内陆河湖防御体系的重要一环。为了保障港口、舰船的安全，通常使用反蛙人声呐对企图潜入港内实施侦察、突袭和破坏的蛙人及无人水下航行器(Unmanned Undersea Vehicles，UUV)等低速小目标实施警戒。反蛙人声呐主要采用主动方式进行探测，但在浅海环境中，海面波浪、卷入水体的气泡等会对小目标探测声呐造成极强的海面和体积混响等杂波干扰，同时水底地形起伏、沉底物等会形成较强的海底混响和杂波干扰，导致反蛙人声呐输出图像中存在着大面积的杂波与干扰，严重影响了声纳的检测性能。为提高反蛙人声纳对低速小目标进行有效探测识别的能力。

混响与目标回波的形成机理存在相似性，在频域上覆盖的区域与发射信号基本重合，时域上与发射信号及目标回波强相关，无法通过常规的时域与频域处理方法得到有效抑制。但由于产生混响的散射体分布位置与目标位置不同，混响与目标回波在空域上不具有相关性。现有文献已经证明，运动的小目标在声呐图像序列间是稀疏的；而对于近岸浅海高混响背景水域区域布置的反蛙人声呐而言，强杂波区域的混响回波在多帧数据间具有一定的低秩结构。因此可利用低秩稀疏矩阵分解将混响与运动的目标分离，在抑制混响的同时提取目标信息。

已有Ge等人利用加速近端梯度法(Accelerated Proximal Gradient，APG)求解低秩稀疏矩阵分解的凸优化问题来实现混响背景中的检测与跟踪(Ge F X,Chen Y,LiW.Target detecton and tracking via structured convex optimization[C]//IcasspIEEE International Conference on Acoustics.IEEE,2017.)。刘冰等人提出基于随机算法(Go Decomposition，GoDec)的抗混响目标探测方法来解决强混响环境下探测移动目标的问题(刘冰,殷敬伟,朱广平,等.基于随机算法的抗混响目标探测方法[J].哈尔滨工程大学学报,2020,041(002):277-281.)。但是相关研究较少关注混响产生机理、混响组成成分等与水声密切相关的先验信息，算法和数据对接过程中重要参数难确定、弱目标易丢失等问题仍然存在。

发明内容

本发明解决的技术问题是：为了解决现有针对水下慢速小目标探测难、易丢失的问题，同时为提高反蛙人声呐的有效探测能力，本发明涉及一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，本方法能在目标回波强度较弱与混响干扰较强的情况下,有效抑制混响，提取运动的小目标信息。本发明具有计算量低、计算快速、通用性强的特点。

本发明的技术方案是：基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，包括以下步骤：

步骤1：获得声纳图像观测矩阵，包括以下子步骤：

子步骤1.1：对主动声呐阵元域接收回波做匹配滤波与常规波束形成处理，幅值归一化得到声呐图像；

子步骤1.2：定义f表示当前帧数，f＝1,2,…，第f帧角度距离二维图像用

子步骤1.3：对多帧结果X作向量化处理，将这些列向量按照帧号顺序组合成一个矩阵D∈R

步骤2：设置分块正则化参数，包括以下步骤：

子步骤2.1：将大小为N

子步骤2.2：计算各子块特征值：

其中，

子步骤2.3：根据得到的各子块特征值得到各子块对应的正则化参数：

步骤3：建立块稀疏模型，包括以下子步骤：

子步骤3.1：建立块稀疏模型：

其中||·||

步骤4：运用非精确拉格朗日乘子法对该模型进行求解，得到稀疏矩阵E的最优解，包括以下子步骤：

步骤4.1：定义增广拉格朗日函数的为：

式中，Y为拉格朗日乘子，μ为一个正标量的惩罚参数；表示矩阵Y与矩阵D-A-E的内积；

步骤4.2：对增广拉格朗日函数中的A、E、Y进行初始化参数设置，E

步骤4.3：固定E

其中S

其中，

步骤4.5：更新Y

Y

步骤4.6：更新μ

μ

步骤4.7：更新迭代次数k＝k+1；

当满足收敛条件||D-A

当不满足收敛条件时，继续下一次循环，即返回到步骤4.3。

步骤5：逆向量化，得到最终输出结果。

将步骤4最终得到的稀疏矩阵E∈R

将F′

本发明进一步的技术方案是：所述子步骤1.2中，最初开始处理的帧数f需要大于等于累积的帧数M，M可取连续的3～10帧。

本发明进一步的技术方案是：所属子步骤2.1中，P×Q为预设的子块大小，且P、Q取值的范围分别为1～3、5～10。

本发明进一步的技术方案是：所属步骤2.2中，在计算各子块特征值时，如果出现ε

本发明进一步的技术方案是：所述步骤4.2中，惩罚参数μ

发明效果

本发明的技术效果在于：本发明提出的基于BS-RPCA的混响抑制方法是一种基于多帧联合处理的抗混响方法，可有效抑制强旁瓣干扰，提高输出信混比。通过分块正则化参数设置提高了低秩稀疏矩阵分解的精度，进一步改善了混响抑制的效果。通过联合当前帧与其之前连续多帧进行处理，可满足声纳系统实时处理的需要。

本方法针对浅海环境下目标强度弱、干扰多的难点，通过混响抑制减小了近距离盲区范围,从而扩大了小目标的有效探测和跟踪距离。

本发明所提的BS-RPCA的混响抑制效果最佳，经本方法处理可降低原始声呐图像的ISLR约35dB，分别较APG、GoDec多降低约25dB、10dB。

附图说明

图1为本发明方法的整体流程示意图。

图2为向量化预处理示意图。

图3为本发明所采用的BS-RPCA算法流程图

图4为第10帧仿真数据处理前后结果图，其中(a)为仿真得到的第10帧声呐图像，(b)、(c)、(d)是分别经APG、GoDec以及本文提出的BS-RPCA方法处理所得的同一帧结果。

图5为分别为原始声呐图像、经APG、GoDec及BS-RPCA方法处理后的声纳图像具体每帧对应的二维积分旁瓣比(ISLR)，单位为dB。

具体实施方式

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

参见图1-图5，本发明的技术方案为：首先将多帧阵元域回波做波束形成、匹配滤波处理，得到角度-距离二维图像序列；然后划分子块，通过计算连续几帧图像间的离散系数，得到各子块对应的特征值，并利用特征值重新分配各子块的正则化参数；接着建立BS-RPCA模型，利用低运算量的非精确拉格朗日乘子法(Inexact Augmented LagrangeMethod，IALM)求解目标函数，得到分离后的稀疏矩阵，实现混响抑制。

本发明的主要内容有：

(1)对阵元域声呐回波做匹配滤波与常规波束形成处理，将获得的波束输出转换为目标的角度距离二维图像。选取连续几帧图像，对其作向量化处理，得到二维观测数据矩阵，将其作为BS-RPCA算法的输入。

(2)RPCA算法中的正则化参数能平衡低秩矩阵与稀疏矩阵间的权重，对混响抑制的性能好坏有着直接的关系。传统的RPCA方法中正则化参数λ为保证目标不丢失，选取固定经验值来算法求解，但这导致结果中也会包含混响中的快变成分。在本发明中，我们将整幅图像划分为大小相同的子块，通过分别设置每个子块的λ值来提高准确性。各子块对应的正则化参数是通过计算各子块张量的离散系数(Coefficient of Variation，CV)进而重新分配λ值的大小来得到的。若CV数值大，则该区域的离散程度大，包含目标的可能性高，λ值相应减小。计算得到的各子块对应的正则化参数作为BS-RPCA算法输入的输入之一。

(3)针对上一步求解的分块正则化参数λ，建立分块参数优化后的BS-RPCA模型。对分离出的前景目标采用块稀疏表示，这是与传统RPCA方法的最大不同之处。利用低运算量的非精确拉格朗日乘子法(Inexact Augmented Lagrange Method)，IALM)求解目标函数，得到分离后的低秩矩阵与稀疏矩阵。

(4)利用计算机数值仿真分别给出了经本发明处理前后的角度-距离二维探测结果，通过计算处理前后积分旁瓣比(Integrated Side Lobe Ratio，ISLR)来定量分析本发明的混响抑制效果。

本发明解决现存问题所采用的技术方案可分为以下5个步骤，整体的流程图如图1所示，下面对每个步骤做出详细说明：

步骤(1)主要关于观测数据矩阵的获得，所涉及的具体内容如下：

步骤1.1：获取声纳图像序列。

首先对主动声呐阵元域接收回波做匹配滤波与常规波束形成处理，将获得的幅值归一化得到目标的角度距离二维图像序列。若用f表示当前帧数，f＝1,2,…。第f帧角度-距离二维图像用

步骤1.2：向量化预处理得到观测矩阵。

选取连续M帧结果构造二维数据观测矩阵。由于主动声呐背景相关性在整个探测过程连续存在，不必将所有获得的数据来做多帧数据联合处理，M可取连续的3～10帧即可。当f≥M时，取第f-M+1到第f帧的结果，记为

步骤(2)主要关于分块正则化参数的设置，所涉及的具体内容如下：

步骤2.1：声纳图像划分子块。

将大小为N

步骤2.2：计算各子块特征值，并由特征值得到各子块对应的正则化参数。

特征值用各分辨单元的离散系数的平方再求和来表示。离散系数定义为方差除均值的平方，各子块的特征值记为ε

其中，

计算整幅图像的平均离散系数，

步骤(3)主要关于目标函数的建立与求解，所涉及的具体内容如下：

对于上一步设置的子块与对应的不同的正则化参数，我们建立BS-RPCA模型求解，来得到精确的目标信息。建立的BS-RPCA模型如下：

式中||·||

式中，Y为拉格朗日乘子，μ为一个正标量的惩罚参数。表示矩阵Y与矩阵D-A-E的内积。采用IALM算法求解式(4)，可得到稀疏矩阵E的最优解。

步骤(4)主要关于利用IALM算法求解增广拉格朗日函数，所涉及的具体内容如下：

步骤4.1：将步骤(1)获得的D与步骤(2)获得的各子块对应的正则化参数λ

步骤4.2：初始化及参数设置。

初始化：E

参数设置：惩罚参数μ

设置迭代次数k＝0。

步骤4.3：固定E

其中S

步骤4.4：固定A

其中，

步骤4.5：更新Y

Y

步骤4.6：更新μ

μ

步骤4.6：更新迭代次数k＝k+1。

当满足收敛条件||D-A

利用IALM算法求解BS-RPCA问题的流程如图3所示。

步骤(5)主要关于稀疏矩阵的后处理，所涉及的具体内容如下：

通过步骤(4)得到稀疏矩阵E∈R

下面结合附图和具体实施例对本发明进行进一步解释说明：

首先建立一个混响模型，利用计算机进行数值仿真，检验本发明所提方法。定义R[f]为角度-距离二维图空间中的第f帧混响序列，则各帧混响可表示为

R[f]＝βR

式中：β是一个用来控制混响变化的剧烈程度；R

若用A

A

式中：0<α<1是一个随机变量，控制快变混响变化的剧烈程度。

在生成的混响数据中添加一个移动目标，目标强度表示为t，用矩阵T

使用上述混响模型生成20帧数据。其中，N

使用该仿真结果对本发明所提的方法进行验证，实验结果如图4所示。图4(a)为仿真数据未作处理的结果，可以看到目标被淹没在混响之中，几乎不能被分辨出来。图4(b)、4(c)、4(d)是分别经APG、GoDec以及本文提出的BS-RPCA方法处理所得的同一帧结果。参见图4，仿真中我们设置目标从(20°，50m)位置开始以每帧移动1°、5m速度向远处(40°，150m)处移动。目标所在的位置采用白色圆圈表示。图4(a)为仿真数据未作处理的结果，可以看到目标被淹没在混响之中，几乎不能被分辨出来。图4(b)、4(c)、4(d)是分别经APG、GoDec以及本文提出的BS-RPCA方法处理所得的同一帧结果。通过图4我们可看出：GoDec算法的效果较APG方法更好一些，BS-RPCA方法的效果优于其余两种方法，证明本发明所提方法确实可提高低秩稀疏矩阵分解的精度，达到更好的混响抑制效果。

此外，在上述仿真条件下，采用的实验系统为Intel Core i7-8565U、主频为1.80GHz的CPU，运算内存为8G的Windows平台。利用软件Matlab R2014b中进行运算，每帧处理所需时间为0.25s左右，满足实时运算处理的需要。

采用雷达成像指标中的积分旁瓣比ISLR来定量分析本发明方法混响抑制的效果。积分旁瓣比的定义为全体旁瓣能量与主瓣能量的比值，用来衡量局部图像对比度，表现为弱信号区域受强信号旁瓣的干扰能力。积分旁瓣比越小，图像质量越高。积分旁瓣比ISLR定义式为

式中，E