CN103076590A - 一种基于频率预估的水声脉冲信号的定位方法 - Google Patents

Abstract

本发明公开了一种基于频率预估的水声脉冲信号的定位方法，包括：步骤10)利用单个水听器接收水声脉冲信号，对该水声脉冲信号的频率进行具有自适应径向高斯核函数的时频分布估计，首先利用单个水听器接收水声脉冲信号，对该水声脉冲信号确定相应的二维频偏—时滞域上的模糊函数，然后设置离散化频偏的总点数和离散时滞的总点数，接着将直角坐标系下的模糊函数转换为极坐标系下的模糊函数，随后采用迭代算法测算最优扩展函数，得到接收信号的时频分布，确定接收信号的中心频率以及相应的瞬时时间；步骤20）利用声场匹配对目标实施定位。该定位方法能够实现利用单个水听器对发射信号未知的水声脉冲信号的参数估计，并对该发射信号精确定位。

Description

一种基于频率预估的水声脉冲信号的定位方法

技术领域

本发明涉及针对水声信号处理技术领域的水下目标定位方法，具体来说，涉及一种基于频率预估的水声脉冲信号的定位方法。

背景技术

水声定位技术在水声技术中扮演着极其重要的角色，是水声技术中的一个重要且基本的问题，同时也是国民经济建设和国防建设的关键技术。在水声定位技术中，被动定位技术一直是水声技术中的重要研究方向。现阶段对于水下被动定位技术多围绕舰船辐射噪声而展开工作，很少有针对各种水声脉冲信号在发射信号未知方式下的被动定位。对于水声脉冲信号而言，其信号形式不同于舰船辐射的噪声，它是由人工所产生的，有比较规则的信号形式，同时在时间上具有不连续性、瞬时性，在带宽上具有窄带特性，不能通过长时间的积分来获得时间增益，这就造成了对于水声脉冲信号的被动定位需要寻求与传统的舰船辐射噪声所不同的方式和手段。

考虑到水声环境的复杂性，为了能够更准确地对水下目标实施定位，需要从声传播的角度出发，借助匹配场定位技术来确定目标声源的位置。传统意义上的匹配场定位技术，一般多采用阵列的处理方式，具有大的孔径，以获得良好的阵增益和分辨性能。但是采用多阵元的大阵列，一方面增加了系统的开销，给基阵的设计带来不便；另一方面，在实际海水中布放时会受到诸如阵倾斜以及阵元失效等问题，增加了对水下目标定位的难度。因此，只借助单个水听器来对水下目标进行定位一直被研究人员所关注，并取得了一些突破和进展。

在这方面，相关技术和研究人员多针对发射信号为已知的情形，借助宽带信号的多频点特性，利用单个水听器对目标实施定位。但是在实际中，所能获得的发射信号的先验信息有限，尤其是对先验信息不足的水声脉冲信号，若对该类信号实施被动定位，由于信号的频率，波形等信息未知，从而进一步的增加了对于此类信号定位的难度。因此，利用单水听器对水声脉冲信号进行定位之前，首先需要对水听器接收信号的参数进行预估，获得目标发射信号的频率，这是进行后续定位工作的前提。由于目标信号的频率未知，若直接采用傅里叶变换的方法在全频率段内进行频率搜索，则计算量较大，且受工作频段的限制，存在频点遗漏的情形，会造成频率估计的偏差较大。另外的一些参数估计技术，如短时傅里叶变换、魏格纳-威利分布等，从时频二维域角度出发，因考虑了水声脉冲信号短时瞬态的非平稳特性，往往具有不错的估计效果。但是这些传统的时频分析方法都是基于固定核函数的时频分析方法，只适应于某种特定的信号，存在时频分辨力低和交叉项干扰等问题，尤其是在低信噪比条件下，会直接影响到参数估计的效果。另外，在对目标信号进行匹配场定位时，需要在搜索域空间计算拷贝场信号，但在发射信号未知的条件下，不能直接求解出所需的拷贝场信号，从而造成了制约单水听器匹配场定位的瓶颈。

针对传统时频分析方法的不足，已有学者提出了一些相关的方法，如采用对核函数进行优化设计的时频分析方法，其核函数的形状能够根据所分析的信号自适应的变化，可提高对于非平稳信号参数估计的性能。但是针对发射信号未知情形下的匹配场定位问题却为数不多。虽然目前已有对此问题的研究报道，但是仍对信号的频率参数做出了已知的假设，在具体求解时，仍需要获知前向声场计算所需的频率。

发明内容

技术问题：本发明所要解决的技术问题是：提供一种基于频率预估的水声脉冲信号的定位方法，该定位方法能够实现利用单个水听器对发射信号未知的水声脉冲信号的参数估计，并对该发射信号精确定位。

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

一种基于频率预估的水声脉冲信号的定位方法，该定位方法包括以下步骤：

步骤10)利用单个水听器接收水声脉冲信号，然后对该水声脉冲信号的频率进行具有自适应径向高斯核函数的时频分布估计，包括步骤101）——步骤104）：

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

$A_d(m,n)=T_s\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{y}^{*}(k{T}_s-n{T}_s)y(k{T}_s+n{T}_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\·\mathrm{ros}\left(\mathrm{\ψ}\right)\end{array}\right.$   式（2）

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

步骤104）采用迭代算法测算最优扩展函数σ，获得接收信号的时频分布，确定接收信号的中心频率以及相应的瞬时时间；

步骤20）利用声场匹配对目标实施定位，包括步骤：201）—步骤203）：

步骤201）确定声源的搜索范围：在观测范围内，对观测范围进行网格点划分，获得网格区域(R，Z)，其中，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围；

步骤202）对步骤201）所划分的网格区域(R,Z)，使用bellhop声场传播模型，将海洋环境参数和步骤104）得到的中心频率，作为声场传播模型的输入值，通过声场传播模型测算，得到各网格区域上的声源在水听器上各声线的传播时间

以及相应的幅度

从而产生相应的拷贝信道脉冲响应

如式（16）所示，

$\hat{h}\left(u\right)=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\hat{a}}_i\mathrm{\δ}(u-{\widehat{\mathrm{\τ}}}_i)$   式（16）

其中，u表示离散时间序列索引，u=0、1、…、M-1，M表示信道脉冲响应的长度；l为在某网格点处，声源传播到达水听器的多途数；

表示u在

处的单位冲击函数；

步骤203）根据声场传播模型测算出的信道脉冲响应和水听器上的接收信号，对目标实施匹配场定位。

进一步，所述的步骤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=1}{\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{\&sigma;}}_q^3\left(k\right)}\underset{p=1}{\overset{P-1}{\mathrm{\&Sigma;}}}{p}^3{\left|A_u(p\&CenterDot;q)\right|}^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）得到的扩展向量σ和径向角ψ扩充到整个模糊域，得到式（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)。

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

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

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

表示峰值所对应的时间t和频率f。

进一步，所述的步骤203）包括如下的步骤：

步骤2031）对各网格区域上测算得到的信道脉冲响应构造如下的拷贝信道卷积矩阵，如式（17）所示，

$\hat{H}={\left[\begin{array}{cccccc}\hat{h}\left(0\right)& \hat{h}\left(0\right)& ...& \hat{h}(M-1)& ...& 0\\ 0& \hat{h}\left(0\right)& ...& \hat{h}(M-2)& ...& 0\\ .& .& & .& & .\\ .& .& ...& .& ...& .\\ .& .& & .& & .\\ 0& 0& ...& \hat{h}\left(0\right)& ...& \hat{h}(M-1)\end{array}\right]}^T$   式（17）

其中，M表示信道脉冲响应的长度，T表示转置。

步骤2032）由单个水听器接收到的水声脉冲信号和步骤2031）得到的拷贝信道卷积矩阵，采用最小二乘算法，测算出单个水听器上的拷贝场信号

其中，

表示矩阵

的逆阵，

表示

的转置；

步骤2033）建立如式（18）所示的误差代价函数L(R,Z)，测算出相应的定位模糊表面：

$L(R,Z)=1/{\left|\right|y-\hat{y}\left|\right|}^2$   式（18）

其中：

表示y与

之间误差的范数平方和，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围；

对模糊表面进行峰值搜索，根据式（19）确定目标的位置：

$({\hat{R}}_0,{\hat{Z}}_0)=\underset{R,Z}{\arg \max }L(R,Z)$   式（19）

其中：

表示定位所得到的距离估计值；

表示定位所得到的深度估计值。

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

（1）时频分辨率高，估计精度高。现有技术多采用固定核函数的时频分析方法，其时频分辨率较低。本发明采用了自适应核函数的时频分析技术来估计信号的频率参数，时频分辨率较高。本发明通过对核函数进行优化设计，使用基于信号的自适应径向高斯核函数的时频分布进行参数估计，充分利用了水声脉冲信号的短时瞬态的非平稳特性，具有较高的时频分辨率，进而提高了估计精度。

（2）定位精确。现有技术多针对多阵元的处理方法以及发射信号已知的情形，而本发明针对发射信号的先验信息不足情形下的单一接收水听器的定位问题。本发明通过采用高分辨率的时频分析技术来获得接收信号的频率参数；采用时域最小二乘的方法测算拷贝场信号，无需已知发射信号的波形信息，体现了算法本身对于信号参数较低的依赖性，具有高的定位精度，同时采用单个水听器避免了阵列处理时的诸多不便。

附图说明

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

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

图3是本发明实施例接收信号的频率曲线图。

图4是本发明实施例的定位模糊表面示意图。

具体实施方式

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

本发明的一种基于频率预估的水声脉冲信号的定位方法，包括步骤10）——步骤30）：

步骤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}^{*}(k{T}_s-n{T}_s)y(k{T}_s+n{T}_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;\mathrm{ros}\left(\mathrm{\&psi;}\right)\end{array}\right.$   式（2）

其中，p表示极坐标系下的离散化频偏索引，q表示极坐标系下的离散化时滞索引，X表示转换后的频偏坐标，Y表示转换后的时滞坐标。二维插值的方法可以使用matlab软件提供的插值函数interp2(·)来实现，并将插值后位于最大频偏τ_max和最大时滞θ_max之外的值设为0。

步骤104）采用迭代算法测算最优扩展函数σ，获得接收信号的时频分布，确定接收信号的中心频率以及相应的瞬时时间。步骤104）具体包括步骤1041）——步骤1046）。

步骤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=1}{\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(p{\mathrm{\&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(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{\&sigma;}}_q^3\left(k\right)}\underset{p=1}{\overset{P-1}{\mathrm{\&Sigma;}}}{p}^3{\left|A_u(p\&CenterDot;q)\right|}^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）得到的扩展向量σ和径向角ψ扩充到整个模糊域，得到式（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)。

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

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

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

表示峰值所对应的时间t和频率f。

步骤20）利用声场匹配对目标实施定位。步骤20）具体包括步骤：201）—步骤203）。

步骤201）确定声源的搜索范围：在观测范围内，对观测范围进行网格点划分，获得网格区域(R，Z)，其中，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围。

步骤202）对步骤201）所划分的网格区域(R,Z)，使用bellhop声场传播模型，将海洋环境参数和步骤104）得到的中心频率，作为声场传播模型的输入值，通过声场传播模型测算，得到各网格区域上的声源在水听器上各声线的传播时间

以及相应的幅度从而产生相应的拷贝信道脉冲响应

如式（16）所示，

$\hat{h}\left(u\right)=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}{\hat{a}}_i\mathrm{\&delta;}(u-{\widehat{\mathrm{\&tau;}}}_i)$   式（16）

其中，u表示离散时间序列索引，u=0、1、…、M-1，M表示信道脉冲响应的长度；l为在某网格点处，声源传播到达水听器的多途数；

表示u在

处的单位冲击函数。

在步骤202）中，海洋环境参数包括声速梯度分布、海水深度、海水密度、海底密度、海底的衰减系数。

步骤203）根据声场传播模型测算出的信道脉冲响应和水听器上的接收信号，对目标实施匹配场定位。步骤203）具体包括步骤2031）——步骤2033）。

步骤2031）对各网格区域上测算得到的信道脉冲响应

构造如下的拷贝信道卷积矩阵，如式（17）所示，

$\hat{H}={\left[\begin{array}{cccccc}\hat{h}\left(0\right)& \hat{h}\left(0\right)& ...& \hat{h}(M-1)& ...& 0\\ 0& \hat{h}\left(0\right)& ...& \hat{h}(M-2)& ...& 0\\ .& .& & .& & .\\ .& .& ...& .& ...& .\\ .& .& & .& & .\\ 0& 0& ...& \hat{h}\left(0\right)& ...& \hat{h}(M-1)\end{array}\right]}^T$   式（17）

其中，M表示信道脉冲响应的长度，T表示转置。

步骤2032）由单个水听器接收到的水声脉冲信号和步骤2031）得到的拷贝信道卷积矩阵，采用最小二乘算法，测算出单个水听器上的拷贝场信号

其中，表示矩阵

的逆阵，

表示

的转置；

步骤2033）建立如式（18）所示的误差代价函数L(R,Z)，测算出相应的定位模糊表面：

$L(R,Z)=1/{\left|\right|y-\hat{y}\left|\right|}^2$   式（18）

其中：

表示y与

之间误差的范数平方和，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围；

对模糊表面进行峰值搜索，根据式（19）确定目标的位置：

$({\hat{R}}_0,{\hat{Z}}_0)=\underset{R,Z}{\arg \max }L(R,Z)$   式（19）

其中：

表示定位所得到的距离估计值；

表示定位所得到的深度估计值。

本发明的基于频率预估的水声脉冲信号的定位方法，对基于信号的自适应径向高斯核函数的时频分布进行参数估计，获得目标信号的频率；采用时域最小二乘的方法估计搜索区域上的拷贝场信号，建立接收信号与拷贝场信号之间的误差函数、网格区域的搜索以及定位结果显示。本发明的技术方案能够实现只借助单一接收水听器即可完成对发射信号未知的水声脉冲信号的定位。

下面例举一实施例。

本实施例的声学环境如图1所示。声学环境的各项参数为：S代表深度为60m的目标信号；R代表深度为80m的单一接收水听器，目标信号与接收水听器之间的距离为5Km；水层声速分布为负跃变层分布，水层深度为110m，底部声速为1700m/s，密度为1.9g/cm3，衰减系数为0.5dB/λ。发射信号为频带范围为100-200Hz的线性调频信号。

图2是本实施例中接收信号的时频表示。图2的横坐标表示时间，单位为秒(s)，纵坐标表示频率，单位赫兹(Hz)。图2采用matlab软件绘制。从图2可以看出：采用自适应径向高斯核函数的时频分布可较好地表征水声脉冲信号，可反映出水声脉冲信号的时频特征。图3为本实施例中对应的频率估计结果。峰值的位置显示了信号的中心频率。图3的横坐标表示频率，单位赫兹(Hz)，纵坐标表示归一化的幅度。按照本发明的方法定位水声脉冲信号，具体包括如下步骤：

（1）测算单一接收水听器接收信号的模糊函数并将其转换到极坐标系下；

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

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

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

（5）对所得的时频分布结果进行峰值搜索，获得关于目标信号的中心频率。

本实施例对目标信号的中心频率的测算结果为149.6Hz，接近目标信号的真实频率150Hz。由此可知，所估计得到的信号频率较为精确。另外，从图2和3中也可以看出：采用自适应径向高斯核函数的时频分布可以较好地表征水声脉冲信号，所得到的频率波形较为尖锐。

图4是本实施例的定位模糊表面图。图4采用matlab软件绘制。图4的横坐标表示距离，单位千米(Km)，纵坐标表示深度，单位米（m）。完成定位模糊表面需包括以下步骤：

(1)对搜索区域进行网格域划分，同时配置好相关声场计算所需要的环境信息。在本实施例中，网格点按照距离为2Km-7Km，步距为100m，深度为5m-105m步距为2.5m的范围进行划分；

(2)选择bellhop声传播模型，将配置好的环境文件以及估计所得到的频率值代入声传播模型中，在所划分的网格区域上进行前向声场的计算，获得每个网格点上的测试声源在接收机处的拷贝信道脉冲响应，产生拷贝信道矩阵；

(3)通过单个水听器所接收的信号和测算出的拷贝信道矩阵，采用时域最小二乘的方法获得搜索网格域上发射信号的估计值，由此产生拷贝场信号；

(4)构造接收信号与拷贝场信号之间的误差函数，测得定位模糊表面，进行网格区域搜索，获得模糊表面的峰值，从而确定目标源的位置。

从图4的最终的结果可以看出，本发明的定位模糊表面的峰值较为清晰，定位结果正确地反映了目标所在的位置。

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

Claims (5)

1.一种基于频率预估的水声脉冲信号的定位方法，其特征在于：该定位方法包括以下步骤：

步骤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)采用迭代算法测算最优扩展函数σ，获得接收信号的时频分布，确定接收信号的中心频率以及相应的瞬时时间；

步骤20)利用声场匹配对目标实施定位，包括步骤：201)-步骤203)：

步骤201)确定声源的搜索范围：在观测范围内，对观测范围进行网格点划分，获得网格区域(R，Z)，其中，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围；

步骤202)对步骤201)所划分的网格区域(R，Z)，使用bellhop声场传播模型，将海洋环境参数和步骤104)得到的中心频率，作为声场传播模型的输入值，通过声场传播模型测算，得到各网格区域上的声源在水听器上各声线的传播时间

以及相应的幅度

从而产生相应的拷贝信道脉冲响应如式(16)所示，

$\hat{h}\left(u\right)=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}{\hat{a}}_i\mathrm{\&delta;}(u-{\widehat{\mathrm{\&tau;}}}_i)$ 式(16)

其中，u表示离散时间序列索引，u＝0、1、...、M-1，M表示信道脉冲响应的长度；l为在某网格点处，声源传播到达水听器的多途数；

表示u在

处的单位冲击函数；

步骤203)根据声场传播模型测算出的信道脉冲响应和水听器上的接收信号，对目标实施匹配场定位。

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)},\&CenterDot;\&CenterDot;\&CenterDot;,\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)

扩展向量约束： $(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)。

步骤1046)通过对时频分布I(t，f)进行峰值搜索，确定接收信号的中心频率以及相应的瞬时时间，如式(15)所示，

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

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

3.按照权利要求1所述的基于频率预估的水声脉冲信号的定位方法，其特征在于，所述的步骤203)包括如下的步骤：

步骤2031)对各网格区域上测算得到的信道脉冲响应构造如下的拷贝信道卷积矩阵，如式(17)所示， $\hat{H}={\left[\begin{array}{cccccc}\hat{h}\left(0\right)& \hat{h}\left(1\right)& ...& \hat{h}(M-1)& ...& 0\\ 0& \hat{h}\left(0\right)& & \hat{h}(M-2)& ...& 0\\ .& .& & .& & .\\ .& .& ...& .& ...& .\\ .& .& & .& & .\\ 0& 0& ...& \hat{h}\left(0\right)& ...& \hat{h}(M-1)\end{array}\right]}^T$ 式(17)

其中，M表示信道脉冲响应的长度，T表示转置。

步骤2032)由单个水听器接收到的水声脉冲信号和步骤2031)得到的拷贝信道卷积矩阵，采用最小二乘算法，测算出单个水听器上的拷贝场信号

其中，

表示矩阵

的逆阵，

表示

的转置；

步骤2033)建立如式(18)所示的误差代价函数L(R，Z)，测算出相应的定位模糊表面：

$L(R,Z)=1/{\left|\right|y-\hat{y}\left|\right|}^2$ 式(18)

其中：

表示y与

之间误差的范数平方和，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围；

对模糊表面进行峰值搜索，根据式(19)确定目标的位置：

$({\hat{R}}_0,{\hat{Z}}_0)=\underset{R,Z}{\arg \max }L(R,Z)$ 式(19)

其中：

表示定位所得到的距离估计值；

表示定位所得到的深度估计值。

4.按照权利要求2所述的基于频率预估的水声脉冲信号的定位方法，其特征在于，所述的步骤1041)中，α的取值范围为1≤α≤5。

5.按照权利要求1所述的基于频率预估的水声脉冲信号的定位方法，其特征在于，所述的步骤202)中，海洋环境参数包括声速梯度分布、海水深度、海水密度、海底密度、海底的衰减系数。

Images