CN109085595A - 一种使用水听器接收信号估计空中运动声源速度的方法 - Google Patents

Abstract

本发明涉及一种使用水听器接收信号估计空中运动声源速度的方法，将PCT与Ferguson和Lo提出的使用单水听器接收信号估计空中目标速度的方法结合在一起。首先，利用PCT估计水听器接收信号在整个时间序列上的瞬时频率，再从瞬时频率中计算4个参数即声源频率f₀、运动速度v、飞行高度h和CPA时间μ_c的初始估计值，最后利用Ferguson和Lo建立的瞬时频率参数模型在最小二乘法准则下迭代拟合瞬时频率估计值求解目标速度。有益效果：估计出了较为准确的信号瞬时频率。并以此为基础利用瞬时频率和声源参数的关系实现了使用水听器接收信号对空中运动声源速度的精确估计。上述方法比利用STFT得到瞬时频率的参数估计方法更为准确。

Description

一种使用水听器接收信号估计空中运动声源速度的方法

技术领域

本发明属于信号处理、时频分析理论和水声测量等领域涉及一种使用水听器接收信号估计空中运动声源速度的方法，利用时频分析方法以及水听器接收信号瞬时频率 与空中声源参数之间的关系实现使用水听器接收信号对空中运动声源的速度估计。

背景技术

估计空中运动声源的速度是水下对空预警技术的一个重要方面。很大一部分飞机是由螺旋桨驱动的，其辐射噪声主要分为宽带连续谱噪声和线谱噪声两部分，线谱成 分呈谐波关系，主要集中在低频段，强度相比于宽带连续谱噪声较高，跨空气-水界面 远距离传播后更容易被探测。飞机在巡航时基本可以认为是保持匀速直线运动的，当 水面上方的飞机相对于水听器运动时，线谱噪声会发生多普勒频移现象。水听器接收 信号的瞬时频率中包含着声源频率、运动速度、飞行高度以及声源经过与水听器的最 近会遇点(ClosestPoint of Approach,CPA,如图1所示)时间4个空中声源参数的信 息。因此，利用水听器接收到声源辐射噪声的瞬时频率就可以对空中声源的速度进行 估计。

Ferguson和Lo利用线谱辐射噪声的多普勒效应，建立了瞬时频率与速度等参数的模型，提出了一种使用单个水听器的接收信号估计空中目标速度的方法，并通过实验 评估了该方法的性能。这种方法首先利用短时傅里叶变换(Short-Time Fourier Transform,STFT)从水听器接收信号中估计出其瞬时频率，然后利用瞬时频率估计空 中声源的飞行速度(B.G.Ferguson,K.W.Lo.Transiting aircraft parameter estimation usingunderwater acoustic sensor data[J].IEEE Journal of Oceanic Engineering,1999,24(4): 424-435.)。然而，在水下环境噪声的干扰下，STFT的时频分辨率和能量聚集度较低， 估计瞬时频率的精度下降从而影响上述速度估计方法的性能。

为了提高瞬时频率的估计精度，需要选取一种能量聚集度高且对噪声有抑制作用的时频分析方法，以准确地估计水听器接收信号的瞬时频率。Z.K.Peng和G.Meng 等人在传统Chirplet变换的基础上提出了多项式Chirplet变换(Polynomial Chirplet Transform,PCT)，PCT具有能量聚集度高的优点，可以有效地抑制噪声的干扰，并且 适用于非线性调频信号(Z.K.Peng,G.Meng,et al.Polynomial chirplet transform with application toinstantaneous frequency estimation[J].IEEE Journal of Oceanic Engineering,2011,60(9):3222-3229.)。而空中运动声源在水下产生的多普勒频移信号正是一种典型的非线性调频信号，所以PCT非常适合用来估计的水听器接收信号的瞬时频率。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种使用水听器接收信号估计空中运动声源速度的方法。

技术方案

一种使用水听器接收信号估计空中运动声源速度的方法，其特征在于步骤如下：

步骤1：放大水听器接收到的信号并记录，记为s_r(t)，声源为亚音速飞机，并且 声波在跨介质入水的过程中入射角小于临界角12.8°；声源某一特定线谱的原始频率 恒为f₀，运动速度为v，高度为h，CPA时间为μ_c；

步骤2：对记录到的信号s_r(t)利用STFT进行预处理，得到接收信号的时频分布 后采用带通滤波，对滤波后的信号进行希尔伯特变换得到相应的复信号z(t)；

步骤3：利用PCT对预处理后的信号z(t)进行时频分析，得到预处理后信号z(t)的时频分布：

PCT的定义式

其中：τ表示对信号做PCT的时刻，f表示信号频率，w_α(t)表示非负、对称且归 一化实窗；

采用高斯函数，α给定了高斯窗的时间跨度，和分别为频率旋转算子和频移算子；

所述β＝(β₁,β₁,,β_N)为多项式系数向量；

通过迭代实现以β为系数的多项式对信号瞬时频率的最优逼近，收敛后得到最优多项式系数，相应的PCT结果将在下一步被用来估计信号z(t)的瞬时频率；

步骤4：在时刻τ_k对信号能量在频率上的分布进行三点二次插值，能量峰值对应的频率即为信号的瞬时频率估计值以此获得信号在采样时间上的瞬时频率估计序 列其中，τ_k表示时间序列中的第k个时刻，K表示时间序列的长度；

步骤5：以水下观测空中运动线谱声源的瞬时频率模型在非线性最小二乘准则下拟合逼近上一步得到的瞬时频率估计序列通过迭代算法估计声源参数， 当算法收敛后即得到声源速度的估计值

有益效果

本发明提出的一种使用水听器接收信号估计空中运动声源速度的方法，将PCT与Ferguson和Lo提出的使用单水听器接收信号估计空中目标速度的方法结合在一起。首 先，利用PCT估计水听器接收信号在整个时间序列上的瞬时频率，再从瞬时频率中计 算4个参数即声源频率f₀、运动速度v、飞行高度h和CPA时间μ_c的初始估计值，最 后利用Ferguson和Lo建立的瞬时频率参数模型在最小二乘法准则下迭代拟合瞬时频 率估计值求解目标速度。

本发明的有益效果是：利用适用于非线性调频信号的PCT对多普勒频移信号进行处理，估计出了较为准确的信号瞬时频率。并以此为基础利用瞬时频率和声源参数的 关系实现了使用水听器接收信号对空中运动声源速度的精确估计。上述方法比利用 STFT得到瞬时频率的参数估计方法更为准确。

附图说明

图1是声音从声源跨空气-水介质到达水听器的传播过程示意图；

图2是使用水听器接收信号估计空中运动声源速度的总体流程图；

图3是PCT中多项式系数优化算法的流程图。

图4是利用PCT得到的信号z(t)的时频分布

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明解决其技术问题所采用的技术方案总体流程如图2所示，具体包括以下步骤：

1)放大水听器接收到的信号并记录，记为s_r(t)。

2)对记录到的信号s_r(t)进行预处理，得到z(t)。

3)利用PCT得到预处理后信号z(t)的时频分布。

4)对信号的时频分布做平滑处理后提取瞬时频率。

5)利用非线性最小二乘法得到声源速度的估计值

以下对本发明的每一个步骤做详细说明：

所述步骤1)的具体实现如下：

声音从声源跨空气-水介质到达水听器的传播过程如图1所示，图中给出了空中声源的飞行轨迹和水听器的空间位置。水听器布放深度为d，目标在巡航状态下以恒定 的高度、方向和速度从水听器上方飞过，经过水听器正上方(CPA)。在这个过程中， 水听器持续接收空中声源辐射的声信号，并将其转换成电信号，电信号通过功率放大 器后用数据采集仪器记录下来，记为s_r(t)。其中0＜t≤T，T表示数据记录的时间长 度。本发明适用的前提是声源为亚音速飞机，并且声波在跨介质入水的过程中入射角 小于临界角(12.8°)。记声源某一特定线谱的原始频率恒为f₀，运动速度为v，高度 为h，CPA时间为μ_c。

所述步骤2)的具体实现如下：

首先利用STFT粗略得到接收信号的时频分布，从中可以大致确定出声源辐射噪声某频率线谱多普勒频移的范围，采用带通滤波器滤除该频率范围以外的干扰，对滤 波后的信号进行希尔伯特变换得到相应的复信号z(t)。

所述步骤3)的具体实现如下：

利用PCT对预处理后的信号z(t)进行时频分析。PCT的定义为

其中τ表示对信号做PCT的时刻，f表示信号频率，w_α(t)表示非负、对称且归一化实窗，通常采用高斯函数，有

w_α(t)＝exp[(logα)t²] (2)

α给定了高斯窗的时间跨度，β＝(β₁,β₁,,β_N)为多项式系数向量，N为多项式阶数，和分别为频率旋转算子和频移算子，表达式为

为了得到最优的多项式系数β，通过迭代用多项式拟合逼近信号的瞬时频率，其流程如图3所示：第一步设定多项式系数初值为β₀＝(0,0,...,0)、并给定多项式阶数N和 判决门限δ；第二步将多项式系数β代入式(1)中得到z(t)关于(τ,f)的时频分布，从中 提取峰值所对应的瞬时频率随τ的变化轨迹IF(τ)；第三步用多项式在最小 二乘准则下拟合逼近前一步得到的瞬时频率，获得多项式系数β_i；第四步计算

其中τ_k表示时间序列中的第k个时刻，K表 示时间序列的长度，i表示当前迭代的次数，i-1表示上一次迭代的次数，若ε≤δ则结 束迭代，否则重复第二步至第三步直到ε≤δ；迭代结束后得到的多项式系数β即为多 项式系数。

所述步骤4)的具体实现如下：

将最优多项式系数β带入到式(1)中得到信号z(t)的时频分布，为了提高频率分 辨率、降低噪声的干扰，对其做平滑处理后提取瞬时频率。在给定时刻τ_k对信号能量 在频率上的分布PCT(τ_k,f)做关于频率的三点二次插值，得到更加光滑的分布 PCT′(τ_k,f)。能量峰值对应的频率为信号的瞬时频率估计值，即 重复这一操作，以此获得信号在采样时间上的瞬时频率估计 序列

所述步骤5)的具体实现如下：

用水下观测空中运动线谱声源的瞬时频率模型在非线性最小二乘准则下拟合逼近 上一步得到的瞬时频率估计序列，估计声源参数，即求解以下最小化问题：

其中，为τ_k时刻水听器瞬时频率的估计值，f′₀,v′,h′,μ′_c为待定的声源参数，

f_k(f′₀,v′,h′,μ′_c)为τ_k时刻水听器接收信号瞬时频率的理论值，满足以下模型：

其中

γ＝c_a/c_w，c_a和c_w分别为空气中和水中的声速，x(μ)＝v′(μ-μ_c)，l_w(μ)为方程x(μ)²+d²(1+γh′/l_w(μ))²＝(l_w(μ)+γh′)²的正实根。

公式(4)给出的是一个非线性最小二乘优化问题，只能通过迭代求解，为了提高迭代的速度，首先选取几个特殊时刻，粗略估计4个声源参数作为迭代的初始值声源频率的初始估计值为

飞行速度的初始估计值为

其中，通过计算瞬时频率序列前M个值的均值和最后M个值的均值分别可以得到和即

飞行高度的估计值为

其中，是τ_k时刻瞬时频率估计值的变化率，k为某一时刻在时间序列中排列的位 置。可以利用有限差分由f(τ)近似得到。最小时τ_k的取值，即

CPA时间的初始估计值为

由于噪声的干扰，瞬时频率估计序列上存在扰动，因此仅利用几个时刻的瞬时频率估计值得到的参数初始估计值的精度不高。为了充分利用完整的瞬时频 率序列信息，把参数初始估计值代入式(4)中，利用通用的非线性最小二乘 方法(如Matlab中的lsqnonlin函数采用的Levenberg-Marquardt算法)求解该问题， 即可得到声源速度的精确估计值

具体实施例子

为了验证利用本发明给出的方法估计空中运动声源参数的有效性，设计仿真实验如下：假设空气中声速c_a＝340m/s，水中声速c_w＝1524m/s，水听器到海平面的垂直距 离d＝50m，空中声源的线谱辐射噪声频率f₀＝68Hz，做速度v＝150m/s的匀速直线运 动、飞行高度h＝300m、CPA时间μ_c＝0s，采样率f_s＝512Hz，观察时间t＝[-1,3]s， 在这段时间内，声波在跨空气-水界面传播时的入射角小于临界角。PCT中多项式阶数N取为20，判决门限δ取为0.001，高斯窗长度取为1024，α取为0.005，计算时 M取为25。由于空中声源是运动状态，它与水听器之间的距离是变化的，所以水听器 接收信号的幅值也随之改变，以致信噪比也是不断改变的。因此假设水听器接收到的 背景噪声大小不变，考虑空中声源辐射噪声的强度按球面扩展规律衰减，在观察时间 两端处信噪比SNR＝0dB，噪声为高斯白噪声。根据以上参数得到仿真的水听器接收信 号，并使用本发明给出的方法对仿真信号进行处理，得到水听器接收信号的时频分布 如图4，空中声源参数估计值如表1。

表1空中运动声源参数估计值

从表1中可以看出利用本发明给出的方法可以较为准确地估计出空中运动声源的速度等参数，证明了这种方法的有效性。

Claims (1)

1.一种使用水听器接收信号估计空中运动声源速度的方法，其特征在于步骤如下：

步骤1：放大水听器接收到的信号并记录，记为s_r(t)，声源为亚音速飞机，并且声波在跨介质入水的过程中入射角小于临界角12.8°；声源某一特定线谱的原始频率恒为f₀，运动速度为v，高度为h，CPA时间为μ_c；

步骤2：对记录到的信号s_r(t)利用STFT进行预处理，得到接收信号的时频分布后采用带通滤波，对滤波后的信号进行希尔伯特变换得到相应的复信号z(t)；

步骤3：利用PCT对预处理后的信号z(t)进行时频分析，得到预处理后信号z(t)的时频分布：

PCT的定义式

其中：τ表示对信号做PCT的时刻，f表示信号频率，w_α(t)表示非负、对称且归一化实窗；

采用高斯函数，α给定了高斯窗的时间跨度，和分别为频率旋转算子和频移算子；

所述β＝(β₁,β₁,...,β_N)为多项式系数向量；

通过迭代实现以β为系数的多项式对信号瞬时频率的最优逼近，收敛后得到最优多项式系数，相应的PCT结果将在下一步被用来估计信号z(t)的瞬时频率；

步骤4：在时刻τ_k对信号能量在频率上的分布进行三点二次插值，能量峰值对应的频率即为信号的瞬时频率估计值以此获得信号在采样时间上的瞬时频率估计序列其中，τ_k表示时间序列中的第k个时刻，K表示时间序列的长度；

步骤5：以水下观测空中运动线谱声源的瞬时频率模型在非线性最小二乘准则下拟合逼近上一步得到的瞬时频率估计序列通过迭代算法估计声源参数，当算法收敛后即得到声源速度的估计值