一种分量数时变的水声信号瞬时频率估计方法

1.本发明涉及水声信号处理技术领域，尤其是一种分量数时变的水声信号瞬时频率估计方法。


背景技术：

2.在海洋实时监测、异常声波事件检测等应用中，通常需要重构特定分量的信息。瞬时频率是时频变换的一个重要属性，通过时频谱和瞬时频率，能重构水声信号所有分量的信息。因此，估计时频谱中的瞬时频率，能精确提取信号中特殊的目标信息。
3.对于多分量的信号，最为直接的方式是经验模式分解或变模式分解，将非平稳信号分解成一系列固有模式函数，所有模式参数可以用于重构分量的瞬时频率。然而，对于具有紧邻甚至频点重叠特性的水声信号，这些方法难以提取出所有分量的瞬时频率。一种替代方法是时频分析方法，如短时傅里叶变换、调频变换和二次时频变换方法，提取信号的时频谱。在时频平面，用时频变换的拟变换重构出分量信息。由于这些方法受噪声干扰，提取的时频分辨率低，估计的瞬时频率误差较大。为了提高时频分辨率，采用重分配方法改进时频谱的可读性，但是重分配方法破坏了时频变换的可逆性。为了解决这个问题，同步压缩变换采用频域重分配方法，在增强时频谱能量集中度的同时，保留时频变换的可逆性。不幸的是，同步压缩变换逆变换受分离条件限制，不能重构出紧邻甚至重叠分量的瞬时频率。固有模式调频分解、固有调频分量分解、核稀疏学习方法能估计具有重叠分量信号的瞬时频率，但是必须设定好信号分量数。
4.由于水声信号具有快时变、低信噪比、紧邻甚至频点重叠特性，至今为止，对于这种水声信号中分量数随时间变化的情况，还没有采用同步压缩变换和自动脊线提取相结合进行瞬时频率估计。


技术实现要素：

5.本发明的目的在于提供一种能够处理分量数随时间变化的信号，克服噪声干扰影响水声信号瞬时频率估计的精度，准确地估计水声信号的瞬时频率的分量数时变的水声信号瞬时频率估计方法。
6.为实现上述目的，本发明采用了以下技术方案：一种分量数时变的水声信号瞬时频率估计方法，该方法包括下列顺序的步骤：
7.(1)输入多分量水声信号s(t)，利用扩展的同步压缩变换esst得到多分量水声信号s(t)的时频谱t(t，ω)；
8.(2)利用自动脊线提取area，追踪到时频谱t(t，ω)中所有可能的脊线，然后利用脊线融合标准，融合脊线；
9.(3)对得到的脊线用傅里叶级数展开，用tikhonov正则化对脊线的傅里叶系数重构，得到重构的瞬时频率。
10.在步骤(1)中，所述扩展的同步压缩变换esst的公式如下：
[0011][0012]
式中，为公共实数集，h(0)为0时刻的高斯窗函数值，δ(
·
)为delta函数，x
[n]
(t，f)为第n阶同步压缩变换sst算子，s
mstfft
(t，f)采用最大化短时间分数傅里叶变换的时频谱，定义为：
[0013][0014]
式中，σ(t，f)和c(t，f)分别是时频变的标准差和啁啾率，h
σ(t，f)，c(t，f)
(
·
)为具有时频变的标准差和啁啾率的高斯窗函数；用k个不同常量的标准差σk和啁啾率ck计算得到k个esst，求k个esst最值作为最大化的s
mstfft(t，f)
。
[0015]
所述步骤(1)具体包括以下步骤：
[0016]
(1a)用短时分数傅里叶变换计算初始的时频谱，估算初始的瞬时频率ω(t，f)，迭代n次，得到ω
[n]
(t，f)；
[0017]
(1b)迭代k次，计算具有标准差σk和啁啾率ck的窗函数和第k个esst的时频谱其中k∈[1，k]；
[0018]
(1c)从k个有限数据集中，筛选出时频谱最大值，得到扩展同步压缩变换的时频谱
[0019]
在步骤(2)中，所述自动脊线提取area具体包括以下步骤：
[0020]
(2a)依据时频谱的幅度系数μ＝min{e(ω)}/max{e(ω)}，其中为时频谱全局能量值，计算自动脊线提取终止的能量阈值δ
t
＝μmax{|f(t，ω)|}；
[0021]
(2b)时频谱能量最大值点作为脊线的起始点在时频谱中提取第k个脊线，对于τ＝t
max
+1，
…
，n
t
时频点进行遍历，选择出符合如下条件的时频点：
[0022]
δω∈[δω-，δω
+
]
[0023]
δω-＝max(1，ιm(τ-1)-δf)
[0024]
δω
+
＝min(nf，ιm(τ-1)+δf)
[0025][0026][0027]
t(t，δω)＝0
[0028]
其中，δω是频带范围，δω-和δω
+
是频带的下界和上界，ιm(τ)为脊线频率索引，为τ时刻脊线能量值，δf为采样的频宽；
[0029]
对于τ＝t
max-1，
…
，1时频点进行遍历，按照步骤(2b)所述的条件，筛选出脊线时频
点；
[0030]
(2c)更新第k个脊线对应的时频谱能为零，即t(t，δω)＝0；
[0031]
(2d)计算第k个脊线的三个属性值：脊线两端点、脊线面积、脊线偏向角θk；
[0032]
(2e)重复步骤(2b)至步骤(2d)，直至max{|t(t，ω)|}＞δ
t
条件不满足，终止迭代，至此，提取时频谱中所有可能的脊线；
[0033]
(2f)根据融合标准，对所有提取到的脊线进行融合。
[0034]
在步骤(2)中，所述融合标准具体是指：
[0035]
对于脊线lk，在所提取的脊线中需要符合时间先后标准，筛去不符合时间顺序的脊线，得到可融合的脊线集：
[0036][0037]
其中，t
ke
为脊线lk的时间结束点值，t
ie
为脊线li的时间结束点值；
[0038]
对符合时间标准的脊线集κk，用角度逼近标准，对脊线集进行滤波，
[0039][0040]
其中，θk是脊线lk的偏向角，θ
ik
是来自于脊线lk和脊线li上距离最近的两点组成的线的偏向角；
[0041]
判断两个脊线是否能融合，必须满足面积-间隙惩罚系数和可融合标准：
[0042][0043]
η
ik
＜＝λ，|t
ke-t
is
|＜δt，|f
ke-f
is
|＜2δf
[0044]
其中，ηk是lk的面积值，t
is
为脊线li的时间开始点值，nk为脊线lk上的采样点个数，f
ke
为脊线lk上频率结束点值，f
is
为脊线li上频率起始点值，δt为采样时宽。
[0045]
所述步骤(3)具体是指：
[0046]
对于第m个分量对应的脊线进行傅里叶级数展开为：
[0047][0048]
其中，l为频率的傅里叶阶数，1≤c≤l，和为傅里叶序列系数，fs为采样频率；
[0049]
用矩阵的形式表示为：
[0050]
其中，系数2l+1阶傅里叶级数表示为；
[0051][0052]
最后，用tikhonov正则化对傅里叶展开系数重构，得到重构的瞬时频率。
[0053]
由上述技术方案可知，本发明的有益效果为：第一，本发明是对自动脊线提取，与
传统重叠分量的瞬时频率估计相比，本发明能够处理分量数随时间变化的信号；第二，本发明结合分数时频谱，采用扩展同步压缩变换提取高能量集中度的时频谱，能够克服噪声干扰影响水声信号瞬时频率估计的精度；第三，本发明方法能够准确的估计水声信号的瞬时频率。
附图说明
[0054]
图1为本发明的方法流程图；
[0055]
图2(a)为采用本发明得到分量数时变的仿真信号的时频谱；
[0056]
图2(b)为采用本发明得到分量数时变的仿真信号的瞬时频率估计；
[0057]
图3为信噪比为[-5，10]db范围下仿真信号瞬时频率估计误差示意图；
[0058]
图4(a)为采用本发明得到鲸鱼信号的时频谱；
[0059]
图4(b)为采用本发明得到鲸鱼信号的瞬时频率；
[0060]
图5(a)为采用本发明得到海豚信号的时频谱；
[0061]
图5(b)为采用本发明得到海豚信号的瞬时频率。
具体实施方式
[0062]
如图1所示，一种分量数时变的水声信号瞬时频率估计方法，该方法包括下列顺序的步骤：
[0063]
(1)输入多分量水声信号s(t)，利用扩展的同步压缩变换esst(extended synchrosqueezing transform)得到多分量水声信号s(t)的时频谱t(t，ω)；esst的含义为扩展的同步压缩变换；
[0064]
(2)利用自动脊线提取area(automatic ridge extraction)，追踪到时频谱t(t，ω)中所有可能的脊线，然后利用脊线融合标准，融合脊线；area的含义为自动脊线提取；
[0065]
(3)对得到的脊线用傅里叶级数展开，用tikhonov正则化对脊线的傅里叶系数重构，得到重构的瞬时频率。
[0066]
在步骤(1)中，所述扩展的同步压缩变换esst的公式如下：
[0067][0068]
式中，为公共实数集，h(0)为0时刻的高斯窗函数值，δ(
·
)为delta函数，ω
[n]
(t，f)为第n阶同步压缩变换sst算子，s
mstfft
(t，f)采用最大化短时间分数傅里叶变换的时频谱，定义为：
[0069][0070]
式中，σ(t，f)和c(t，f)分别是时频变的标准差和啁啾率，h
σ(t，f)，c(t，f)
(
·
)为具有时频变的标准差和啁啾率的高斯窗函数；用k个不同常量的标准差σk和啁啾率ck计算得到k个esst，求k个esst最值作为最大化的s
mstfft
(t，f)。
[0071]
所述步骤(1)具体包括以下步骤：
[0072]
(1a)用短时分数傅里叶变换计算初始的时频谱，估算初始的瞬时频率ω(t，f)，迭
代n次，得到ω
[n]
(t，f)；
[0073]
(1b)迭代k次，计算具有标准差σk和啁啾率ck的窗函数和第k个esst的时频谱其中k∈[1，k]；
[0074]
(1c)从k个有限数据集中，筛选出时频谱最大值，得到扩展同步压缩变换的时频谱
[0075]
在步骤(2)中，所述自动脊线提取area具体包括以下步骤：
[0076]
(2a)依据时频谱的幅度系数μ＝min{e(ω)}/max{e(ω)}，其中为时频谱全局能量值，计算自动脊线提取终止的能量阈值δ
t
＝μmax{|t(t，ω)|}；
[0077]
(2b)时频谱能量最大值点作为脊线的起始点在时频谱中提取第k个脊线，对于τ＝t
max
+1，
…
，n
t
时频点进行遍历，选择出符合如下条件的时频点：
[0078]
δω∈[δω-，δω
+
]
[0079]
δω-＝max(1，ιm(τ-1)-δf)
[0080]
δω
+
＝min(nf，ιm(τ-1)+δf)
[0081][0082][0083]
t(t，δω)＝0
[0084]
其中，δω是频带范围，δω-和δω
+
是频带的下界和上界，ιm(τ)为脊线频率索引，为τ时刻脊线能量值，δf为采样的频宽；
[0085]
对于τ＝t
max-1，
…
，1时频点进行遍历，按照步骤(2b)所述的条件，筛选出脊线时频点；
[0086]
(2c)更新第k个脊线对应的时频谱能为零，即t(t，δω)＝0；
[0087]
(2d)计算第k个脊线的三个属性值：脊线两端点、脊线面积、脊线偏向角θk；
[0088]
(2e)重复步骤(2b)至步骤(2d)，直至max{|t(t，ω)|}＞δ
t
条件不满足，终止迭代，至此，提取时频谱中所有可能的脊线；
[0089]
(2f)根据融合标准，对所有提取到的脊线进行融合。
[0090]
在步骤(2)中，所述融合标准具体是指：
[0091]
对于脊线lk，在所提取的脊线中需要符合时间先后标准，筛去不符合时间顺序的脊线，得到可融合的脊线集：
[0092][0093]
其中，t
ke
为脊线lk的时间结束点值，t
ie
为脊线li的时间结束点值；
[0094]
对符合时间标准的脊线集κk，用角度逼近标准，对脊线集进行滤波，
[0095]
[0096]
其中，θk是脊线lk的偏向角，θ
ik
是来自于脊线lk和脊线li上距离最近的两点组成的线的偏向角；
[0097]
判断两个脊线是否能融合，必须满足面积-间隙惩罚系数和可融合标准：
[0098][0099]
η
ik
＜＝λ，|t
ke-t
is
|＜δt，|f
ke-f
is
|＜2δf
[0100]
其中，ηk是lk的面积值，t
is
为脊线li的时间开始点值，nk为脊线lk上的采样点个数，f
ke
为脊线lk上频率结束点值，f
is
为脊线li上频率起始点值，δt为采样时宽。
[0101]
所述步骤(3)具体是指：
[0102]
对于第m个分量对应的脊线进行傅里叶级数展开为：
[0103][0104]
其中，l为频率的傅里叶阶数，1≤c≤l，和为傅里叶序列系数，fs为采样频率；
[0105]
用矩阵的形式表示为：
[0106]
其中，系数2l+1阶傅里叶级数表示为；
[0107][0108]
最后，用tikhonov正则化对傅里叶展开系数重构，得到重构的瞬时频率。
[0109]
图2(a)中，时长为1秒的仿真信号加了高斯白噪声，信噪比为-3db，采样率为3khz，窗函数宽的参数为0.05。用自动脊线提取area，追踪到时频谱中所有分量的瞬时频率，如图2(b)所示。从图2(b)可以看出，本发明能估计分量数时变信号的瞬时频率。仿真信号在不同信噪比[-5，10]db条件下，计算估计的瞬时频率和真实的瞬时频率的均方误差。
[0110]
如图3所示，本发明与固有模式调频分解(intrinsic mode chirp decomposition，imcd)、固有调频分量分解(intrinsic chirp component decomposition，iccd)、核稀疏学习(kernel sparse learning，ksl)和快瞬时频率(fast instantaneous frequency，fast-if)进行估计结果对比，结果表明，本发明对应的曲线area在不同信噪比下能估计具有分量数时变信号的瞬时频率。
[0111]
本发明的扩展同步压缩变换结合多个短时分数傅里叶变换，最大化地增强沿脊线的时频能量集中度，获取高分辨率的时频谱，并降低噪声干扰，有利于降低瞬时频率估计的均方误差。
[0112]
图4(a)中，采样点为20002，采样频率为22.05khz，用自动脊线提取和脊线重构估计时频谱中所有分量数时变的瞬时频率，如图4(b)所示。图5(a)为实际海豚信号的时频谱，采样点为50001，采样频率为44.1khz。图5(b)为最终估计到的瞬时频率。
[0113]
综上所述，本发明是对自动脊线提取，与传统重叠分量的瞬时频率估计相比，本发明能够处理分量数随时间变化的信号；本发明结合分数时频谱，采用扩展同步压缩变换提取高能量集中度的时频谱，能够克服噪声干扰影响水声信号瞬时频率估计的精度；本发明能够准确的估计水声信号的瞬时频率。