CN108646248A - 一种针对低速运动声源的被动声学测速测距方法 - Google Patents

Abstract

本发明涉及一种针对低速运动声源的被动声学测速测距方法，根据传感器的接收信号s(n)的STFT分析结果确定待分析谱线所在的频率范围，该范围应当将这根谱线由于多普勒效应所造成的最大频率偏移和最小频率偏移包含在内，输入待定参数，通过定义D‑C‑CZT方法得到s(n)的时频分布，从时频分布中提取瞬时频率变化，基于最小二乘准则得到待定参数估计值，将待定参数估计值作为下次迭代的待定参数输入值，重复计算直至得到准确的速度估计值。本发明可以精确的提取出传感器接收信号的瞬时频率变化，尤其适用于传感器接收信号的多普勒频移特征较为微弱时的瞬时频率变化分析，能够获得较高精度的声源运动速度等参数的估计结果。

Description

一种针对低速运动声源的被动声学测速测距方法

技术领域

本发明属于被动声学测速测距方法，涉及一种针对低速运动声源的被动声学测速测距方法，特别是涉及一种匀低速直线运动条件下利用单个声学传感器对运动噪声源的速度及其与传感器之间的距离进行窄带被动估计的算法。

背景技术

经过检索，发现两个与本申请相关的专利。专利(CN106054159A)“一种多普勒信号的瞬时频率提取方法”发明了一种可以对各类激光速度干涉仪的输出信号进行瞬时频率提取的算法，算法中首先使用直接正交法提取离散信号瞬时频率，然后采用多项式拟合离散时刻的瞬时频率，在使用最小二乘准则确定多项式系数后，各个时刻瞬时频率的估计值都可以根据多项式得到。本申请在发明目的、技术方法和所利用的目标物理特征上均有别于专利(CN106054159A)。专利(CN1588124A)“利用多普勒变换的被动测速测距法及装置”发明了一种应用于水下和空中动目标检测和定位的一种利用多普勒变换的被动测速测距法和装置，该发明利用基于匀速运动模型的Dopplerlet变换，分析空中(包括水下和陆上)运动目标辐射出的声音信号，通过定义Dopplerlet基函数集，采用自适应匹配投影塔形分解法，搜索出一组与信号的主要成分最佳匹配的Dopplerlet基函数，进而准确地计算出目标的距离、速度以及声波传播速度等参数。本申请在发明目的和所利用的目标物理特征上与专利(CN1588124A)基本一致，在技术方法上与之相异。

汽车、飞机、船只等一类运动目标都可以看成是辐射噪声源，当其相对于接收系统运动时，接收系统输出信号的频率会因多普勒效应产生相应的变化，目标的速度、距离等信息就蕴含在瞬时频率变化之中。Reid等人在研究中指出估计速度等参数的关键在于精确估计观测声信号的瞬时频率(David C.Reid,Abdelhak M.Zoubir,and BoualemBoashash.Aircraft flight parameter estimation based on passive acoustictechniques using the polynomial Wigner–Ville distribution[J].The Journal ofthe Acoustical Society of America,1997,102(1):207-223.)。徐灵基等人使用PCT(Polynomial Chirplet Transform)、WVD(Wigner-Ville Distribution)和STFT(ShortTime Fourier Transform)提取瞬时频率变化曲线，成功地从一段匀速直线运动卡车的实测数据中获得了卡车的速度、距离等参数，研究中发现利用PCT方法较利用WVD方法和STFT方法可获得更精确的速度估计结果(Xu L,Yang Y,Yu S.Analysis of moving sourcecharacteristics using polynomial chirplet transform[J].Journal of theAcoustical Society of America,2015,137(4):EL320-EL326.)。但是基于PCT的测速测距方法存在计算过程耗时较长的问题。此外，当目标以较低速度运动时，接收信号的多普勒频移较为微弱，PCT方法不能精确的提取出接收信号的瞬时频率(牺牲计算效率增加时间窗长度可以在一定程度上提高PCT的频率分辨能力，考虑到PCT的计算效率本来就低，牺牲计算效率的做法不可取)。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种针对低速运动声源的被动声学测速测距方法，针对低速运动声源的接收信号提供一种具有较高计算效率的高精度瞬时频率估计方法，以便于快速准确地获取目标的速度估计和距离估计。

技术方案

一种针对低速运动声源的被动声学测速测距方法，其特征在于步骤如下：

步骤1、获取传感器的接收信号s(n)：目标在匀低速直线运动过程中自身噪声信号被传感器捕获，传感器将声信号转换为电信号记为s(n)；

步骤2、根据STFT确定待分析谱线所在的频率范围：对s(n)做短时傅里叶变换STFT得到时频分布，从中选择一条能量最强的瞬时频率变化曲线，确定其最大频率偏移和最小频率偏移的频率范围；

步骤3、设定待定参数的初始值：

频移模型

式中n＝0,1,…,N-1指示离散时间变化，N表示离散信号s(n)的长度，声速c和采样间隔Δt为已知常量，待定的参数共有4个，分别是原始线谱的频率f₀、噪声源的速度v、CPA时刻τ_c和CPA距离R_c；

初始待定参数中频率f₀ ⁰设定为步骤2得到的频率范围的中心频率；

噪声源的初始速度v⁰、CPA时刻τ_c ⁰和CPA距离R_c ⁰设定为大于0的任意值；

步骤4、输入待定参数，通过D-C-CZT方法得到s(n)的时频分布：

D-C-CZT定义：

所述

式中k＝0,1,…,N_w-1指示离散频率变化，N_w代表加窗的长度，D＝(f₀,v,τ_c,R_c)代表迭代过程中不断更新的4个参数，是旋转算子，是频移算子，两个算子起提高时频聚集度的作用，w_σ是窗函数，A＝exp(j2πf_b)，W＝exp(j2πΔf)，其中f_b代表频率细化的起始频率，Δf代表频率细化的间隔；

所述频率细化的范围即为步骤2)中待分析谱线所在的频率范围；

高斯窗函数w_σ(n)＝exp(-0.5(n/σ)²)

式中，σ＝(N_w-1)/5表示高斯窗的标准差，|n|＞(N_w-1)/2时w_σ(n)＝0；

步骤5、从时频分布中提取瞬时频率变化：从步骤4得到的s(n)时频分布中提取谱峰得到一组不同时刻下的谱峰频率f(n)，即为瞬时频率变化；

步骤6、基于最小二乘准则得到待定参数估计值：基于最小二乘准则，利用频移模型拟合瞬时频率变化f(n)获取频移模型中待定参数的估计值：频率f₀ ^s、噪声源速度v^s、CPA时刻τ_c ^s和CPA距离R_c ^s，其中上角标s代表D-C-CZT迭代次数；

步骤7：将待定参数估计值作为下次迭代的待定参数输入值，重复步骤4～步骤6，当满足终止迭代的判定条件时，即得到准确的速度估计值；

所述终止迭代的判定条件：

所述步骤5中选取n在N_w/2到N-N_w/2-1之间的瞬时频率变化。

所述阈值δ为

有益效果

本发明提出的一种针对低速运动声源的被动声学测速测距方法，根据传感器的接收信号s(n)的STFT分析结果确定待分析谱线所在的频率范围，该范围应当将这根谱线由于多普勒效应所造成的最大频率偏移和最小频率偏移包含在内，输入待定参数，通过定义D-C-CZT方法得到s(n)的时频分布，从时频分布中提取瞬时频率变化，基于最小二乘准则得到待定参数估计值，将待定参数估计值作为下次迭代的待定参数输入值，重复计算直至得到准确的速度估计值。

本发明可以精确的提取出传感器接收信号的瞬时频率变化，尤其适用于传感器接收信号的多普勒频移特征较为微弱时的瞬时频率变化分析，能够获得较高精度的声源运动速度等参数的估计结果。

使用本发明的方法成功地从一段卡车的记录噪声中估计出了其行驶的速度为6.15m/s，CPA距离为37.89m，与真实的速度和CPA距离吻合。

附图说明

图1：匀低速直线运动的示意图

图2：典型的噪声功率谱示意图

图3：D-C-CZT方法迭代流程图

具体实施方式

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

如图1所示，目标以恒定速度v循直线经过接收传感器，该过程中运动声源与接收传感器最接近的点称为CPA(Closest Point of Approach)，定义CPA时刻指声源经过CPA的时刻，定义CPA距离指CPA与接收传感器之间的距离。根据多普勒效应，目标在接近传感器并远离传感器的过程中，传感器接收信号的瞬时频率会产生一个非线性的由高到低的变化。只要能够准确地提取出接收信号的瞬时频率变化信息，就可以得到蕴含在其中的速度、距离等信息。本发明使用一种基于CZT(Chirp Z-transform)的联合时频分析方法D-C-CZT(Doppler Chirplet CZT)提取瞬时频率变化曲线，实现声源速度、CPA距离的估计。

1)获取传感器的接收信号s(n)

目标在匀低速直线运动过程中自身噪声信号被传感器捕获，传感器将声信号转换为电信号，经过放大滤波处理后采样并记录，记为s(n)。

2)根据STFT确定待分析谱线所在的频率范围：对s(n)做短时傅里叶变换STFT得到时频分布，从中选择一条能量最强的瞬时频率变化曲线，确定其最大频率偏移和最小频率偏移的频率范围；这个频率范围就是CZT的分析频段。

如图2所示，目标噪声源典型的噪声功率谱一般由低频离散线谱(500Hz以下)和宽带连续谱构成，其中低频线谱由于能量远高于宽带连续谱因而更易于检测。测速测距之前根据信号s(n)的STFT分析结果确定待分析谱线所在的频率范围，该范围应当将这根谱线由于多普勒效应所造成的最大频率偏移和最小频率偏移包含在内，并且将相邻谱线排除在外。

3)设定待定参数的初始值：

根据假设的匀速直线运动条件以及相关的物理机理可以建立一个频移模型，用以描述声源匀速直线运动时传感器接收信号的瞬时频率变化。除去时间变量外，该模型还与有5个常量参数有关，这些参数中声速c一般已知，其余为待定的参数，分别是原始线谱的频率f₀、噪声源的速度v、CPA时刻τ_c和CPA距离R_c。一般情况下，初始待定参数中频率f₀ ⁰可设定为步骤2)中频率范围的中心频率，速度v⁰、CPA时刻τ_c ⁰和CPA距离R_c ⁰可根据先验知识设定为大于0的某个合适值，如无先验知识，也可设定为大于0的任意值。

定义f_D如下

上式为频移模型(离散形式)，描述了声源匀速直线运动时传感器接收信号的瞬时频率变化。式中n＝0,1,…,N-1指示离散时间变化，N表示离散信号s(n)的长度，声速c和采样间隔Δt为已知常量，待定的参数共有4个，分别是原始线谱的频率f₀、噪声源的速度v、CPA时刻τ_c和CPA距离R_c。一般情况下，初始待定参数中频率f₀ ⁰可设定为步骤2)中频率范围的中心频率，速度v⁰、CPA时刻τ_c ⁰和CPA距离R_c ⁰可根据先验知识设定为大于0的某个合适值，如无先验知识，也可设定为大于0的任意值。

如果声速未知，也可以将声速作为待定参数。此时声速的初值c⁰可以根据常识设置，一般水中声速为1500m/s，空气中声速为340m/s。

4)输入待定参数，通过D-C-CZT方法得到s(n)的时频分布

D-C-CZT定义如下

其中

式中k＝0,1,…,N_w-1指示离散频率变化，N_w代表加窗的长度，D＝(f₀,v,τ_c,R_c)代表迭代过程中不断更新的4个参数，是旋转算子，是频移算子，两个算子起提高时频聚集度的作用，w_σ是窗函数，A＝exp(j2πf_b)，W＝exp(j2πΔf)，其中f_b代表频率细化的起始频率，Δf代表频率细化的间隔，频率细化的范围即为步骤2)中待分析谱线所在的频率范围。

窗函数一般可以选择为高斯窗，定义如下

w_σ(n)＝exp(-0.5(n/σ)²)

式中，σ＝(N_w-1)/5表示高斯窗的标准差，|n|＞(N_w-1)/2时w_σ(n)＝0。

根据D-C-CZT的定义式即可得到s(n)的时频分布。计算中初始待定参数由步骤3)给定，后续迭代中待定参数由上次迭代的估计结果给出，D-C-CZT的迭代流程如图3所示。

5)从时频分布中提取瞬时频率变化

从步骤4)得到的s(n)时频分布中提取谱峰得到一组不同时刻下的谱峰频率f(n)。考虑到信号的两端存在截断效应，所以实际只选取n在N_w/2到N-N_w/2-1之间的瞬时频率变化。

6)基于最小二乘准则得到待定参数估计值

基于最小二乘准则，利用频移模型拟合瞬时频率变化f(n)的问题可写作如下形式

该优化问题可以采用高斯牛顿迭代方法解决，迭代公式如下

上式中s代表迭代次数， 代表雅可比矩阵。当高斯牛顿迭代算法收敛后可得到待定参数估计值。

7)将待定参数估计值作为下次迭代的待定参数输入值，重复步骤4)～6)直至得到准确的速度估计值

把步骤6)解出的运动噪声源的参数估计值作为下次迭代的待定参数输入值，重复步骤4)～6)，直至D-C-CZT迭代收敛，即可得到准确的速度估计值。终止迭代的判定条件可设定如下

式中，δ为一阈值。

实例的数据来源于International Student Challenge Problem in AcousticSignal Processing 2014。数据为一辆卡车以恒定速度通过一条直道时路边的麦克风所记录的一段时长为30秒的卡车辐射噪声，记录时空气中声速为347m/s，采样率为12000Hz。分析中首先降采样至1200Hz，使用STFT分析发现能量最强的一条线谱大约在118Hz左右，确定分析范围为118±5Hz，窗长为256点，设定初始待定参数为f₀ ⁰＝118,v⁰＝10,R_c ⁰＝10,τ_c ⁰＝15，迭代3次后输出结果如表1所示。参考值为官方在挑战结束后公布的各个参数的真值估计结果，估计值为本发明的参数估计结果。

表1 D-C-CZT参数估计结果

Claims (3)

1.一种针对低速运动声源的被动声学测速测距方法，其特征在于步骤如下：

步骤1、获取传感器的接收信号s(n)：目标在匀低速直线运动过程中自身噪声信号被传感器捕获，传感器将声信号转换为电信号记为s(n)；

步骤2、根据STFT确定待分析谱线所在的频率范围：对s(n)做短时傅里叶变换STFT得到时频分布，从中选择一条能量最强的瞬时频率变化曲线，确定其最大频率偏移和最小频率偏移的频率范围；

步骤3、设定待定参数的初始值：

频移模型

式中n＝0,1,…,N-1指示离散时间变化，N表示离散信号s(n)的长度，声速c和采样间隔Δt为已知常量，待定的参数共有4个，分别是原始线谱的频率f₀、噪声源的速度v、CPA时刻τ_c和CPA距离R_c；

初始待定参数中频率f₀ ⁰设定为步骤2得到的频率范围的中心频率；

噪声源的初始速度v⁰、CPA时刻τ_c ⁰和CPA距离R_c ⁰设定为大于0的任意值；

步骤4、输入待定参数，通过D-C-CZT方法得到s(n)的时频分布：

D-C-CZT定义：

所述

式中k＝0,1,…,N_w-1指示离散频率变化，N_w代表加窗的长度，D＝(f₀,v,τ_c,R_c)代表迭代过程中不断更新的4个参数，是旋转算子，是频移算子，两个算子起提高时频聚集度的作用，w_σ是窗函数，A＝exp(j2πf_b)，W＝exp(j2πΔf)，其中f_b代表频率细化的起始频率，Δf代表频率细化的间隔；

所述频率细化的范围即为步骤2)中待分析谱线所在的频率范围；

高斯窗函数w_σ(n)＝exp(-0.5(n/σ)²)

式中，σ＝(N_w-1)/5表示高斯窗的标准差，|n|＞(N_w-1)/2时w_σ(n)＝0；

步骤5、从时频分布中提取瞬时频率变化：从步骤4得到的s(n)时频分布中提取谱峰得到一组不同时刻下的谱峰频率f(n)，即为瞬时频率变化；

步骤6、基于最小二乘准则得到待定参数估计值：基于最小二乘准则，利用频移模型拟合瞬时频率变化f(n)获取频移模型中待定参数的估计值：频率f₀ ^s、噪声源速度v^s、CPA时刻τ_c ^s和CPA距离R_c ^s，其中上角标s代表D-C-CZT迭代次数；

步骤7：将待定参数估计值作为下次迭代的待定参数输入值，重复步骤4～步骤6，当满足终止迭代的判定条件时，即得到准确的速度估计值；

所述终止迭代的判定条件：

2.根据权利要求1所述针对低速运动声源的被动声学测速测距方法，其特征在于：所述步骤5中选取n在N_w/2到N-N_w/2-1之间的瞬时频率变化。

3.根据权利要求1所述针对低速运动声源的被动声学测速测距方法，其特征在于：所述阈值δ为