CN103344961A - 舰船速度和距离联合测量的被动声多普勒相位方法 - Google Patents

Abstract

本发明提供了一种舰船速度和距离联合测量的被动声多普勒相位方法，测量被测舰船匀速直线接近测量水听器并离开测量水听器水听器的过程中接收舰船辐射的声信号，并做快速傅立叶变换分析频谱，找到信号中存在的信噪比大于或等于5dB的低频线谱，通过带通滤波器分离出各个低频线谱多普勒信号，复数化得到多普勒解析信号；获得多普勒信号的时变幅度估计值；通过多线谱合成的加权非线性最小二乘估计法得到舰船航行速度的估计值、正横距离的估计值、各线谱声源经过测量水听器时的正横时刻估计值、各线谱声源的频率估计值和多普勒信号的初始相位估计值。本发明实施简单，测量精度较高，计算量相对较小。

Description

舰船速度和距离联合测量的被动声多普勒相位方法

技术领域

本发明涉及一种联合测量舰船航行速度和正横距离的被动声多普勒相位方法，涉及信号处理、水声探测等领域。

背景技术

速度和距离是水下运动目标的关键参数，所以速度和距离的测量对于水声探测具有重要意义，同时也是声呐的一项重要任务。水下测速测距方法包括有源测速测距和无源测速测距，这两类方法在技术上有着本质上的不同，因而测速测距的精度也不大相同。然而，无论何种测速测距方法都是利用速度和距离不同引起的信号的各种变化来进行间接测量。与有源测速测距相比，无源测速测距具有隐蔽性，其也称为被动测速测距。而利用单点被动测速测距实施简易，有着广泛的需求，可应用于声呐浮标测速测距、水声试验中目标的运动分析和水下目标声源级测量等。

单点被动测速测距利用的是水声信道或水下目标的运动特征。杨娟（杨娟，惠俊英，江磊，等.基于STFT-Hough变换的目标运动分析[J].哈尔滨工程大学学报,2007,28(2):137-142.）根据舰船辐射的连续谱信号在水声信道中多途传播形成的时频域干涉谱，提出一种联合估计水下运动目标速度和正横距离的方法。该方法依赖于具体的水声场环境，在浅海条件下效果才明显，其适用范围受限。基于目标与接收点的相对运动引起的多普勒效应，Ferguson（Ferguson B G.A ground-based narrow-band passiveacoustic technique for estimating the altitude and speed of a propeller-driven aircraft[J].Journal of the Acoustical Society of America,1992,92(3):1403-1407.）、Reid（Reid D C,Zoubir A M,Boashash B.Aircraft flight parameter estimation based on passive acoustictechniques using the polynomial Wigner-Ville distribution[J].Journal of the AcousticalSociety of America,1997,102(1):207-223.）、李宏（李宏，王昭，赵俊渭.进化谱估计及其在飞机定位中的应用[J].西北工业大学学报,1998,16(4):511-515.）发展了一类利用接收到的线谱声信号的多普勒频移来估计旋翼式飞机飞行参数（速度、高度和距离）的方法。该类方法也适用于水下运动目标测速测距，但由于其需要通过瞬时频率估计算法提取接收信号的多普勒频移曲线，计算量较大。吴国清（吴国清，陈永强.线谱非平稳性与Dopplerlet测距测速[J].声学学报,2003,28(2):130-136.）提出利用Dopperlet算法估计舰船航行速度和正横距离的方法，但Dopperlet算法参数搜索过程繁琐，且速度和正横距离的估计精度较低。吴国清（吴国清，马力.利用辐射噪声多线谱的多普勒进行距离估计[J].声学学报,2006,31(2):140-145.）后又提出了改进，根据多线谱的多普勒频移估计舰船正横距离，该方法利用了更多的多普勒频移信息，虽然提高了正横距离估计的精度，但本质上还是使用Dopplerlet算法。严光洪（YAN Guang-hong,CHEN Zhi-fei,and SUN Jin-cai.Using a linear array to estimate the velocity of underwatermoving targets[J].Journal of Marine Science and Application,2009,8(4):343-347.）利用阵列的高分辨测向技术来测量水下运动目标的速度，理论上得到了较精确的结果，然而该方法对阵列布放的形状要求严格，较难应用于实际的测量且也只能测量水下运动目标的航行速度。

发明内容

为了克服现有技术测量精度较低且计算量较大的不足，本发明提供一种舰船速度和距离联合测量的被动声多普勒相位方法，利用舰船与测量水听器相对运动形成的多普勒信号的相位特性，在仅使用单个水听器的简易测量条件下，可以实现对舰船的航行速度和正横距离的联合精确测量。

本发明解决其技术问题所采用的技术方案包括以下步骤：

1）测量被测舰船匀速直线接近测量水听器并离开测量水听器水听器的过程中接收舰船辐射的声信号，并将其转换成电压信号s(t),0＜t≤T，通过前放后用数据采集仪记录下来，t表示记录时间，T表示数据记录的时间长度；

2）把采集得到的水听器接收信号s(t)做快速傅立叶变换分析频谱，找到信号中存在的信噪比大于或等于5dB的低频线谱，其频率为

m＝1,2,…,M，M表示低频线谱的个数；通过带通滤波器分离出存在线谱的频带宽度为

各个小区域信号，得到各个低频线谱多普勒信号s_m(t)；将分离出的线谱多普勒信号s_m(t)复数化得到多普勒解析信号z_m(t)，其中z_m(t)＝s_m(t)+jH[s_m(t)]，H[s_m(t)]表示信号s_m(t)的希尔伯特变换；

3）多普勒信号的表达式为z_m(t)＝A_m(t)exp[jφ_m(t)]，其中，时变幅度 $A_m\left(t\right)=\frac{A_{m0}(c^2-v^2)}{c[\sqrt{{c^2{v}^2(t-t_{m0})}^2+R_0^2(c^2-v^2)}-v^2(t-t_{m0})]},$ 瞬时相位 ${\mathrm{\φ}}_m\left(t\right)=2\mathrm{\π}{f}_{m0}\frac{c^2}{c^2-v^2}[t-\frac{v^2}{c^2}{t}_{m0}-\frac{1}{c^2}\sqrt{v^2{c}^2{(t-t_{m0})}^2+R_0^2(c^2-v^2)}]+{\mathrm{\φ}}_{m0},$ A_m0表示舰船各辐射线谱声信号的幅度，f_m0表示各辐射线谱声信号的频率，t_m0为各线谱声源通过正横位置离测量水听器最近时的正横时刻，φ_m0为多普勒信号z_m(t)的初始相位；将多普勒信号z_m(t)表示为阶数为I的高阶多项式相位信号

其中a_mi,i＝0,1,…,I为多普勒信号z_m(t)的多项式相位系数；利用高阶模糊函数算法估计出多项式相位系数

i＝0,1,…,I,进而得到多普勒信号z_m(t)的瞬时相位估计值

将接收到的多普勒信号z_m(t)乘上

对乘积的实部进行低通滤波，获得多普勒信号z_m(t)的时变幅度估计值

4）根据步骤(3)中多普勒信号z_m(t)的瞬时相位表达式，并利用获得的瞬时相位估计值

和时变幅度估计值

通过多线谱合成的加权非线性最小二乘估计法得到舰船航行速度的估计值、正横距离的估计值、各线谱声源经过测量水听器时的正横时刻估计值、各线谱声源的频率估计值和多普勒信号的初始相位估计值： $\{\hat{v},{\hat{R}}_0,{\hat{f}}_{m0},{\hat{t}}_{m0},{\widehat{\mathrm{\φ}}}_{m0}\}=\underset{(v,R_0,f_{m0},t_{m0},{\mathrm{\φ}}_{m0})}{\arg \min }\left\{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\∫}_0^T{\hat{A}}_m^2{[{\mathrm{\φ}}_m\left(t\right)-{\widehat{\mathrm{\φ}}}_m\left(t\right)]}^2\mathrm{dt}\right\}.$

本发明的有益效果是：在现有的单个水听器测量水下运动目标参数的平台上，根据测量时舰船与测量水听器之间相对运动产生的多普勒效应，利用高阶模糊函数方法，得到多普勒信号的瞬时相位和时变幅度，然后使用多线谱合成的加权非线性最小二乘法估计出多普勒信号的参数，实现了舰船航行速度和正横距离的联合测量，该方法实施简单，测量精度较高。而且，由于可以采用Levenberg-Marquardt迭代算法实现多线谱合成的加权非线性最小二乘法估计舰船的运动参数，计算量相对较小。

附图说明

图1是舰船运动模型示意图；

图2是舰船航行速度和正横距离联合测量方法的总体流程框图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明包括以下步骤：

1）舰船航行速度测量时，水中声速c可提前测得，视为已知量。被测舰船以速度ν做匀速直线航行，由远及近接近测量水听器并由近及远离开测量水听器，测量水听器与舰船运动轨迹的正横距离为R₀，在这个过程中测量水听器接收舰船辐射的声信号，并将其转换成电压信号s(t),0＜t≤T，通过前放后用数据采集仪记录下来，T表示数据记录的时间长度。

2）由于被测舰船与测量水听器之间存在相对运动，会引起多普勒效应，测量水听器接收到的线谱声信号称为多普勒信号。测量水听器记录的信号需要预处理，目的为：一是确定存在的强线谱声信号；二是滤波去除多普勒信号的带外噪声；三是复数化得到解析形式的多普勒信号。先把采集得到的水听器接收信号s(t)做快速傅立叶变换（FFT）分析频谱，找到信号中存在的信噪比（线谱信号的功率与其所在频率点的噪声功率的比值）大于或等于5dB的低频线谱，其频率一般选择为m＝1,2,…,M，M表示低频线谱的个数。通过带通滤波器分离出存在线谱的频带宽度为

各个小区域信号，这便得到各个低频线谱多普勒信号s_m(t),m＝1,2,…,M，t表示时间。将分离出的线谱多普勒信号s_m(t)复数化得到多普勒解析信号z_m(t)，其中z_m(t)＝s_m(t)+jH[s_m(t)]，H[s_m(t)]表示信号s_m(t)的希尔伯特变换。

3）多普勒信号的表达式为

其中，时变幅度 $A_m\left(t\right)=\frac{A_{m0}(c^2-v^2)}{c[\sqrt{{c^2{v}^2(t-t_{m0})}^2+R_0^2(c^2-v^2)}-v^2(t-t_{m0})]},$ 瞬时相位 ${\mathrm{\φ}}_m\left(t\right)=2\mathrm{\π}{f}_{m0}\frac{c^2}{c^2-v^2}[t-\frac{v^2}{c^2}{t}_{m0}-\frac{1}{c^2}\sqrt{v^2{c}^2{(t-t_{m0})}^2+R_0^2(c^2-v^2)}]+{\mathrm{\φ}}_{m0},$ A_m0表示舰船各辐射线谱声信号的幅度，f_m0表示各辐射线谱声信号的频率，t_m0为各线谱声源通过正横位置离测量水听器最近时的正横时刻，φ_m0为多普勒信号z_m(t)的初始相位。根据Stone-Weierstrass理论，多普勒信号z_m(t)可以表示为阶数为I的高阶多项式相位信号：

其中a_mi,i＝0,1,…,I为多普勒信号z_m(t)的多项式相位系数。利用高阶模糊函数（High-order Ambiguity Function，HAF）算法估计出多项式相位系数

i＝0,1,…,I,进而得到多普勒信号z_m(t)的瞬时相位估计值

由于多普勒信号的幅度是慢变的，将接收到的多普勒信号z_m(t)乘上

对乘积的实部进行低通滤波，可以获得多普勒信号z_m(t)的时变幅度估计值

4）根据步骤(3)中多普勒信号z_m(t)的瞬时相位表达式，并利用上一步骤已获得的瞬时相位估计值

和时变幅度估计值

通过多线谱合成的加权非线性最小二乘估计法得到舰船航行速度的估计值、正横距离的估计值、各线谱声源经过测量水听器时的正横时刻估计值、各线谱声源的频率估计值和多普勒信号的初始相位估计值： $\{\hat{v},{\hat{R}}_0,{\hat{f}}_{m0},{\hat{t}}_{m0},{\widehat{\mathrm{\φ}}}_{m0}\}=\underset{(v,R_0,f_{m0},t_{m0},{\mathrm{\φ}}_{m0})}{\arg \min }\left\{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\∫}_0^T{\hat{A}}_m^2{[{\mathrm{\φ}}_m\left(t\right)-{\widehat{\mathrm{\φ}}}_m\left(t\right)]}^2\mathrm{dt}\right\}.$

图1为舰船运动模型示意图，测量水听器位于S处，被测舰船匀速航行，航行速度为ν，轨迹为MN，正横位置O为其运行轨迹上离测量水听器距离最近点，正横距离为R₀。被测舰船通过测量水听器的过程中，测量水听器接收舰船辐射的声信号，并将其转换成电压信号s(t)，t为时间，通过前置放大后用数据采集仪以采样率F_S记录时间长度为T的数据，得到信号s(t)的离散时间序列s(n),n＝1,…,N，N＝TF_S，n与t的对应关系为t＝nΔ，Δ＝1F_S。

图2为舰船航行速度和正横距离联合测量方法的总体流程框图，具体实施如下：

1）对水听器接收信号预处理，先将采集到的信号s(n)做快速傅立叶变换（FFT）分析其频谱，确定信号中存在的信噪比（线谱信号的功率与其所在频率点的噪声功率的比值）大于或等于5dB的低频线谱，其频率一般选择为

m＝1,2,…,M，M表示低频线谱的个数。通过FIR带通滤波器分离出存在线谱且频带宽度为

的各个小区域信号，c为水声传播速度，这便得到各个线谱多普勒信号s_m(n),m＝1,2,…,M。将分离出的线谱多普勒信号s_m(n)复数化得到z_m(n)，其中z_m(n)＝s_m(n)+jH[s_m(n)]，H[s_m(n)]表示信号s_m(n)的希尔伯特变换。

2）对已得到的多普勒解析信号z_m(n)进行截取，截取后的多普勒信号时长为T′且以预估计的正横时刻为截取时间段的中点，其中T′∈[20,40]秒，典型值取为30秒。利用通过特性法可获得正横时刻预估计值

按公式对序列z_m(n)截取便得到总点数为N′的多普勒信号时间序列zz_m(n),n＝1,2,…,N′，N′＝T′F_S。多普勒信号zz_m(n)等效的高阶多项式相位信号的阶数I为整数，取值范围为[3,5]，典型值取为4，利用HAF算法估计信号zz_m(n)的瞬时相位和时变幅度，具体步骤实施如下：(a)构造高阶瞬时矩 ${\mathrm{DP}}_I[{\mathrm{zz}}_m\left(n\right),\mathrm{\τ}]=\underset{q=0}{\overset{I-1}{\mathrm{\Π}}}{\left[{\mathrm{zz}}_m^{{\$}_q}(n-\mathrm{q\τ})\right]}^{\left(\begin{array}{c}I-1\\ q\end{array}\right)},$ 其中 $\mathrm{\τ}=\frac{N^{\′}}{I}$ 为时间延迟，I为信号的阶数，

(b)对步骤(a)中的高阶瞬时矩进行离散傅立叶变换，得到信号zz_m(n)的高阶模糊函数：

ω为数字角频率，通过频域的一维峰值搜索可得多项式相位系数 ${\hat{a}}_m=\frac{1}{I!{\left(\mathrm{\τ\Δ}\right)}^{I-1}}\underset{\mathrm{\ω}}{\arg \max }\left\{\right|\mathrm{HAF}\left|\right\};$ (c)令 ${\mathrm{zz}}_m\left(n\right)={\mathrm{zz}}_m\left(n\right)\exp \{-j{\hat{a}}_m{\left(\mathrm{n\Δ}\right)}^I\},$ 并使I＝I-1，如果I≥1，返回步骤(a)，直至估计得到全部的多项式相位系数

i＝1,2,…I;(d)计算

的相位角可以得到

这样，多普勒信号zz_m(n)的瞬时相位估计值

其中时间

(e)将接收到的多普勒信号z_m(n)乘上对乘积序列的实部进行低通滤波，可以获得多普勒信号zz_m(n)的时变幅度估计值 ${\hat{A}}_m\left(t_m\right),t_m\∈(t_{m0}^{\′}-T^{\′}/2,t_{m0}^{\′}+T^{\′}/2].$

3）多普勒信号zz_m(n)的瞬时相位模型为 ${\mathrm{\φ}}_m\left(t\right)=2\mathrm{\π}{f}_{m0}\frac{c^2}{c^2-v^2}[t-\frac{v^2}{c^2}{t}_{m0}-\frac{1}{c^2}\sqrt{v^2{c}^2{(t-t_{m0})}^2+R_0^2(c^2-v^2)}]+{\mathrm{\φ}}_{m0},$ 通过多线谱合成的加权非线性最小二乘算法可以得到舰船的航行速度估计值

和舰船通过测量水听器时的正横距离估计值

$\{\hat{v},{\hat{R}}_0,{\hat{f}}_{m0},{\hat{t}}_{m0},{\widehat{\mathrm{\φ}}}_{m0}\}=\underset{(v,R_0,f_{m0},t_{m0},{\mathrm{\φ}}_{m0})}{\arg \min }\left\{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{k=1}{\overset{N^{\′}}{\mathrm{\Σ}}}{\hat{A}}_m^2{[{\mathrm{\φ}}_m(t_{\mathrm{mk}};v,R_0,f_{m0},t_{m0},{\mathrm{\φ}}_{m0})-{\widehat{\mathrm{\φ}}}_m\left(t_{\mathrm{mk}}\right)]}^2\right\},$ k＝1,2,…,N′为时间t_m的标记序号。令向量变量x＝[v;R₀;f_m0;t_m0;φ_m0]，假设最小化目标函数 $Q\left(x\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{k=1}{\overset{N^{\′}}{\mathrm{\Σ}}}{\hat{A}}_m^2{[{\mathrm{\φ}}_m(t_{\mathrm{mk}};v,R_0,f_{m0},t_{m0},{\mathrm{\φ}}_{m0})-{\widehat{\mathrm{\φ}}}_m\left(t_{\mathrm{mk}}\right)]}^2,$ 下面给出

和的具体计算步骤如下：(a)设定步骤(1)中得到的为各线谱声源频率的初值

(b)将通过特性法获得的正横时刻预估计值

作为正横时刻的初值

(c)由于初始相位为常数，其取值范围为[0,2π]，其值对迭代实现加权非线性最小二乘法的收敛影响不大，设定其初值

为π；(d)根据先验知识可知舰船航行速度范围一般为4kn到40kn，则航行速度的初始值取为中间值v⁰＝22kn；(e)对于步骤(a)、(b)、(c)和(d)给定的

和v⁰，确定舰船通过水听器时的正横距离初值

(f)容易获得最小化目标函数Q(x)对向量变量x的一阶偏导数为

再利用Levenberg-Marquardt迭代算法求解得到各参数的估计值

和

其迭代关系式为 $x^{(p+1)}=x^{\left(p\right)}-\frac{{\left[J\left(x^{\left(p\right)}\right)\right]}^{T}Q\left(x^{\left(p\right)}\right)}{{\left[J\left(x^{\left(p\right)}\right)\right]}^{T}J\left(x^{\left(p\right)}\right)+\mathrm{\μI}},$ x^(p)表示第p次迭代变量x的值，μ为收敛因子，I为单位矩阵，“T”为矩阵转置符号，迭代初始值

为步骤(a)、(b)、(c)、(d)和(e)所得。

Claims (1)

1.一种舰船速度和距离联合测量的被动声多普勒相位方法，其特征在于包括下述步

骤：

1）测量被测舰船匀速直线接近测量水听器并离开测量水听器水听器的过程中接收舰船辐射的声信号，并将其转换成电压信号s(t),0＜t≤T，通过前放后用数据采集仪记录下来，t表示记录时间，T表示数据记录的时间长度；

2）把采集得到的水听器接收信号s(t)做快速傅立叶变换分析频谱，找到信号中存在的信噪比大于或等于5dB的低频线谱，其频率为

m＝1,2,…,M，M表示低频线谱的个数；通过带通滤波器分离出存在线谱的频带宽度为

各个小区域信号，得到各个低频线谱多普勒信号s_m(t)；将分离出的线谱多普勒信号s_m(t)复数化得到多普勒解析信号z_m(t)，其中z_m(t)＝s_m(t)+jH[s_m(t)]，H[s_m(t)]表示信号s_m(t)的希尔伯特变换；

3）多普勒信号的表达式为z_m(t)＝A_m(t)exp[jφ_m(t)]，其中，时变幅度 $A_m\left(t\right)=\frac{A_{m0}(c^2-v^2)}{c[\sqrt{{c^2{v}^2(t-t_{m0})}^2+R_0^2(c^2-v^2)}-v^2(t-t_{m0})]},$ 瞬时相位 ${\mathrm{\φ}}_m\left(t\right)=2\mathrm{\π}{f}_{m0}\frac{c^2}{c^2-v^2}[t-\frac{v^2}{c^2}{t}_{m0}-\frac{1}{c^2}\sqrt{v^2{c}^2{(t-t_{m0})}^2+R_0^2(c^2-v^2)}]+{\mathrm{\φ}}_{m0},$ A_m0表示舰船各辐射线谱声信号的幅度，f_m0表示各辐射线谱声信号的频率，t_m0为各线谱声源通过正横位置离测量水听器最近时的正横时刻，φ_m0为多普勒信号z_m(t)的初始相位；将多普勒信号z_m(t)表示为阶数为I的高阶多项式相位信号

其中a_mi,i＝0,1,…,I为多普勒信号z_m(t)的多项式相位系数；利用高阶模糊函数算法估计出多项式相位系数

i＝0,1,…,I，进而得到多普勒信号z_m(t)的瞬时相位估计值

将接收到的多普勒信号z_m(t)乘上

对乘积的实部进行低通滤波，获得多普勒信号z_m(t)的时变幅度估计值

4）根据步骤(3)中多普勒信号z_m(t)的瞬时相位表达式，并利用获得的瞬时相位估计值

和时变幅度估计值

通过多线谱合成的加权非线性最小二乘估计法得到舰船航行速度的估计值、正横距离的估计值、各线谱声源经过测量水听器时的正横时刻估计值、各线谱声源的频率估计值和多普勒信号的初始相位估计值： $\{\hat{v},{\hat{R}}_0,{\hat{f}}_{m0},{\hat{t}}_{m0},{\widehat{\mathrm{\φ}}}_{m0}\}=\underset{(v,R_0,f_{m0},t_{m0},{\mathrm{\φ}}_{m0})}{\arg \min }\left\{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\∫}_0^T{\hat{A}}_m^2{[{\mathrm{\φ}}_m\left(t\right)-{\widehat{\mathrm{\φ}}}_m\left(t\right)]}^2\mathrm{dt}\right\}.$

Images