CN103076604A - 一种基于频散特征的低频水声脉冲信号距离的测量方法 - Google Patents

Abstract

本发明公开了一种基于频散特征的低频水声脉冲信号距离的测量方法，包括步骤10）单个水听器接收低频水声脉冲信号，然后对该低频水声脉冲信号的频散特征进行具有自适应径向高斯核函数的时频分布进行表征；步骤20）提取时频分布中传播模式的频散曲线，得到模式间的到达时间差，步骤30）根据各频点下提取出的模式之间的时间差，测算声源的距离：首先使用简正波模型测算出所提取出的模式，在各自频点下所对应的理论群速度值；然后在模式m_c和模式n_c相同的频率范围内，估计各频点下声源的距离；最后对各频点所估计出的距离进行算数平均值测算，得到声源的距离。该测量方法能够利用单个水听器实现浅海低频水声脉冲信号距离的准确测量。

Description

一种基于频散特征的低频水声脉冲信号距离的测量方法

技术领域

本发明涉及针对水声信号处理技术领域的水下目标距离的测量方法，具体来说，涉及一种基于频散特征的低频水声脉冲信号距离的测量方法。

背景技术

测定水下目标的距离在水声技术中扮演着极其重要的角色，是水声技术中的一个重要且基本的问题，同时也是国民经济建设和国防建设的关键技术。传统的水声测距技术一般采用阵列的处理方式，具有大的孔径，以获得良好的阵增益和分辨性能。但是采用多阵元的大阵列，一方面增加了系统的开销，给基阵的设计带来不便；另一方面，在实际海水中布放时会受到诸如阵倾斜以及阵元失效等问题。此外，在某些应用方面，由于受到安装平台的限制，使得多阵元的布放无法实现。因此，若能利用单一水听器来测定水下目标的距离，就能很好地解决阵列处理时所带来的诸多不便。

对于在浅海环境中传播的低频水声脉冲信号而言，受到海洋媒质的影响，会发生频散现象。频散现象使得信号的不同频率分量会以不同的群速度传播，同时受到海洋吸收和衰减的影响，造成接收波形的失真。频散一方面对发射信号进行了更为复杂的变化，引入了具有非线性时频形状的多分量结构；但另一方面，频散本身也蕴含了信号的距离信息。通过分析频散波导中所接收到的信号，可以有助于我们获取关于信号的距离信息。

对于海洋波导中频散现象的研究，多集中于时频分析的方法。传统的时频分析方法如基于线性时频表示的短时傅里叶变换受不确定原理的影响，其时频分辨率较低；而基于二次型时频表示的魏格纳-威利分布等存在严重的交叉项干扰问题。针对传统时频分析方法的不足，有学者提出了基于信号的自适应径向高斯核函数的时频分析方法，能够很好地反映水声脉冲信号的短时瞬态的非平稳特性，避免了各信号分量之间的交叉项，提高了对于非平稳水声脉冲信号的分辨能力。

发明内容

技术问题：本发明所要解决的技术问题是：提供一种基于频散特征的低频水声脉冲信号距离的测量方法，该测量方法能够利用单个水听器实现浅海低频水声脉冲信号距离的准确测量。

技术方案：为解决上述技术问题，本发明采用的技术方案是：

一种基于频散特征的低频水声脉冲信号距离的测量方法，该测量方法包括以下步骤：

步骤10)单个水听器接收低频水声脉冲信号，然后对该低频水声脉冲信号的频散特征进行具有自适应径向高斯核函数的时频分布进行表征，包括如下的步骤：

步骤101)对单个水听器接收到的低频水声脉冲信号，根据式(1)确定相应的二维频偏-时滞域上的模糊函数A_d(m，n)；

$A_d(m,n)=T_s\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{y}^{*}({\mathrm{kT}}_s-{\mathrm{nT}}_s)y({\mathrm{kT}}_s+{\mathrm{nT}}_s)e^{i2\mathrm{\πmk}/N}$ 式(1)

其中，N表示接收信号的点长；m表示直角坐标系下的离散化的频偏索引；n表示直角坐标系下的离散化的时滞索引；T_s表示采样时间间隔；y表示单水听器所接收到的水声脉冲信号；y^*表示y的共轭；i表示虚数单位，

步骤102)首先设置离散化频偏的总点数P和离散时滞的总点数Q，其中，

Q＝N；然后从r＝0至

产生P个离散化的极径，从ψ＝0至π产生Q个离散化的径向角，r表示极径，ψ表示径向角；

步骤103)以离散化的时滞

和离散化的频偏

采用二维插值的方法，通过式(2)的坐标变换公式，将直角坐标系下的模糊函数A_d(m，n)转换为极坐标系下的模糊函数A_u(p，q)：

$\left\{\begin{array}{c}X=r\·\cos \left(\mathrm{\ψ}\right)\\ Y=r\·\sin \left(\mathrm{\ψ}\right)\end{array}\right.$ 式(2)

其中，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，X表示转换后的频偏坐标，Y表示转换后的时滞坐标；

步骤104)采用迭代算法测算最优扩展函数σ，确定低频水声脉冲信号的时频分布I(t，f)，其中，t表示时频分布的时间，f表示时频分布的频率；

步骤20)提取时频分布中传播模式的频散曲线，得到模式间的到达时间差，包括以下步骤：

步骤201)根据步骤104)所得到的时频分布I(t，f)进行峰值搜索，确定接收信号的中心频率以及相应的瞬时时间，如式(15)所示，

$[t_0,f_0]=\underset{t,f}{\arg \max }I(t,f)$ 式(15)

其中，f₀表示峰值位置处所对应的中心频率，t₀表示相应的瞬时时间；

表示I(t，f)最大值时的时间和频率；

步骤202)由步骤201)确定的峰值所对应的时间t₀和频率f₀，得到时频分布I(t，f)上最高的能量模式m_c，并由m_c在时频分布I(t，f)上所显示的范围，确定其所对应的最低下限频率

和对应时间

以及最高上限频率

和对应时间

步骤203)在模式m_c的时频显示范围内，沿频率方向对时频分布I(t，f)进行最大值搜索，确定每个搜索频率所对应的瞬时时间，从而得到模式m_c所对应的频散曲线的估计值

其中，

表示模式m_c的瞬时时间，

表示模式m_c的瞬时频率；在此过程中，如果满足

则对得到的瞬时时间进行校正，校正如式(16)所示：

${\hat{t}}_{j+1}=\arg \mathrm{sub}\max \left(I(t_l^{m_c}:t_h^{m_c},f_{j+1})\right)$

其中，

表示第j次最大值搜索所得到的瞬时时间，

表示第j+1次最大值搜索所得到的瞬时时间，

表示

的绝对值，S表示时频域上模式间的时间判断门限，

表示时频分布I(t，f)沿频率f_j+1方向从

到所得到的次大值所对应的时间；

步骤204)由模式m_c频散曲线的估计值

确定二值掩模滤波矩阵M(t，f)；

步骤205)将二值掩模滤波矩阵M(t，f)进行膨胀处理，并将膨胀之后的掩模滤波矩阵M′(t，f)与时频分布矩阵I(t，f)相乘，得到掩模后的时频分布矩阵I′(t，f)；

步骤206)对掩模后的时频分布I′(t，f)进行峰值搜索，确定掩模后的时频分布I′(t，f)峰值所对应的时间和频率，从而得到时频分布I(t，f)上次高的能量模式n_c，对掩模后的时频分布I′(t，f)执行同步骤203)相同的方法，得到模式m_c所对应频散曲线的估计值

其中，

表示模式n_c的瞬时时间，

表示模式n_c的瞬时频率；

步骤207)测算模式m_c和模式n_c估计得到的频散曲线各频点下的时间差ΔT(f)，若 ${\hat{t}}^{m_c}>{\hat{t}}^{n_c},$ 则 $\mathrm{\&Delta;T}\left(f\right)={\hat{t}}^{m_c}-{\hat{t}}^{n_c};$ 若 ${\hat{t}}^{m_c}<{\hat{t}}^{n_c},$ 则 $\mathrm{\&Delta;T}\left(f\right)={\hat{t}}^{n_c}-{\hat{t}}^{m_c};$

步骤30)根据步骤20)得到的各频点下提取出的模式之间的时间差，测算声源的距离，包括如下的步骤：

步骤301)使用简正波模型kraken测算出所提取出的模式，在各自频点下所对应的理论群速度值

和

其中，

表示模式m_c在频率f处的群速度值，

表示模式n_c在频率f处的群速度值；

步骤302)在模式m_c和模式_nc相同的频率范围内，采用式(16)估计各频点下声源的距离：

$\hat{r}\left(f\right)=\frac{v_g^{m_c}\left(f\right)v_g^{n_c}\left(f\right)}{|v_g^{n_c}\left(f\right)-v_g^{m_c}\left(f\right)|}\&CenterDot;\mathrm{\&Delta;T}\left(f\right)$ 式(16)

其中，

表示频点f处的声源距离的估计值，

表示与

之差的绝对值；

步骤303)对各频点所估计出的距离

进行算数平均值测算，得到声源的距离

进一步，所述的步骤104)包括以下步骤：

步骤1041)设置常量参数及初值：收敛因子ε、最大迭代次数λ、初始步长μ(0)、初始扩展向量σ(0)、步长控制因子Δ、常数ρ₁和ρ₂，0＜ρ₁＜ρ₂＜1，Δ＞0，

其中，α表示核函数体积控制参数，Q表示离散时滞的总点数，Δ_ψ表示径向角分辨率，[1，…，1]表示1构成的向量，T表示转置；

步骤1042)根据式(3)确定初始代价函数f(σ(0))，根据式(4)确定初始代价函数的梯度向量

$f\left(\mathrm{\&sigma;}\left(0\right)\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}p{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left(p{\mathrm{\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(0\right)}}$ 式(3)

其中，P表示离散频偏的总点数，Q表示离散时滞的总点数，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，Δ_r表示极径分辨率，σ_q(0)表示初始扩展向量σ(0)的第q个采样值，且σ_q(0)＝σ(qΔ_ψ)，σ表示最优扩展函数；

$\&dtri;f\left(0\right)={[\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&sigma;}}_1\left(0\right)},...,\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&sigma;}}_Q\left(0\right)}]}^T$ 式(4)

其中，

表示偏导数，上标T表示向量的转置；

步骤1043)进行迭代操作，测算k+1次迭代的扩展向量，如式(5)所示，并进行归一化处理，如式(6)所示

$\mathrm{\&sigma;}(k+1)=\mathrm{\&sigma;}\left(k\right)+\mathrm{\&mu;}\left(k\right)\&dtri;f\left(k\right)$ 式(5)

$\mathrm{\&sigma;}(k+1)\&LeftArrow;\frac{\mathrm{\&sigma;}(k+1)}{\left|\right|\mathrm{\&sigma;}(k+1)\left|\right|}\sqrt{\frac{2\mathrm{\&pi;\&alpha;}}{{\mathrm{\&Delta;}}_{\mathrm{\&psi;}}}}$ 式(6)

其中，σ(k+1)表示第k+1次迭代的扩展向量，σ(k)表示第k次迭代的扩展向量，μ(k)表示第k次的迭代步长，

表示第k次迭代的梯度向量，←表示更新，||σ(k+1)||表示第k+1次迭代的扩展向量σ(k+1)的范数；

根据式(7)更新每次迭代的代价函数f(σ(k+1))，根据式(8)确定每次迭代的梯度向量

$f\left(\mathrm{\&sigma;}(k+1)\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}p{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(k\right)}}$ 式(7)

$\&dtri;f\left(k\right)=\frac{{2\mathrm{\&Delta;}}_r^2}{{\mathrm{\&sigma;}}_q^3\left(k\right)}\underset{p=1}{\overset{P-1}{\mathrm{\&Sigma;}}}{p}^3{\left|A_u(p.q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(k\right)}}$ 式(8)

其中，k表示迭代次数，σ_q(k)为第k次迭代时扩展向量σ(k+1)的第q个采样值；更新迭代步长，若满足 $f\left(\mathrm{\&sigma;}(k+1)\right)-f\left(\mathrm{\&sigma;}\left(k\right)\right)\&le;{\mathrm{\&rho;}}_1\mathrm{\&mu;}\left(k\right){\&dtri;f}^T\left(k\right)(\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)),$ 则步长为μ(k+1)＝μ(k)/Δ；若满足 $f\left(\mathrm{\&sigma;}(k+1)\right)-f\left(\mathrm{\&sigma;}\left(k\right)\right)>{\mathrm{\&rho;}}_2\mathrm{\&mu;}\left(k\right){\&dtri;f}^T\left(k\right)(\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)),$ 则步长为μ(k+1)＝μ(k)Δ；

如果迭代次数达到预先的最大迭代次数λ，或者同时满足式(9)和式(10)，则终止迭代，得到的扩展向量σ；

代价函数约束：f(σ(k+1))-f(σ(k))＜εf(σ(k))   式(9)

扩展向量约束： $\left|\right|\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)\left|\right|<\mathrm{\&epsiv;}\sqrt{2\mathrm{\&pi;\&alpha;}/{\mathrm{\&Delta;}}_{\mathrm{\&psi;}}}$ 式(10)；

步骤1044)根据式(11)确定直角坐标系下的径向角ψ_d：

${\mathrm{\&psi;}}_d=\arctan \frac{\mathrm{\&tau;}}{\mathrm{\&theta;}}$ 式(11)

其中，τ表示离散化的时滞，θ表示离散化的频偏；

将步骤1043)得到的扩展向量σ和径向角ψ扩充到整个模糊域，得到式(12)和式(13)

ψ←[ψ(1:Q)，ψ(2:Q)+π]  式(12)

σ←[σ(1:Q)，σ(2:Q)]     式(13)

由ψ和σ采用插值的方式，测算出直角坐标系下径向角ψ_d对应的扩展函数σ_d，并由直角坐标系下的极径公式得到直角坐标系下的径向高斯核函数Φ_d(m，n)，如式(14)所示，

${\mathrm{\&Phi;}}_d(m,n)=e^{-r^2/{2\mathrm{\&sigma;}}^2\left(\mathrm{\&psi;}\right)}$ 式(14)

步骤1045)将模糊函数A_d与径向高斯核函数Φ_d相乘，并利用二维傅里叶变换，得到信号的时频分布I(t，f)。

进一步，所述的步骤204)包括以下过程：当

且

时，M(t，f)＝0；当且

时，M(t，f)＝1；其中，表示时频分布I(t，f)中时间t位于的区间，

表示时频分布I(t，f)中频率f位于的区间，

表示时间t不位于区间

内，

表示频率f不位于区间内。

有益效果：与现有技术相比，本发明具有以下显著优点：

(1)时频分辨率高，各阶传播模式的辨识更加容易。现有技术多采用固定核函数的时频分析方法，其时频分辨率较低，频散表征不是十分理想。本发明采用了自适应径向高斯核函数的时频分析技术来表征低频水声脉冲信号中的频散特征，具有较高的时频分辨率，能够很好地反映低频水声脉冲信号自身的频散特征，使得传播的各阶模式能够更好地在时频域进行辨识。

(2)测距简单，运算量低。现有技术多采用多阵元的处理方法且常采用匹配场搜索的方式，而本发明只利用了单个水听器来求解低频水声脉冲信号的距离。本发明通过在时频域中提取低频水声脉冲信号中的频散特征，获得传播模式的到达时间差的方法来估计声源的距离，避免了多阵元处理时的诸多不便以及匹配场算法所需的网格搜索计算，降低了运算量。

附图说明

图1是本发明实施例采用的声学环境示意图。

图2是本发明实施例接收信号的时频表示图。

图3是本发明实施例估计出的频散曲线图。

具体实施方式

下面将结合本发明实施例和附图，对本发明的技术方案进行清楚、完整地描述。

本发明的基于频散特征的低频水声脉冲信号距离的测量方法，包括以下步骤：

步骤10)单个水听器接收低频水声脉冲信号，然后对该低频水声脉冲信号的频散特征进行具有自适应径向高斯核函数的时频分布进行表征。步骤10)具体包括步骤101)至步骤104)：

步骤101)对单个水听器接收到的低频水声脉冲信号，根据式(1)确定相应的二维频偏-时滞域上的模糊函数A_d(m，n)；

$A_d(m,n)=T_s\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{y}^{*}({\mathrm{kT}}_s-{\mathrm{nT}}_s)y({\mathrm{kT}}_s+{\mathrm{nT}}_s)e^{i2\mathrm{\&pi;mk}/N}$ 式(1)

其中，N表示接收信号的点长；m表示直角坐标系下的离散化的频偏索引；n表示直角坐标系下的离散化的时滞索引；T_s表示采样时间间隔；y表示单水听器所接收到的水声脉冲信号；y^*表示y的共轭；i表示虚数单位，

步骤102)首先设置离散化频偏的总点数P和离散时滞的总点数Q，其中，Q＝N；然后从r＝0至

产生P个离散化的极径，从ψ＝0至π产生Q个离散化的径向角，r表示极径，ψ表示径向角。

步骤103)以离散化的时滞

和离散化的频偏

采用二维插值的方法，通过式(2)的坐标变换公式，将直角坐标系下的模糊函数A_d(m，n)转换为极坐标系下的模糊函数A_u(p，q)：

$\left\{\begin{array}{c}X=r\&CenterDot;\cos \left(\mathrm{\&psi;}\right)\\ Y=r\&CenterDot;\sin \left(\mathrm{\&psi;}\right)\end{array}\right.$ 式(2)

其中，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，X表示转换后的频偏坐标，Y表示转换后的时滞坐标。

步骤104)采用迭代算法测算最优扩展函数σ，确定低频水声脉冲信号的时频分布I(t，f)，其中，t表示时频分布的时间，f表示时频分布的频率。步骤104)具体包括以下的步骤1041)-步骤1045)。

步骤1041)设置常量参数及初值：收敛因子ε、最大迭代次数λ、初始步长μ(0)、初始扩展向量σ(0)、步长控制因子Δ、常数ρ₁和ρ₂，0＜ρ₁＜ρ₂＜1，Δ＞0，

其中，α表示核函数体积控制参数，α的取值范围优选1≤α≤5，Q表示离散时滞的总点数，Δ_ψ表示径向角分辨率，[1，…，1]表示1构成的向量，T表示转置。

步骤1042)根据式(3)确定初始代价函数f(σ(0))，根据式(4)确定初始代价函数的梯度向量

$f\left(\mathrm{\&sigma;}\left(0\right)\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}p{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left(p{\mathrm{\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(0\right)}}$ 式(3)

其中，P表示离散频偏的总点数，Q表示离散时滞的总点数，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，Δ_r表示极径分辨率，σ_q(0)表示初始扩展向量σ(0)的第q个采样值，且σ_q(0)＝σ(qΔ_ψ)，σ表示最优扩展函数；

$\&dtri;f\left(0\right)={[\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&sigma;}}_1\left(0\right)},...,\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&sigma;}}_Q\left(0\right)}]}^T$ 式(4)

其中，

表示偏导数，上标T表示向量的转置。

在初始梯度向量中各项的计算式为：

$\frac{\&PartialD;f}{{\&PartialD;\mathrm{\&sigma;}}_q\left(0\right)}=\frac{{2\mathrm{\&Delta;}}_r^2}{{\mathrm{\&sigma;}}_q^3\left(0\right)}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}{p}^3{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(0\right)}}$

其中，q＝1、2、…、Q。

步骤1043)进行迭代操作，测算k+1次迭代的扩展向量，如式(5)所示，并进行归一化处理，如式(6)所示：

$\mathrm{\&sigma;}(k+1)=\mathrm{\&sigma;}\left(k\right)+\mathrm{\&mu;}\left(k\right)\&dtri;f\left(k\right)$ 式(5)

$\mathrm{\&sigma;}(k+1)\&LeftArrow;\frac{\mathrm{\&sigma;}(k+1)}{\left|\right|\mathrm{\&sigma;}(k+1)\left|\right|}\sqrt{\frac{2\mathrm{\&pi;\&alpha;}}{{\mathrm{\&Delta;}}_{\mathrm{\&psi;}}}}$ 式(6)

其中，σ(k+1)表示第k+1次迭代的扩展向量，σ(k)表示第k次迭代的扩展向量，μ(k)表示第k次的迭代步长，

表示第k次迭代的梯度向量，←表示更新，||σ(k+1)||表示第k+1次迭代的扩展向量σ(k+1)的范数；

根据式(7)更新每次迭代的代价函数f(σ(k+1))，根据式(8)确定每次迭代的梯度向量

$f\left(\mathrm{\&sigma;}(k+1)\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}p{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(k\right)}}$ 式(7)

$\&dtri;f\left(k\right)=\frac{{2\mathrm{\&Delta;}}_r^2}{{\mathrm{\&sigma;}}_q^3\left(k\right)}\underset{p=1}{\overset{P-1}{\mathrm{\&Sigma;}}}{p}^3{\left|A_u(p.q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(k\right)}}$ 式(8)

其中，k表示迭代次数，σ_q(k)为第k次迭代时扩展向量σ(k+1)的第q个采样值；

更新迭代步长，若满足 $f\left(\mathrm{\&sigma;}(k+1)\right)-f\left(\mathrm{\&sigma;}\left(k\right)\right)\&le;{\mathrm{\&rho;}}_1\mathrm{\&mu;}\left(k\right){\&dtri;f}^T\left(k\right)(\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)),$ 则步长为μ(k+1)＝μ(k)/Δ；若满足 $f\left(\mathrm{\&sigma;}(k+1)\right)-f\left(\mathrm{\&sigma;}\left(k\right)\right)>{\mathrm{\&rho;}}_2\mathrm{\&mu;}\left(k\right){\&dtri;f}^T\left(k\right)(\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)),$ 则步长为μ(k+1)＝μ(k)Δ；

如果迭代次数达到预先的最大迭代次数λ，或者同时满足式(9)和式(10)，则终止迭代，得到的扩展向量σ；

代价函数约束：f(σ(k+1))-f(σ(k))＜εf(σ(k))   式(9)

扩展向量约束： $\left|\right|\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)\left|\right|<\mathrm{\&epsiv;}\sqrt{2\mathrm{\&pi;\&alpha;}/{\mathrm{\&Delta;}}_{\mathrm{\&psi;}}}$ 式(10)。

步骤1044)根据式(11)确定直角坐标系下的径向角ψ_d：

${\mathrm{\&psi;}}_d=\arctan \frac{\mathrm{\&tau;}}{\mathrm{\&theta;}}$ 式(11)

其中，τ表示离散化的时滞，θ表示离散化的频偏；

将步骤1043)得到的扩展向量σ和径向角ψ扩充到整个模糊域，得到式(12)和式(13)

ψ←[ψ(1:Q)，ψ(2:Q)+π]   式(12)

σ←[σ(1:Q)，σ(2:Q)]      式(13)

由ψ和σ采用插值的方式，测算出直角坐标系下径向角ψ_d对应的扩展函数σ_d，并由直角坐标系下的极径公式

得到直角坐标系下的径向高斯核函数Φ_d(m，n)，如式(14)所示，

${\mathrm{\&Phi;}}_d(m,n)=e^{-r^2/{2\mathrm{\&sigma;}}^2\left(\mathrm{\&psi;}\right)}$ 式(14)

步骤1045)将模糊函数A_d与径向高斯核函数Φ_d相乘，并利用二维傅里叶变换，得到信号的时频分布I(t，f)。

步骤20)提取时频分布中传播模式的频散曲线，得到模式间的到达时间差。步骤20)具体包括以下步骤201)-步骤207)。

步骤201)根据步骤104)所得到的时频分布I(t，f)进行峰值搜索，确定接收信号的中心频率以及相应的瞬时时间，如式(15)所示，

$[t_0,f_0]=\underset{t,f}{\arg \max }I(t,f)$ 式(15)

其中，f₀表示峰值位置处所对应的中心频率，t₀表示相应的瞬时时间；

表示I(t，f)最大值时的时间和频率。

步骤202)由步骤201)确定的峰值所对应的时间t₀和频率f₀，得到时频分布I(t，f)上最高的能量模式m_c，并由m_c在时频分布I(t，f)上所显示的范围，确定其所对应的最低下限频率

和对应时间

以及最高上限频率

和对应时间

步骤203)在模式m_c的时频显示范围内，沿频率方向对时频分布I(t，f)进行最大值搜索，确定每个搜索频率所对应的瞬时时间，从而得到模式m_c所对应的频散曲线的估计值

其中，

表示模式m_c的瞬时时间，

表示模式m_c的瞬时频率；在此过程中，如果满足

则对得到的瞬时时间进行校正，校正如式(16)所示：

${\hat{t}}_{j+1}=\arg \mathrm{sub}\max \left(I(t_l^{m_c}:t_h^{m_c},f_{j+1})\right)$

其中，

表示第j次最大值搜索所得到的瞬时时间，

表示第j+1次最大值搜索所得到的瞬时时间，表示

的绝对值。S表示时频域上模式间的时间判断门限。S根据时频分布所显示出的模式选择所需的时间判断门限，S可为采样时间间隔T_s的10至50倍。

表示时频分布I(t，f)沿频率f_j+1方向从

到

所得到的次大值所对应的时间。

步骤204)由模式m_c频散曲线的估计值

确定二值掩模滤波矩阵M(t，f)。

$M(t,f)=\left\{\begin{array}{cc}0,& t\&Element;[t_l^m,t_h^m],f\&Element;[f_l^m,f_h^m]\\ 1,& t\&NotElement;[t_l^m,t_h^m],f\&NotElement;[f_l^m,f_h^m]\end{array}\right.$

在步骤204)中，确定二值掩模滤波矩阵的过程：当

且

时，M(t，f)＝0；当

且

时，M(t，f)＝1；其中，

表示时频分布I(t，f)中时间t位于的区间，

表示时频分布I(t，f)中频率f位于的区间，

表示时间t不位于区间

内，

表示频率f不位于区间

内。

步骤205)将二值掩模滤波矩阵M(t，f)进行膨胀处理，并将膨胀之后的掩模滤波矩阵M′(t，f)与时频分布矩阵I(t，f)相乘，得到掩模后的时频分布矩阵I′(t，f)。

在步骤205)中，膨胀处理为现有技术，具体过程可参见matlab图像工具箱下的函数imdilate()进行求解。

步骤206)对掩模后的时频分布I′(t，f)进行峰值搜索，确定掩模后的时频分布I′(t，f)峰值所对应的时间和频率，从而得到时频分布I(t，f)上次高的能量模式n_c，对掩模后的时频分布I′(t，f)执行同步骤203)相同的方法，得到模式n_c所对应频散曲线的估计值

其中，表示模式n_c的瞬时时间，

表示模式n_c的瞬时频率。

步骤207)测算模式m_c和模式n_c估计得到的频散曲线各频点下的时间差ΔT(f)，若 ${\hat{t}}^{m_c}>{\hat{t}}^{n_c},$ 则 $\mathrm{\&Delta;T}\left(f\right)={\hat{t}}^{m_c}-{\hat{t}}^{n_c};$ 若 ${\hat{t}}^{m_c}<{\hat{t}}^{n_c},$ 则 $\mathrm{\&Delta;T}\left(f\right)={\hat{t}}^{n_c}-{\hat{t}}^{m_c}.$

步骤30)根据步骤20)得到的各频点下提取出的模式之间的时间差，测算声源的距离。步骤30)包括如下步骤301)-步骤303)。

步骤301)使用简正波模型kraken测算出所提取出的模式，在各自频点下所对应的理论群速度值

和

其中，

表示模式m_c在频率f处的群速度值，

表示模式n_c在频率f处的群速度值。

步骤302)在模式m_c和模式n_c相同的频率范围内，采用式(16)估计各频点下声源的距离：

$\hat{r}\left(f\right)=\frac{v_g^{m_c}\left(f\right)-v_g^{n_c}\left(f\right)}{|v_g^{n_c}\left(f\right)-v_g^{m_c}\left(f\right)|}\&CenterDot;\mathrm{\&Delta;T}\left(f\right)$ 式(16)

其中，

表示频点f处的声源距离的估计值，表示

与之差的绝对值。

步骤303)对各频点所估计出的距离进行算数平均值测算，得到声源的距离

本发明的基于频散特征的低频水声脉冲信号距离的测量方法包括：首先确定基于信号的自适应径向高斯核函数的时频分布；通过时频分析提取低频水声脉冲信号特征参数；通过测算模式的群速度和模式之间的到达时间差，估计出声源的距离。本发明技术方案能够实现只借助单一接收水听器完成对低频水声脉冲信号的距离估计。

下面例举一实施例。

如图1所示，为本实施例所用到的声学环境。其中，声源位于水下30m，单一接收水听器位于水下100m，声源与接收机之间的距离为15km，海水深度为120m。水层声速为1500m/s，海水密度为1.0g/cm³，底部声速为1800m/s，密度为1.8g/cm³。发射信号为频带范围为40-120Hz的线性调频信号。图2是本实施例的接收信号的时频分布图。图2的横坐标表示时间，单位为秒(s)，纵坐标表示频率，单位赫兹(Hz)。从图2可以看出：采用自适应径向高斯核函数的时频分布可以在时频图上显示出低频水声脉冲信号中所蕴含的传播模式，其中的箭头表示所要提取的两个传播模式，分别为模式2和模式3。图3为本实施例中模式2和模式3的频散曲线估计结果线条图。图3的横坐标表示时间，单位为秒(s)，纵坐标表示频率，单位赫兹(Hz)。按照本发明的方法测量低频水声脉冲信号距离，具体包括如下步骤：

(1)计算接收信号的模糊函数并将其转换到极坐标系下。

(2)使用优化迭代算法测算最优核函数，本实施例中选择的参数ρ₁＝0.1，ρ₂＝0.9，α＝2，步长控制因子Δ＝10，初始迭代步长为μ(0)＝1，最大迭代次数为λ＝50，收敛因子ε＝10^-5。

(3)由最优核函数测算出极坐标系下的径向高斯核函数，再通过坐标变换转换到直角坐标系下。

(4)将模糊函数与径向高斯核函数的乘积执行二维傅里叶变换，从而得到径向高斯核函数所表示的接收信号的时频分布显示。

从图3的时频分布图中可以看出，径向高斯核函数具有较高的时频分辨率，可准确地反映出低频水声脉冲信号中的传播模式。

(5)对所得的时频分布结果进行最大值搜索，获得确定时频分布关于最大值所在模式，在本实施例中为模式2，如图3所示。

(6)对模式2在其所显示的范围内，按频率方向进行搜索，获得相应的频散曲线的估计值。

(7)利用二值掩模滤波矩阵，对模式2进行过滤，对过滤之后的时频图进行峰值搜索，确定出能量次高的模式，本实施例中为模式3，如图3所示。

(8)采用同提取模式2相同的方法，确定模式3所对应的频散曲线的估计值。

(9)在模式2和模式3相同的公共频带内，测算出两个模式的到达时间差；并由声场理论测算出两个模式各自的群速度，根据本发明所提供的公式测算出各频点下的距离值，并对结果求平均值，得到距离的估计值为15.2732km。

由距离估计结果可以看出，本发明所采用的自适应高斯径向核函数的时频分布能够较为准确地反映低频水声脉冲信号的频散特征，具有较高的时频分辨率，所得到的距离估计结果较为准确。

上述实施例对本发明的基于频散特征的低频水声脉冲信号距离的测量方法进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，不应理解为对本发明的限制。

Claims (5)

1.一种基于频散特征的低频水声脉冲信号距离的测量方法，其特征在于：该测量方法包括以下步骤：

步骤10）单个水听器接收低频水声脉冲信号，然后对该低频水声脉冲信号的频散特征进行具有自适应径向高斯核函数的时频分布进行表征，包括如下的步骤：

步骤101）对单个水听器接收到的低频水声脉冲信号，根据式（1）确定相应的二维频偏—时滞域上的模糊函数A_d(m,n)；

$A_d(m,n)=T_s\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{y}^{*}({\mathrm{kT}}_s-{\mathrm{nT}}_s)y({\mathrm{kT}}_s+{\mathrm{nT}}_s)e^{i2\mathrm{\&pi;mk}/N}$ 式（1）

其中，N表示接收信号的点长；m表示直角坐标系下的离散化的频偏索引；n表示直角坐标系下的离散化的时滞索引；T_S表示采样时间间隔；y表示单水听器所接收到的水声脉冲信号；y^*表示y的共轭；i表示虚数单位，

步骤102）首先设置离散化频偏的总点数P和离散时滞的总点数Q，其中，

Q=N；然后从r=0至

产生P个离散化的极径，从y=0至π产生Q个离散化的径向角，r表示极径，y表示径向角；

步骤103)以离散化的时滞和离散化的频偏

采用二维插值的方法，通过式（2）的坐标变换公式，将直角坐标系下的模糊函数A_d(m,n)转换为极坐标系下的模糊函数A_u(p,q)：

$\left\{\begin{array}{c}X=r\&CenterDot;\cos \left(\mathrm{\&psi;}\right)\\ Y=r\&CenterDot;\sin \left(\mathrm{\&psi;}\right)\end{array}\right.$ 式（2）

其中，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，X表示转换后的频偏坐标，Y表示转换后的时滞坐标；

步骤104）采用迭代算法测算最优扩展函数σ，确定低频水声脉冲信号的时频分布I(t,f)，其中，t表示时频分布的时间，f表示时频分布的频率；

步骤20）提取时频分布中传播模式的频散曲线，得到模式间的到达时间差，包括以下步骤：

步骤201）根据步骤104）所得到的时频分布I(t,f)进行峰值搜索，确定接收信号的中心频率以及相应的瞬时时间，如式（15）所示，

$[t_0,f_0]=\underset{t,f}{\arg \max }I(t,f)$ 式（15）

其中，f₀表示峰值位置处所对应的中心频率，t₀表示相应的瞬时时间；

表示I(t,f)最大值时的时间和频率；

步骤202）由步骤201）确定的峰值所对应的时间t₀和频率f₀，得到时频分布I(t,f)上最高的能量模式m_c，并由m_c在时频分布I(t,f)上所显示的范围，确定其所对应的最低下限频率

和对应时间以及最高上限频率

和对应时间

步骤203）在模式m_c的时频显示范围内，沿频率方向对时频分布I(t,f)进行最大值搜索，确定每个搜索频率所对应的瞬时时间，从而得到模式m_c所对应的频散曲线的估计值

其中，

表示模式m_c的瞬时时间，

表示模式m_c的瞬时频率；在此过程中，如果满足

则对得到的瞬时时间进行校正，校正如式（16）所示：

${\hat{t}}_{j+1}=\arg \mathrm{sub}\max \left(I(t_l^{m_c}:t_h^{m_c},f_{j+1})\right)$

其中，表示第j次最大值搜索所得到的瞬时时间，

表示第j+1次最大值搜索所得到的瞬时时间，表示

的绝对值，S表示时频域上模式间的时间判断门限，

表示时频分布I(t,f)沿频率f_j+1方向从

到

所得到的次大值所对应的时间；

步骤204）由模式m_c频散曲线的估计值

确定二值掩模滤波矩阵M(t,f)；

步骤205）将二值掩模滤波矩阵M(t,f)进行膨胀处理，并将膨胀之后的掩模滤波矩阵M'(t,f)与时频分布矩阵I(t,f)相乘，得到掩模后的时频分布矩阵I′(t,f)；

步骤206）对掩模后的时频分布I′(t,f)进行峰值搜索，确定掩模后的时频分布I′(t,f)峰值所对应的时间和频率，从而得到时频分布I(t,f)上次高的能量模式n_c，对掩模后的时频分布I′(t,f)执行同步骤203）相同的方法，得到模式n_c所对应频散曲线的估计值

其中，表示模式n_c的瞬时时间，

表示模式n_c的瞬时频率；

步骤207）测算模式m_c和模式n_c估计得到的频散曲线各频点下的时间差ΔT(f)，若 ${\hat{t}}^{m_c}>{\hat{t}}^{n_c},$ 则 $\mathrm{\&Delta;T}\left(f\right)={\hat{t}}^{m_c}-{\hat{t}}^{n_c};$ 若 ${\hat{t}}^{m_c}<{\hat{t}}^{n_c},$ 则 $\mathrm{\&Delta;T}\left(f\right)={\hat{t}}^{n_n}-{\hat{t}}^{m_c};$

步骤30）根据步骤20）得到的各频点下提取出的模式之间的时间差，测算声源的距离，包括如下的步骤：

步骤301）使用简正波模型kraken测算出所提取出的模式，在各自频点下所对应的理论群速度值

和

其中，

表示模式m_c在频率f处的群速度值，

表示模式n_c在频率f处的群速度值；

步骤302）在模式m_c和模式n_c相同的频率范围内，采用式（16）估计各频点下声源的距离：

$\hat{r}\left(f\right)=\frac{v_g^{m_c}\left(f\right)v_g^{n_c}\left(f\right)}{|v_g^{n_c}\left(f\right)-v_g^{m_c}\left(f\right)|}\&CenterDot;\mathrm{\&Delta;T}\left(F\right)$ 式（16）

其中，

表示频点f处的声源距离的估计值，表示

与

之差的绝对值；

步骤303）对各频点所估计出的距离

进行算数平均值测算，得到声源的距离

2.按照权利要求1所述的基于频散特征的低频水声脉冲信号距离的测量方法，其特征在于，所述的步骤104）包括以下步骤：

步骤1041）设置常量参数及初值：收敛因子ε、最大迭代次数λ、初始步长μ(0)、初始扩展向量σ(0)、步长控制因子Δ、常数ρ₁和ρ₂，0<ρ₁<ρ₂<1，Δ>0，

其中，α表示核函数体积控制参数，Q表示离散时滞的总点数，Δ_ψ表示径向角分辨率，[1,…,1]表示1构成的向量，T表示转置；

步骤1042）根据式（3）确定初始代价函数f(σ(0))，根据式（4）确定初始代价函数的梯度向量

$f\left(\mathrm{\&sigma;}\left(0\right)\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}p{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left(p{\mathrm{\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(0\right)}}$ 式（3）

其中，P表示离散频偏的总点数，Q表示离散时滞的总点数，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，Δ_r表示极径分辨率，σ_q(0)表示初始扩展向量σ(0)的第q个采样值，且σ_q(0)＝σ(qΔ_ψ)，σ表示最优扩展函数；

$\&dtri;f\left(0\right)={[\frac{\&PartialD;f}{\&PartialD;{\mathrm{\&sigma;}}_1\left(0\right)},...,\frac{\&PartialD;f}{\&PartialD;{\mathrm{\&sigma;}}_Q\left(0\right)}]}^T$ 式（4）

其中，

表示偏导数，上标T表示向量的转置；

步骤1043）进行迭代操作，测算k+1次迭代的扩展向量，如式（5）所示，并进行归一化处理，如式（6）所示

$\mathrm{\&sigma;}(k+1)=\mathrm{\&sigma;}\left(k\right)+\mathrm{\&mu;}\left(k\right)\&dtri;f\left(k\right)$ 式（5）

$\mathrm{\&sigma;}(k+1)\&LeftArrow;\frac{\mathrm{\&sigma;}(k+1)}{\left|\right|\mathrm{\&sigma;}(k+1)\left|\right|}\sqrt{\frac{2\mathrm{\&pi;\&alpha;}}{{\mathrm{\&Delta;}}_{\mathrm{\&psi;}}}}$ 式（6）

其中，σ(k+1)表示第k+1次迭代的扩展向量，σ(k)表示第k次迭代的扩展向量，μ(k)表示第k次的迭代步长，

表示第k次迭代的梯度向量，←表示更新，||σ(k+1)||表示第k+1次迭代的扩展向量σ(k+1)的范数；

根据式（7）更新每次迭代的代价函数f(σ(k+1))，根据式（8）确定每次迭代的梯度向量

$f\left(\mathrm{\&sigma;}(k+1)\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}\underset{p=2}{\overset{P}{\mathrm{\&Sigma;}}}p{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left(p{\mathrm{\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(k\right)}}$ 式（7）

$\&dtri;f\left(k\right)=\frac{2{\mathrm{\&Delta;}}_r^2\mathrm{}}{{\mathrm{\&sigma;}}_q^3\left(k\right)}\underset{p=1}{\overset{P-1}{\mathrm{\&Sigma;}}}{p}^3{|A_u(p.q)}^2{e}^{-\frac{{\left(p{\mathrm{\&Delta;}}_r\right)}^2}{{\mathrm{\&sigma;}}_q^2\left(k\right)}}$ 式（8）

其中，k表示迭代次数，σ_q(k)为第k次迭代时扩展向量σ(k+1)的第q个采样值；

更新迭代步长，若满足 $f\left(\mathrm{\&sigma;}(k+1)\right)-f\left(\mathrm{\&sigma;}\left(k\right)\right)\&le;{\mathrm{\&rho;}}_1\mathrm{\&mu;}\left(k\right)\&dtri;f^T\left(k\right)(\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)),$ 则步长为μ(k+1)=μ(k)/Δ；若满足 $f\left(\mathrm{\&sigma;}(k+1)\right)-f\left(\mathrm{\&sigma;}\left(k\right)\right)>{\mathrm{\&rho;}}_2\mathrm{\&mu;}\left(k\right)\&dtri;f^T\left(K\right)(\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)),$ 则步长为μ(k+1)=μ(k)Δ；

如果迭代次数达到预先的最大迭代次数λ，或者同时满足式（9）和式（10），则终止迭代，得到的扩展向量σ；

代价函数约束：f(σ(k+1))-f(σ(k))<εf(σ(k))  式（9）

扩展向量约束： $\left|\right|\mathrm{\&sigma;}(k+1)-\mathrm{\&sigma;}\left(k\right)\left|\right|<\mathrm{\&epsiv;}\sqrt{2\mathrm{\&pi;\&alpha;}/{\mathrm{\&Delta;}}_{\mathrm{\&psi;}}}$ 式（10）；

步骤1044）根据式（11）确定直角坐标系下的径向角ψ_d：

${\mathrm{\&psi;}}_d=\arctan \frac{\mathrm{\&tau;}}{\mathrm{\&theta;}}$ 式（11）

其中，τ表示离散化的时滞，θ表示离散化的频偏；

将步骤1043）得到的扩展向量σ和径向角y扩充到整个模糊域，得到式（12）和式（13）

y←[ψ(1:Q),y(2:Q)+π]  式（12）

σ←[σ(1:Q),σ(2:Q)]  式（13）

由y和σ采用插值的方式，测算出直角坐标系下径向角ψ_d对应的扩展函数σ_d，并由直角坐标系下的极径公式

得到直角坐标系下的径向高斯核函数Φ_d(m,n)，如式（14）所示，

${\mathrm{\&Phi;}}_d(m,n)=e^{-r^2/2{\mathrm{\&sigma;}}^2\left(\mathrm{\&psi;}\right)}$ 式（14）

步骤1045）将模糊函数A_d与径向高斯核函数Φ_d相乘，并利用二维傅里叶变换，得到信号的时频分布I(t,f)。

3.按照权利要求1所述的基于频散特征的低频水声脉冲信号距离的测量方法，其特征在于，所述的步骤203）中，S为采样时间间隔T_S的10至50倍。

4.按照权利要求1所述的基于频散特征的低频水声脉冲信号距离的测量方法，其特征在于，所述的步骤204）包括以下过程：当且

时，M(t,f)=0；当且

时，M(t,f)=1；其中，

表示时频分布I(t,f)中时间t位于的区间，

表示时频分布I(t,f)中频率f位于的区间，

表示时间t不位于区间

内，

表示频率f不位于区间

内。

5.按照权利要求2所述的基于互相关的水声脉冲信号双阵元定位的方法，其特征在于，所述的步骤1041）中，α的取值范围为1≤α≤5。

Images