CN113050098B - 基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法 - Google Patents
基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法 Download PDFInfo
- Publication number
- CN113050098B CN113050098B CN202110249031.0A CN202110249031A CN113050098B CN 113050098 B CN113050098 B CN 113050098B CN 202110249031 A CN202110249031 A CN 202110249031A CN 113050098 B CN113050098 B CN 113050098B
- Authority
- CN
- China
- Prior art keywords
- matrix
- block
- frame
- sparse
- sonar
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 52
- 230000001629 suppression Effects 0.000 title claims abstract description 23
- 238000000513 principal component analysis Methods 0.000 title claims abstract description 11
- 239000011159 matrix material Substances 0.000 claims abstract description 66
- 238000012545 processing Methods 0.000 claims abstract description 20
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 11
- 238000001914 filtration Methods 0.000 claims abstract description 6
- 230000008569 process Effects 0.000 claims description 7
- 230000003190 augmentative effect Effects 0.000 claims description 5
- 239000013598 vector Substances 0.000 claims description 4
- 230000008602 contraction Effects 0.000 claims description 3
- 238000013139 quantization Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000003595 spectral effect Effects 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 8
- 238000002592 echocardiography Methods 0.000 abstract description 2
- 238000002156 mixing Methods 0.000 abstract description 2
- 238000004422 calculation algorithm Methods 0.000 description 15
- 238000001514 detection method Methods 0.000 description 13
- 238000004088 simulation Methods 0.000 description 9
- 230000006870 function Effects 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 238000005457 optimization Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000000638 solvent extraction Methods 0.000 description 2
- 241000287196 Asthenes Species 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 210000001503 joint Anatomy 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
- G01S15/88—Sonar systems specially adapted for specific applications
- G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
- G01S7/539—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computer Networks & Wireless Communication (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Acoustics & Sound (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明涉及一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，首先将多帧阵元域回波做波束形成、匹配滤波处理，得到角度‑距离二维声呐图像序列；然后划分子块，通过计算连续几帧图像间的离散系数，得到各子块对应的特征值，并利用特征值重新分配各子块的正则化参数；接着建立BS‑RPCA模型，利用低运算量的非精确拉格朗日乘子法求解目标函数，得到分离后的稀疏矩阵，实现混响抑制。本方法可有效抑制强旁瓣干扰，提高输出信混比。通过分块正则化参数设置提高了低秩稀疏矩阵分解的精度，进一步改善了混响抑制的效果。通过联合当前帧与其之前连续多帧进行处理，可满足声纳系统实时处理的需要。
Description
技术领域
本发明属于水声探测领域，特别涉及一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法。
背景技术
水下入侵小目标的探测是岛礁、海岸线、港口、内陆河湖防御体系的重要一环。为了保障港口、舰船的安全，通常使用反蛙人声呐对企图潜入港内实施侦察、突袭和破坏的蛙人及无人水下航行器(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帧角度距离二维图像用表示，其中Nθ表示波束数，Nt时间采样点数；选取连续M帧结果构造二维数据观测矩阵；当f≥M时，取第f-M+1到第f帧的结果，记为张量/>则X＝(F[f-M+1],F[f-M+2],...,F[f])，其中F[·]表示括号内帧数的声呐图像；
子步骤1.3：对多帧结果X作向量化处理，将这些列向量按照帧号顺序组合成一个矩阵D∈RN×M，其中N＝Nθ×Nr，M为选取的帧数；
步骤2：设置分块正则化参数，包括以下步骤：
子步骤2.1：将大小为Nr×Nθ的声呐图像分为m×n个子块，子块大小为P×Q，P、Q为初始设定值；记WINi,j表示第i行、第j列子块的位置，Gi,j(p,q,m)表示在张量X中WINi,j内(p,q,m)位置的数值，其中p＝1,2,…,P，q＝1,2,…,Q，m＝1,2,…,M。
子步骤2.2：计算各子块特征值：
其中，表示遍历(p,q)对所有元素求和；
子步骤2.3：根据得到的各子块特征值得到各子块对应的正则化参数：
步骤3：建立块稀疏模型，包括以下子步骤：
子步骤3.1：建立块稀疏模型：
其中||·||F为矩阵的Frobenius范数，Pi,j表示将E中的每一列还原成矩阵的操作，返回WINi,j位置稀疏块在每一帧图像内的矩阵表示；
步骤4：运用非精确拉格朗日乘子法对该模型进行求解，得到稀疏矩阵E的最优解，包括以下子步骤：
步骤4.1：定义增广拉格朗日函数的为：
式中，Y为拉格朗日乘子，μ为一个正标量的惩罚参数；<Y,D-A-E>表示矩阵Y与矩阵D-A-E的内积；
步骤4.2：对增广拉格朗日函数中的A、E、Y进行初始化参数设置，E0＝0，A0＝0，Y0＝D/max(||D||2,λ0 -1||D||∞)，其中||·||2为矩阵的谱范数，称为行和范数；定义惩罚参数μ0，迭代步长ρ，最大迭代次数为kmax，设置迭代次数k＝0；
步骤4.3：固定Ek、Yk，更新低秩矩阵Ak+1
其中Sε[·]表示常规软阈值算子，步骤4.4：固定Ak+1、Yk，更新稀疏矩阵Ek+1
表示块收缩算子，/>
其中，
步骤4.5：更新Yk+1
Yk+1＝Yk+μk(D-Ak+1-Ek+1)
步骤4.6：更新μk+1
μk+1＝min(ρμk,107μ0)
步骤4.7：更新迭代次数k＝k+1；
当满足收敛条件||D-Ak-Ek||F/||D||F<10-7或达到最大迭代次数kmax时，输出此时的Ak+1、Ek+1，此时的Ak即为最终分解得到的低秩矩阵，属于被抑制掉的混响部分；Ek即为最终分解得到的稀疏矩阵，包含运动的目标部分，是需要的部分。
当不满足收敛条件时，继续下一次循环，即返回到步骤4.3。
步骤5：逆向量化，得到最终输出结果。
将步骤4最终得到的稀疏矩阵E∈RN×M中的每列重组为原始单帧数据维度，用F′表示逆向量化后的结果，可得到张量[F′1,F′2,...,F′M]，F′M即为当前f帧输出的最终结果。至此，我们完成了第f帧数据的混响抑制处理，通过低秩稀疏矩阵分解在获得运动目标信息的同时抑制了慢变的混响成分。
将F′M输出后返回到步骤1进行第f+1帧数据处理，继续步骤1～5。
本发明进一步的技术方案是：所述子步骤1.2中，最初开始处理的帧数f需要大于等于累积的帧数M，M可取连续的3～10帧。
本发明进一步的技术方案是：所属子步骤2.1中，P×Q为预设的子块大小，且P、Q取值的范围分别为1～3、5～10。
本发明进一步的技术方案是：所属步骤2.2中，在计算各子块特征值时，如果出现εi,j＝0，则步骤2.3中λi,j直接记为λ0，λ0的取值建议为其中N＝Nθ×Nr。
本发明进一步的技术方案是：所述步骤4.2中，惩罚参数μ0需为正值，迭代步长ρ需大于1；最大迭代次数kmax设置为10次左右。
发明效果
本发明的技术效果在于：本发明提出的基于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帧角度-距离二维图像用表示，其中Nθ表示波束数，Nt时间采样点数。
步骤1.2：向量化预处理得到观测矩阵。
选取连续M帧结果构造二维数据观测矩阵。由于主动声呐背景相关性在整个探测过程连续存在，不必将所有获得的数据来做多帧数据联合处理，M可取连续的3～10帧即可。当f≥M时，取第f-M+1到第f帧的结果，记为则X＝(F[f-M+1],F[f-M+2],...,F[f])，其中F[·]表示括号内帧数的角度-距离二维图像。对多帧结果X作向量化处理，即每一帧图像按行拉伸为一个列向量，将这些列向量按照帧号顺序组合成一个矩阵D∈RN×M，其中N＝Nθ×Nr，M为选取的帧数。向量化预处理步骤示意图如图2所示。
步骤(2)主要关于分块正则化参数的设置，所涉及的具体内容如下：
步骤2.1：声纳图像划分子块。
将大小为Nr×Nθ的声纳图像分为m×n个子块，子块大小为P×Q，P、Q为预设。记WINi,j表示第i行、第j列子块的位置，Gi,j(p,q,m)表示在张量X中WINi,j内(p,q,m)位置的数值，其中p＝1,2,…,P，q＝1,2,…,Q，m＝1,2,…,M。
步骤2.2：计算各子块特征值，并由特征值得到各子块对应的正则化参数。
特征值用各分辨单元的离散系数的平方再求和来表示。离散系数定义为方差除均值的平方，各子块的特征值记为εi,j，
其中，表示遍历(p,q)对所有元素求和。
计算整幅图像的平均离散系数， 表示遍历(i,j)所有元素求平均。当前帧数为f时，记WINi,j位置对应的正则化参数为λi,j，
步骤(3)主要关于目标函数的建立与求解，所涉及的具体内容如下：
对于上一步设置的子块与对应的不同的正则化参数，我们建立BS-RPCA模型求解，来得到精确的目标信息。建立的BS-RPCA模型如下：
式中||·||F为矩阵的Frobenius范数，即矩阵元素的平方和开根号；Pi,j表示将E中的每一列还原成矩阵的操作，返回WINi,j位置稀疏块在每一帧图像内的矩阵表示。这里对前景目标的块稀疏表示，也是与传统RPCA方法的最大不同之处。上式仍是一个凸优化问题，出于计算速度与计算精度的考虑，可通过IALM算法求解，其增广拉格朗日函数的定义为:
式中，Y为拉格朗日乘子，μ为一个正标量的惩罚参数。<Y,D-A-E>表示矩阵Y与矩阵D-A-E的内积。采用IALM算法求解式(4)，可得到稀疏矩阵E的最优解。
步骤(4)主要关于利用IALM算法求解增广拉格朗日函数，所涉及的具体内容如下：
步骤4.1：将步骤(1)获得的D与步骤(2)获得的各子块对应的正则化参数λi,j作为IALM算法的输入。
步骤4.2：初始化及参数设置。
初始化：E0＝0，A0＝0，Y0＝D/max(||D||2,λ0 -1||D||∞)，其中||·||2为矩阵的谱范数，即矩阵中最大的正奇异值，(i＝1,2,…)称为行和范数。
参数设置：惩罚参数μ0＞0，迭代步长ρ＞1，最大迭代次数为kmax。
设置迭代次数k＝0。
步骤4.3：固定Ek、Yk，更新低秩矩阵Ak+1
其中Sε[·]表示常规软阈值算子，
步骤4.4：固定Ak+1、Yk，更新稀疏矩阵Ek+1
表示块收缩算子，/>
其中，
步骤4.5：更新Yk+1
Yk+1＝Yk+μk(D-Ak+1-Ek+1) (8)
步骤4.6：更新μk+1
μk+1＝min(ρμk,107μ0) (9)
步骤4.6：更新迭代次数k＝k+1。
当满足收敛条件||D-Ak-Ek||F/||D||F<10-7或达到最大迭代次数kmax时，输出此时的Ak+1、Ek+1。
利用IALM算法求解BS-RPCA问题的流程如图3所示。
步骤(5)主要关于稀疏矩阵的后处理，所涉及的具体内容如下：
通过步骤(4)得到稀疏矩阵E∈RN×M，将稀疏矩阵E中的每列重组为原始单帧数据维度，用F′表示逆向量化后的结果，可得到张量[F′1,F′2,...,F′M]，F′M即为当前f帧输出的最终结果。输出结果后进行第f+1帧数据处理。
下面结合附图和具体实施例对本发明进行进一步解释说明：
首先建立一个混响模型，利用计算机进行数值仿真，检验本发明所提方法。定义R[f]为角度-距离二维图空间中的第f帧混响序列，则各帧混响可表示为
R[f]＝βR0+(1-β)A[f] (10)
式中：β是一个用来控制混响变化的剧烈程度；R0表示慢变混响，由实际水域中散射体的分布与性质决定，主要由海底混响组成；A[f]表示第f帧快变混响，由海面波浪、卷入水体的气泡等引起，会对目标探测结果造成干扰。
若用Ai,j[f]表示A[f]的第i行、第j列的元素，则
Ai,j[f]＝αAi,j[f-1]+(1-α)v (11)
式中：0<α<1是一个随机变量，控制快变混响变化的剧烈程度。
在生成的混响数据中添加一个移动目标，目标强度表示为t，用矩阵Tf表示第f帧混响中添加的目标，Tf(m,n)表示Tf矩阵中第m行、第n列的元素，则
使用上述混响模型生成20帧数据。其中，Nr＝500，Nθ＝60，距离分辨率为1m，角度分辨率为1°。β＝0.8，α＝0.2。R0是一个服从二维瑞利分布的随机矩阵，其尺度参数设置为1。v服从期望为0、方差为1的高斯分布。目标强度为t＝0.8。BS-RPCA算法参数设置为：多帧累计数M＝5，迭代步长ρ＝1.5,，正则化参数，惩罚因子μ0＝1.25/||D||2，最大迭代次数kmax＝10，子块大小P×Q＝2×5。
使用该仿真结果对本发明所提的方法进行验证，实验结果如图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定义式为
式中，Es为所有旁瓣能量，Em为主瓣能量，二维情况下主瓣按照矩形来近似。图5为原始声呐图像、经APG、GoDec及本发明所提BS-RPCA方法处理后具体每帧所对应的二维积分旁瓣比(ISLR)，单位为dB，ISLR的数值越小，代表目标能量所占整幅图像能量比重越大，混响抑制效果更佳。通过图5我们可看出：本发明所提的BS-RPCA的混响抑制效果最佳，经本方法处理可降低原始声呐图像的ISLR约35dB，分别较APG、GoDec多降低约25dB、10dB。本发明提出的方法可显著降低声呐图像的积分旁瓣比，且相比现有的多帧混响抑制方法(APG、GoDec)效果更好，可有效提高目标的可检测性。
Claims (5)
1.一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，其特征在于，包括以下步骤：
步骤1：获得声纳图像观测矩阵，包括以下子步骤：
子步骤1.1：对主动声呐阵元域接收回波做匹配滤波与常规波束形成处理，幅值归一化得到声呐图像；
子步骤1.2：定义f表示当前帧数，f＝1,2,…，第f帧角度距离二维图像用表示，其中Nθ表示波束数，Nt时间采样点数；选取连续M帧结果构造二维数据观测矩阵；当f≥M时，取第f-M+1到第f帧的结果，记为张量/>则X＝(F[f-M+1],F[f-M+2],...,F[f])，其中F[·]表示括号内帧数的声呐图像；
子步骤1.3：对多帧结果X作向量化处理，将这些列向量按照帧号顺序组合成一个矩阵D∈RN×M，其中N＝Nθ×Nr，M为选取的帧数；
步骤2：设置分块正则化参数，包括以下步骤：
子步骤2.1：将大小为Nr×Nθ的声呐图像分为m×n个子块，子块大小为P×Q，P、Q为初始设定值；记WINi,j表示第i行、第j列子块的位置，Gi,j(p,q,m)表示在张量X中WINi,j内(p,q,m)位置的数值，其中p＝1,2,…,P，q＝1,2,…,Q，m＝1,2,…,M；
子步骤2.2：计算各子块特征值：
其中，表示遍历(p,q)对所有元素求和；
子步骤2.3：根据得到的各子块特征值得到各子块对应的正则化参数：
步骤3：建立块稀疏模型，包括以下子步骤：
子步骤3.1：建立块稀疏模型：
其中||·||F为矩阵的Frobenius范数，Pi,j表示将E中的每一列还原成矩阵的操作，返回WINi,j位置稀疏块在每一帧图像内的矩阵表示；
步骤4：运用非精确拉格朗日乘子法对该模型进行求解，得到稀疏矩阵E的最优解，包括以下子步骤：
步骤4.1：定义增广拉格朗日函数的为：
式中，Y为拉格朗日乘子，μ为一个正标量的惩罚参数；<Y,D-A-E>表示矩阵Y与矩阵D-A-E的内积；
步骤4.2：对增广拉格朗日函数中的A、E、Y进行初始化参数设置，E0＝0，A0＝0，Y0＝D/max(||D||2,λ0 -1||D||∞)，其中||·||2为矩阵的谱范数，称为行和范数；定义惩罚参数μ0，迭代步长ρ，最大迭代次数为kmax，设置迭代次数k＝0；
步骤4.3：固定Ek、Yk，更新低秩矩阵Ak+1
其中Sε[·]表示常规软阈值算子，
步骤4.4：固定Ak+1、Yk，更新稀疏矩阵Ek+1
表示块收缩算子，/>
其中，；
步骤4.5：更新Yk+1
Yk+1＝Yk+μk(D-Ak+1-Ek+1)
步骤4.6：更新μk+1
μk+1＝min(ρμk,107μ0)
步骤4.7：更新迭代次数k＝k+1；
当满足收敛条件||D-Ak-Ek||F/||D||F<10-7或达到最大迭代次数kmax时，输出此时的Ak+1、Ek+1，此时的Ak即为最终分解得到的低秩矩阵，属于被抑制掉的混响部分；Ek即为最终分解得到的稀疏矩阵，包含运动的目标部分，是我们需要的部分；
当不满足收敛条件时，继续下一次循环，即返回到步骤4.3；
步骤5：逆向量化，得到最终输出结果；
将步骤4最终得到的稀疏矩阵E∈RN×M中的每列重组为原始单帧数据维度，用F′表示逆向量化后的结果，可得到张量[F1′,F′2,...,F′M]，F′M即为当前f帧输出的最终结果；至此，我们完成了第f帧数据的混响抑制处理，通过低秩稀疏矩阵分解在获得运动目标信息的同时抑制了慢变的混响成分；
将F′M输出后返回到步骤1进行第f+1帧数据处理，继续步骤1～5。
2.如权利要求1所述的一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，其特征在于，所述子步骤1.2中，最初开始处理的帧数f需要大于等于累积的帧数M，M取连续的3～10帧。
3.如权利要求1所述的一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，其特征在于，所属子步骤2.1中，P×Q为预设的子块大小，且P、Q取值的范围分别为1～3、5～10。
4.如权利要求1所述的一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，其特征在于，所属步骤2.2中，在计算各子块特征值时，如果出现εi,j＝0，则步骤2.3中λi,j直接记为λ0，λ0的取值建议为其中N＝Nθ×Nr。
5.如权利要求1所述的一种基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法，其特征在于，所述步骤4.2中，惩罚参数μ0为正值，迭代步长ρ大于1；最大迭代次数kmax设置为10次。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110249031.0A CN113050098B (zh) | 2021-03-08 | 2021-03-08 | 基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110249031.0A CN113050098B (zh) | 2021-03-08 | 2021-03-08 | 基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113050098A CN113050098A (zh) | 2021-06-29 |
| CN113050098B true CN113050098B (zh) | 2024-04-16 |
Family
ID=76510400
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110249031.0A Active CN113050098B (zh) | 2021-03-08 | 2021-03-08 | 基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113050098B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117250604B (zh) * | 2023-11-17 | 2024-02-13 | 中国海洋大学 | 一种目标反射信号与浅海混响的分离方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2016136229A (ja) * | 2015-01-14 | 2016-07-28 | 本田技研工業株式会社 | 音声処理装置、音声処理方法、及び音声処理システム |
| CN110287819A (zh) * | 2019-06-05 | 2019-09-27 | 大连大学 | 动态背景下基于低秩及稀疏分解的动目标检测方法 |
-
2021
- 2021-03-08 CN CN202110249031.0A patent/CN113050098B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2016136229A (ja) * | 2015-01-14 | 2016-07-28 | 本田技研工業株式会社 | 音声処理装置、音声処理方法、及び音声処理システム |
| CN110287819A (zh) * | 2019-06-05 | 2019-09-27 | 大连大学 | 动态背景下基于低秩及稀疏分解的动目标检测方法 |
Non-Patent Citations (3)
| Title |
|---|
| 一种加权稀疏约束稳健Capon波束形成方法;刘振;孙超;刘雄厚;郭祺丽;物理学报;20161231(第010期);全文 * |
| 基于RPCA的小孔径垂直阵辐射噪声测量方法;蒋国庆;孙超;刘雄厚;蒋光禹;哈尔滨工程大学学报;20201231;第41卷(第010期);全文 * |
| 被动时反混响抑制方法研究;郭国强;杨益新;孙超;;自然科学进展;20090115(第01期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113050098A (zh) | 2021-06-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106248340B (zh) | 一种基于三维超声成像技术的风洞模型3d冰形在线测量方法 | |
| CN109116311B (zh) | 基于知识辅助稀疏迭代协方差估计的杂波抑制方法 | |
| CN104035095B (zh) | 基于空时最优处理器的低空风切变风速估计方法 | |
| CN107576961B (zh) | 一种互质降采样间歇合成孔径雷达稀疏成像方法 | |
| CN110780271A (zh) | 基于卷积神经网络的空间目标多模式雷达分类方法 | |
| CN106338713B (zh) | 一种基于波束零陷权的矢量阵目标左右舷分辨方法 | |
| CN108646247B (zh) | 基于伽马过程线性回归的逆合成孔径雷达成像方法 | |
| CN109541585B (zh) | 一种基于峰度评估的人体穿墙检测成像方法 | |
| CN101876715A (zh) | 一种拖曳声学阵列的拖船噪声抑制方法 | |
| CN108020817A (zh) | 基于配准的机载前视阵雷达杂波抑制方法 | |
| CN113589287B (zh) | 合成孔径雷达稀疏成像方法、装置、电子设备及存储介质 | |
| CN104656073A (zh) | 三维成像声纳波束形成方法及在多核处理器上的实现方法 | |
| CN113050098B (zh) | 基于块稀疏稳健主成分分析的反蛙人声呐混响抑制方法 | |
| CN108761417B (zh) | 基于知识辅助最大似然的机载雷达杂波抑制方法 | |
| Ye et al. | Mobilenetv3-yolov4-sonar: Object detection model based on lightweight network for forward-looking sonar image | |
| CN112215832B (zh) | Sar尾迹图像质量评估及自适应探测参数调整方法 | |
| CN110850421A (zh) | 基于混响对称谱的空时自适应处理的水下目标检测方法 | |
| CN116403100A (zh) | 一种基于矩阵分解的声呐图像小目标检测方法 | |
| CN113176573B (zh) | 复数域结构化sar舰船目标动态仿真与速度估计方法 | |
| CN113640793B (zh) | 基于mrf的实孔径扫描雷达超分辨成像方法 | |
| CN113406630B (zh) | 一种星载多通道sar动目标成像的误差估计和补偿算法 | |
| CN114720944A (zh) | 一种基于svd迭代熵值的墙体杂波抑制方法 | |
| CN111414580B (zh) | 一种低信混比条件下的混响抑制方法 | |
| CN111273299B (zh) | 一种针对组网声纳的水下分布式压制性干扰布置方法 | |
| CN113589300A (zh) | 一种基于压缩感知的合成孔径声纳沉底目标成像增强方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |