一种基于组稀疏结构的水声目标辐射噪声调制谱重构方法与流程

本发明属于声纳信号处理领域，具体涉及一种基于组稀疏结构的水声目标辐射噪声调制谱重构方法。

背景技术

舰船辐射噪声是良好的水下声源。水下目标的机械振动噪声和螺旋桨噪声是水面舰船、潜艇、鱼雷等水声目标的主要噪声源,螺旋桨空化噪声会产生幅度调制,具有鲜明的节奏感,通过解调处理后在频谱图上表现为线状特征的谱线,称之为“线谱”。调制线谱在频谱图上的位置较为稳定,能够反映诸如螺旋桨转速、叶片数目等一些目标特性,是舰船辐射噪声的不变特征,在一定信噪比下由于它的线状特征也较容易发现、分离。这个特性使之成为目标检测和识别的重要因素之一。因此辐射噪声解调谱分析是舰船噪声自动目标识别的重要分析手段之一。

对于舰船辐射噪声来说,包络调制使辐射噪声往往呈现明显的节奏感，调制频率中包含螺旋桨的轴频、叶频以及它们的谐波频率分量等信息，在宽带解调中通过滤波和检波可以获得其解调谱。对于这种信号的解调通常用绝对值低通解调,平方低通解调和希尔伯特幅值解调方法,所谓绝对值低通解调就是对宽带噪声先做取绝对值非线性运算,再通过低通滤波器获得调制信号成分,而平方低通解调是对宽带噪声先做平方运算,再通过低通滤波器获得调制信号成分。希尔伯特幅值解调方法是让信号频谱的幅度不发生改变,引起频谱变化的只是其相位,从而降低信号的采样率。

但是，螺旋桨节拍是由舰船的某些周期性振动的调幅效应引起的辐射噪声响度周期性变化形成的，它是对全频域水声目标辐射噪声进行调制，其中螺旋桨节拍对其辐射的空化噪声有明显的振幅调制作用，其调制频率及调制深度与螺旋桨转速、桨叶数及舰船航速等有关。舰船辐射噪声中所表现出的常见节拍类型主要包括：叶片频率、轴频、轴频加叶片频、轻重节奏(往复)、轻重节奏加叶片频率、内燃机汽缸频率加轻重节奏等。往往在不同的频带呈现不同的节奏类型，且各频带上的调制深度和信噪比是不同的，这是因为在不同频带中引起不同节奏的噪声源的辐射强度不同。因此，舰船辐射噪声的调制方式并不是单一的，在频域上也不均匀，它是宽带中各种节奏的混合。因此，常规宽频带解调法往往很难获得清晰的调制谱。对这类信号通常采用窄带解调技术。通常的做法是将信号通过带通滤波器，得到不同频段的子带信号，然后进行解调。

传统的解调谱估计一般以fft为基础,运算速度快,实现简单。但fft方法得到的信号线谱分辨率受窗长影响，就可能使比较重要的谱线不能被检测出来，导致在水声测量中信号的重要频域信息被遗漏，另一方面，传统基于多子带分频技术调制谱分析方法各子频带处理均是相互独立的，未利用子频带间子带调制谱的关联性，因此需要寻找新的解决方法进行调制谱分析，得到高精度的调制谱。

技术实现要素：

发明目的：本发明旨在提供一种水声目标辐射噪声调制谱重构方法，该方法可以自动学习出调制谱的稀疏度，同时充分利用了水声目标辐射噪声子带间调制谱位置的关联性，实现了水声目标辐射噪声调制谱的高分辨率重构。

技术方案：本发明采用如下技术方案：

一种基于组稀疏结构的水声目标辐射噪声调制谱重构方法，包括如下步骤：

(1)模拟水声目标辐射噪声中的连续谱分量rc(t)和线谱分量rl(t)，构成水声目标辐射噪声r(t)，r(t)＝rc(t)+rl(t)；

(2)对水声目标辐射噪声r(t)进行幅度调制，得到水声目标辐射噪声调制信号x(t)；

(3)对水声目标辐射噪声调制信号x(t)归一化，利用带通滤波器获取l个不相邻子频带的噪声调制信号yl(t)，其中l∈[1,...,l]；

(4)估计各子带上噪声调制信号yl(t)的幅度调制数据其中l＝1,2,...,l；

(5)对估计出的幅度调制数据进行离散采样并用稀疏的频率系数表示；通过设计稀疏的频率系数(l＝1,2,...,l)的先验分布，构建基于组稀疏结构的高分辨调制谱生成模型；

(6)基于期望最大化方法推导稀疏的频率系数的后验分布；

(7)利用参数估计公式，迭代求解估计出高分辨稀疏调制谱。

步骤(1)中所述水声目标辐射噪声中的连续谱分量rc(t)的获取步骤如下：

(a.1)采用三参数模型法模拟平稳连续谱的功率谱gxf(ωt)：

其中ωm，ωc和λ为三参数模型的三个参数，决定了该连续谱的形状；ωt为频率，ωm为尖锐度因子，决定谱锋的尖锐程度和高度，ωc决定谱锋的位置，λ决定功率谱高、低频端幅度的相对比例，σ表示平稳连续谱信号的能量；

(a.2)建立p阶ar滤波器，其yule-walker方程为：

其中a[q]，q∈{1,2,…,p}和p0为p阶ar滤波器系数，δ[k]为冲击函数；rx[k]为gxf(ωt)的自相关函数rc(τ)的采样值；

(a.3)采用levison-durbin算法求解式(2)方程，得到p阶ar滤波器系数；高斯白噪声通过该ar滤波器后得到的信号，即为水声目标辐射噪声中的平稳连续谱分量rc(t)。

步骤(1)中所述水声目标辐射噪声中的线谱分量rl(t)的获取步骤如下：

(b.1)采用k个正弦信号来模拟目标信号的线谱分量，其中ak为正弦信号幅度，fk为正弦信号的频率，t∈[0,t]为观测时间；

(b.2)在线谱位置fk处计算平稳连续谱分量rc(t)的能量pik，k＝1,2,...,k；

(b.3)根据已知的信号干扰比计算出各个正弦信号幅度ak，即得到水声目标辐射噪声中的线谱分量rl(t)。

步骤(2)中对水声目标辐射噪声r(t)进行幅度调制，得到水声目标辐射噪声调制信号x(t)，x(t)可表示为：

x(t)＝a[1+αs(t)]r(t)

其中a为信号的幅值，为无量纲比例因子。s(t)为调制函数。

步骤(3)包括如下步骤：

(3.1)对水声目标辐射噪声调制信号x(t)进行归一化处理；归一化后的水声目标辐射噪声调制信号y(t)可表示为：

y(t)＝[1+αs(t)]r(t)

(3.2)利用带通滤波器将归一化后全频带的水声目标辐射噪声调制信号y(t)分解为l个子频带的信号yl(t)，其中l∈[1,...,l]。

步骤(4)包括如下步骤：

(4.1)通过绝对值检波器获取水声目标幅度调制信息；

(4.2)通过低通滤波以及直流抑制估计出幅度调制数据考虑载波为单频信号的解调情况，假设r(t)＝cosωt，ω为载频，经过绝对值检波器的信号可表示为：

经过低通滤波以及直流抑制后，可估计出水声目标辐射噪声信号的调制频率成分：

其中，低通滤波的截止频率flpf满足ω<flpf<2ω-ω。

步骤(5)包括如下步骤：

(5.1)将估计出的幅度调制数据进行离散采样后得到采样后的长度为n，可稀疏表示为：

式(3)中εl为加性噪声，服从均值为0、方差为β0in的复高斯分布：

其中in为n×n维单位矩阵，β0为加性噪声方差，服从超参数为c0＝d0＝0的gamma分布：

式(4)中为不同频率点向量构成的复感知矩阵，为复数标志，n×m表示该复感知矩阵的维数，且为：

其中，表示第i个频率点，频率间隔m为重构的稀疏向量维数。假设水声目标辐射噪声信号的观测时间为t,采样频率为fs，则观测向量的点数n＝fs×t，重构的稀疏向量维数m满足m＞＞n；

(5.2)设计稀疏的频率系数(l＝1,2,...,l)的先验分布；服从均值为0、方差为diag(γ)的复高斯分布：

式(5)中为复随机变量x服从均值为a，方差为b的复高斯分布，diag(γ)表示向量γ的元素构成的对角矩阵，γ^-1满足gamma分布：

其中，γ(a0)为gamma函数，不失一般性，超参数a0＝b0＝0。

步骤(6)包括如下步骤：

(6.1)基于期望最大化方法推导出频率系数的后验分布：

p为的概率分布函数，式中

γ＝diag(γ)(10)

(6.2)通过最大化边缘似然函数的方法估计出参数γ和β0：

式中，γi是参数向量γ的第i个参数，σii为矩阵σ第i行第i列上的元素，μil是第l个均值向量的第i个元素，是前一次迭代过程中的β0值，c为n维矩阵c∈r^n×n,可以表示为：

c＝β0in+φγφ^h(13)

步骤(7)包括如下步骤：

(7.1)迭代计算式(8)、式(9)、式(11)和式(12)直至收敛，计算出参数γ和β0；

(7.2)通过式(8)、式(9)以及估计出的参数γ和β0，得到稀疏系数的均值μl和方差σ，即稀疏系数服从均值为μl，方差为σ的复高斯分布；

(7.3)对l组估计的具有组稀疏结构的系数进行能量累积获得高分辨率调制谱。

有益效果：本发明公开的基于组稀疏结构的水声目标辐射噪声调制谱重构方法，首先模拟出调制后的水声目标辐射噪声信号；然后对辐射噪声信号进行频带分解，在每个子频带上分别估计出调制函数；其次设计稀疏的频率系数先验分布，同时构建基于组稀疏结构的高分辨调制谱生成模型；接着基于期望最大化的方法推导稀疏频率系数的后验分布；最后，通过稀疏系数的先验分布和后验分布，迭代求解估计出对应子带的高分辨调制谱。与现有技术相比，本发明公开的方法具有以下优点：该方法克服了传统的调制谱获取方法中频率分辨率不够的问题，实现高分辨率的调制谱重构；同时该方法在非参数化的稀疏贝叶斯框架下，可以自动学习出调制谱的稀疏度，又充分利用了水声目标辐射噪声子带间调制谱位置的关联性，使调制谱的重构有较好的质量和效果。

附图说明

图1为本发明方法的流程图；

图2为实施例1中目标辐射噪声信号的频谱图；

图3为实施例1中目标辐射噪声信号的调制谱真值示意图；

图4为实施例1中在第一个子频带上采用传统傅里叶变换法估计出的调制频谱图；

图5为实施例1中在第二个子频带上采用传统傅里叶变换法估计出的调制频谱图；

图6为实施例1中在第三个子频带上采用传统傅里叶变换法估计出的调制频谱图；

图7为实施例1中在第四个子频带上采用传统傅里叶变换法估计出的调制频谱图；

图8为实施例1中基于组稀疏结构重构得到的调制谱。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明的具体实施案例做说明。

一种基于组稀疏结构的水声目标辐射噪声调制谱重构方法，如图1所示，包括如下步骤：

步骤1、模拟水声目标辐射噪声中的连续谱分量rc(t)和线谱分量rl(t)，构成水声目标辐射噪声r(t)，r(t)＝rc(t)+rl(t)；

水声目标辐射噪声中的连续谱分量rc(t)的获取步骤如下：

(a.1)采用三参数模型法模拟平稳连续谱的功率谱gxf(ωt)：

其中ωm，ωc和λ为三参数模型的三个参数，决定了该连续谱的形状；ωt为频率，ωm为尖锐度因子，决定谱锋的尖锐程度和高度，ωc决定谱锋的位置，λ决定功率谱高、低频端幅度的相对比例，σ表示平稳连续谱信号的能量；

(a.2)根据wiener-khinchin定理，式(1)的逆傅里叶变换即为平稳连续谱信号的自相关函数rc(τ)，可写成：

rc(τ)＝σexp(-ωm|τ|)[cosωcτ+λsin(ωc|τ|)]

假设以fs为采样率对时域信号进行等间隔采样，则上述自相关函数可以写成离散形式为：

rc(kts)＝σexp(-ωm|kts|)[cosωckts+λsin(ωc|kts|)]

其中ts＝1/fs；根据式(1)建立p阶ar滤波器，其yule-walker方程为：

其中a[q]，q∈{1,2,…,p}和p0为p阶ar滤波器系数，δ[k]为冲击函数；rx[k]为gxf(ωt)的自相关函数rc(τ)的采样值；

(a.3)采用levison-durbin算法求解式(2)方程，得到p阶ar滤波器系数；高斯白噪声通过该ar滤波器后得到的信号，即为水声目标辐射噪声中的平稳连续谱分量rc(t)。

水声目标辐射噪声中的线谱分量rl(t)的获取步骤如下：

(b.1)采用k个正弦信号来模拟目标信号的线谱分量，其中ak为正弦信号幅度，fk为正弦信号的频率，t∈[0,t]为观测时间；

(b.2)在线谱位置fk处计算平稳连续谱分量rc(t)的能量pik，k＝1,2,...,k；

(b.3)根据已知的信号干扰比计算出各个正弦信号幅度ak，即得到水声目标辐射噪声中的线谱分量rl(t)。

步骤2、对水声目标辐射噪声r(t)进行幅度调制，得到水声目标辐射噪声调制信号x(t)；x(t)可表示为：

x(t)＝a[1+αs(t)]r(t)

其中a为信号的幅值，为无量纲比例因子。s(t)为调制函数，不失一般性，可以假设为s(t)＝spsin(ωt)，即若s(t)为周期函数，可以分解为各次谐波之和，其中sp为调制深度，ω为调制频率。

步骤3、对水声目标辐射噪声调制信号x(t)归一化，利用带通滤波器获取l个不相邻子频带的噪声调制信号yl(t)，其中l∈[1,...,l]；包括如下步骤：

(3.1)对水声目标辐射噪声调制信号x(t)进行归一化处理；归一化后的水声目标辐射噪声调制信号y(t)可表示为：

y(t)＝[1+αs(t)]r(t)

(3.2)利用带通滤波器将归一化后全频带的水声目标辐射噪声调制信号y(t)分解为l个子频带的信号yl(t)，其中l∈[1,...,l]。

步骤4、估计各子带上噪声调制信号yl(t)的幅度调制数据其中l＝1,2,...,l；包括如下步骤：

(4.1)通过绝对值检波器获取水声目标幅度调制信息；

(4.2)通过低通滤波以及直流抑制估计出幅度调制数据

考虑载波为单频信号的解调情况，假设r(t)＝cosωt，ω为载频，经过绝对值检波器的信号可表示为：

经过低通滤波以及直流抑制后，可估计出水声目标辐射噪声信号的调制频率成分：

其中，低通滤波的截止频率flpf满足ω<flpf<2ω-ω。

步骤5、对估计出的幅度调制数据进行离散采样并用稀疏的频率系数表示；通过设计稀疏的频率系数(l＝1,2,...,l)的先验分布，构建基于组稀疏结构的高分辨调制谱生成模型；包括如下步骤：

(5.1)将估计出的幅度调制数据进行离散采样后得到采样后的长度为n，可稀疏表示为：

式(3)中εl为加性噪声，服从均值为0、方差为β0in的复高斯分布：

其中in为n×n维单位矩阵，β0为加性噪声方差，服从超参数为c0＝d0＝0的gamma分布：

式(4)中为不同频率点向量构成的复感知矩阵，为复数标志，n×m表示该复感知矩阵的维数，且为：

其中，表示第i个频率点，频率间隔m为重构的稀疏向量维数。假设水声目标辐射噪声信号的观测时间为t,采样频率为fs，则观测向量的点数n＝fs×t，重构的稀疏向量维数m满足m＞＞n；

(5.2)设计稀疏的频率系数(l＝1,2,...,l)的先验分布；服从均值为0、方差为diag(γ)的复高斯分布：

式(5)中为复随机变量x服从均值为a，方差为b的复高斯分布，diag(γ)表示向量γ的元素构成的对角矩阵，γ^-1满足gamma分布：

其中，γ(a0)为gamma函数，不失一般性，超参数a0＝b0＝0。

步骤6、基于期望最大化方法推导稀疏的频率系数的后验分布；包括如下步骤：

(6.1)基于期望最大化方法推导出频率系数的后验分布：

p为的概率分布函数，式中

γ＝diag(γ)(10)

(6.2)通过最大化边缘似然函数的方法估计出参数γ和β0：

式中，γi是参数向量γ的第i个参数，σii为矩阵σ第i行第i列上的元素，μil是第l个均值向量的第i个元素，是前一次迭代过程中的β0值，c为n维矩阵c∈r^n×n,可以表示为：

c＝β0in+φγφ^h(13)

步骤7、利用参数估计公式，迭代求解估计出高分辨稀疏调制谱；包括如下步骤：

(7.1)迭代计算式(8)、式(9)、式(11)和式(12)直至收敛，计算出参数γ和β0；

(7.2)通过式(8)、式(9)以及估计出的参数γ和β0，得到稀疏系数的均值μl和方差σ，即稀疏系数服从均值为μl，方差为σ的复高斯分布；

(7.3)对l组估计的具有组稀疏结构的系数进行能量累积获得高分辨率调制谱。

实施例1：

本实施例中，采样频率fs＝4khz。利用三参数模型法模拟水声目标辐射噪声的平稳连续谱的功率谱gxf，仿真过程中三参数设置如下：尖锐度因子ωm＝200hz，谱峰中心位置因子ωc＝500hz，谱高、低频段幅度的相对比例影响因子λ＝0，平稳连续谱信号能量σ＝500。

模拟目标辐射噪声的3个线谱分量：设置正弦信号的频率fi分别为20hz，50hz，100hz及相应的4次谐波线谱。观测时间为t＝10s。将平稳连续谱分量以及线谱分量加起来，获得目标辐射噪声信号r(t)。目标辐射噪声信号的频谱如图2所示。

本实施例中，调制谱s(t)由轴频、叶频以及对应的谐波分量组成，其中轴频设置为0.8hz线谱，且具有四次谐波分量。螺旋桨为7叶，其叶频为5.6hz，具有3次谐波分量。轴频及其谐波调制深度服从均值为0.1，方差为0.01的高斯分布；叶频深度同样服从高斯分布，其均值为轴频调制深度的两倍，其叶频谐波调制深度随谐波次数呈指数衰减。具体调制谱如图3所示。通过阵列波束形成信号增强后的加性噪声信号信噪比为10db。

在此实施例中，把舰船辐射噪声信号分解为4个子频带分别解调，4个子频带的频率范围分别为0hz～100hz，100hz～500hz，500hz～1000hz，1000hz～2000hz。图4至图7分别给出了基于传统傅里叶变换的方法得出的各个子带的调制谱估计结果。从图中可以看出，在某些子频带中轴频线谱和叶频线谱频率的调制幅度较低，用传统的门限检测方法容易被遗漏。

图8给出了用本发明的方法得到的调制谱估计结果。从图中可以看出，与传统方法相比，本发明公开的方法能够有效地检测出包含轴频和叶频及其谐波分量的所有调制谱。