CN103076594B - 一种基于互相关的水声脉冲信号双阵元定位的方法 - Google Patents

Abstract

本发明公开了一种基于互相关的水声脉冲信号双阵元定位的方法：步骤10）对水声脉冲信号的频率进行具有自适应径向高斯核函数的时频分布估计，包括如下的步骤：步骤101）确定二维频偏—时滞域上的模糊函数；步骤102）设置离散频偏的总点数和离散时滞的总点数，步骤103）将直角坐标系下的模糊函数转化为极坐标系下的模糊函数：步骤14）确定接收信号的中心频率以及相应的瞬时时间；步骤20）对双阵元接收的水声脉冲信号实施定位，包括步骤201）确定声源的搜索范围；步骤202）确定相应的拷贝信道脉冲响应；步骤203）对目标实施匹配场定位。该方法利用两个垂直接收水听器对发射信号未知的水声脉冲信号实现精确定位。

Description

一种基于互相关的水声脉冲信号双阵元定位的方法

技术领域

本发明涉及针对水声信号处理技术领域的水下目标定位方法，具体来说，涉及一种基于互相关的水声脉冲信号双阵元定位的方法。

背景技术

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

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

使用较少的阵元个数进行目标位置估计的一个难点在于空间信息的缺乏。多数研究人员借助宽带信号的多频点特性，从假设发射信号为已知的情形出发，采用“频点换孔径”的思想来对目标信号实施定位。但在被动定位技术中，由于所获得的发射信号的先验信息有限，往往需要在发射信号未知的情形下进行。由于信号的波形、频率等信息未知，从而进一步增加了定位的难度。另外，传统匹配场处理技术当阵元个数变少时，其定位结果输出的模糊表面旁瓣较高，无法对目标进行精确定位，也成为了困扰此项技术的一个瓶颈。

因此，在对发射信号未知的水声脉冲信号进行定位之前，首先需要获得接收信号的相关参数。对于发射信号未知的水声脉冲信号而言，准确获知目标信号的频率至关重要，这也是后续定位工作的前提。由于目标信号的频率未知，若直接采用傅里叶变换的方法在全频率段内进行频率搜索，则计算量较大，且受工作频段的限制，存在频点遗漏的情形，会造成频率估计的偏差较大。另外的一些参数估计技术，如短时傅里叶变换、魏格纳-威利分布等，从时频二维域角度出发，因考虑了水声脉冲信号短时瞬态的非平稳特性，往往具有不错的估计效果。但是这些传统的时频分析方法都是基于固定核函数的时频分析方法，只适应于某种特定的信号，存在时频分辨力低和交叉项干扰等问题，尤其是在低信噪比条件下，会直接影响到参数估计的效果。针对传统时频分析方法的不足，已有学者提出了一些相关的方法，如采用对核函数进行优化设计的时频分析方法，其核函数的形状能够根据所分析的信号自适应的变化，可提高对于非平稳信号参数估计的性能。但是针对发射信号未知情形下的匹配场定位问题却为数不多。虽然目前已有对此问题的研究报道，但是仍对信号的频率参数做出了已知的假设，在具体求解时，仍需要获知前向声场计算所需的频率。

发明内容

技术问题：本发明所要解决的技术问题是：提供一种基于互相关的水声脉冲信号双阵元定位的方法，该方法能够实现利用两个垂直接收水听器对发射信号未知的水声脉冲信号的精确定位。

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

一种基于互相关的水声脉冲信号双阵元定位的方法，该方法包括以下步骤：

步骤10)双阵元水听器接收水声脉冲信号，然后对该水声脉冲信号的频率进行具有自适应径向高斯核函数的时频分布估计，包括如下的步骤：

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

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

其中，N表示接收信号的点长；m表示直角坐标系下的离散化的频偏索引；n表示直角坐标系下的离散化的时滞索引；T_s表示采样时间间隔；y_j表示双阵元水听器所接收到的水声脉冲信号，j＝1和2；表示y_j的共轭；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)采用迭代算法测算最优扩展函数σ，获得接收信号的时频分布，确定出接收信号的中心频率以及相应的瞬时时间；

步骤20)对双阵元水听器接收的水声脉冲信号实施定位，包括如下的步骤：

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

步骤202)对步骤201)所划分的网格区域(R，Z)，使用bellhop声场传播模型，将海洋环境参数和步骤1046)得到的中心频率，作为声场传播模型的输入值，通过声场传播模型测算，得到各网格区域上的声源在每个水听器上各声线的传播时间以及相应的幅度从而产生相应的拷贝信道脉冲响应如式(16)所示，

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

其中，u表示离散时间序列索引，u＝0、1、...、M-1，j表示接收机，j＝1和2，1为在某网格点处，声源传播到达接收机j的多途数；表示u在处的单位冲击函数；

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

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

步骤1041)设置常量参数及初值：收敛因子ε、最大迭代次数λ、初始步长μ(0)、初始扩展向量σ(0)、步长控制因子Δ、常数ρ₁和ρ₂，0＜ρ₁＜ρ₂＜1，Δ＞0，其中，α表示核函数体积控制参数，Q表示离散时滞的总点数，Δ_ψ表示径向角分辨率，[1，...，1]表示1构成的向量，T表示转置；

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

$f\left(\mathrm{\σ}\left(0\right)\right)=\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}\underset{p=2}{\overset{p}{\mathrm{\Σ}}}{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\Δ}}_r\right)}^2}{{\mathrm{\σ}}_q^2\left(0\right)}}$ 式(3)

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

$\▿f\left(0\right)={[\frac{\∂f}{\∂{\mathrm{\σ}}_1\left(0\right)},\·\·\·,\frac{\∂f}{\∂{\mathrm{\σ}}_Q\left(0\right)}]}^T$ 式(4)

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

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

$\mathrm{\σ}(k+1)=\mathrm{\σ}\left(k\right)+\mathrm{\μ}\left(k\right)\▿f\left(k\right)$ 式(5)

$\mathrm{\σ}(k+1)\←\frac{\mathrm{\σ}(k+1)}{\left|\right|\mathrm{\σ}(k+1)\left|\right|}\sqrt{\frac{2\mathrm{\π\α}}{{\mathrm{\Δ}}_{\mathrm{\ψ}}}}$ 式(6)

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

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

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

$\▿f\left(k\right)=\frac{{2\mathrm{\Δ}}_r^2}{{\mathrm{\σ}}_q^3\left(k\right)}\underset{p=1}{\overset{P-1}{\mathrm{\Σ}}}{p}^3{\left|A_u(p\·q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\Δ}}_r\right)}^2}{{\mathrm{\σ}}_q^2\left(k\right)}}$ 式(8)

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

更新迭代步长，若满足 $f\left(\mathrm{\σ}(k+1)\right)-f\left(\mathrm{\σ}\left(k\right)\right)\≤{\mathrm{\ρ}}_1\mathrm{\μ}\left(k\right)\▿f^T\left(k\right)(\mathrm{\σ}(k+1)-\mathrm{\σ}\left(k\right)),$ 则步长为μ(k+1)＝μ(k)/Δ；若满足 $f\left(\mathrm{\σ}(k+1)\right)-f\left(\mathrm{\σ}\left(k\right)\right)>{\mathrm{\ρ}}_2\mathrm{\μ}\left(k\right)\▿f^T\left(k\right)(\mathrm{\σ}(k+1)-\mathrm{\σ}\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 }(t,f)$ 式(15)

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

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

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

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

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

步骤2032)由双阵元水听器接收的水声脉冲信号和步骤2031)得到的拷贝信道卷积矩阵，采用最小二乘算法，测算出每个接收阵元上的拷贝场信号其中，表示矩阵的逆阵，表示的转置，y_j表示双阵元水听器所接收到的水声脉冲信号；

步骤2033)根据式(18)测算双阵元水听器接收的水声脉冲信号的互相关函数ρ₁₂，根据式(19)测算拷贝场信号的互相关函数

式(18)

式(19)

其中，□表示互相关运算，y₁表示双阵元水听器中的一个水听器接收的水声脉冲信号，y₂表示双阵元水听器中的另一个水听器接收的水声脉冲信号，表示双阵元水听器中的一个拷贝场信号，表示双阵元水听器中的另一个拷贝场信号，v表示接收信号的离散时间索引，m＝-N+1，…，N-1；

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

$L(R,Z)=1/{\left|\right|{\mathrm{\&rho;}}_{12}-{\widehat{\mathrm{\&rho;}}}_{12}\left|\right|}^2$ 式(20)

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

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

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

其中：表示定位所得到的距离估计值；表示定位所得到的深度估计值。

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

(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}_j^{*}({\mathrm{kT}}_s-{\mathrm{nT}}_s)y_j({\mathrm{kT}}_s+{\mathrm{nT}}_s)e^{i2\mathrm{\&pi;mk}/N}$ 式(1)

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

步骤102)首先设置离散频偏的总点数P和离散时滞的总点数Q，其中，Q＝N；然后从产生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表示转换后的时滞坐标。二维插值的方法可以使用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=2}{\overset{p}{\mathrm{\&Sigma;}}}{\left|A_u(p,q)\right|}^2{e}^{-\frac{{\left({\mathrm{p\&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表示向量的转置。

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

$\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\&CenterDot;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)。

步骤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声场传播模型，将海洋环境参数和步骤1046)得到的中心频率，作为声场传播模型的输入值，通过声场传播模型测算，得到各网格区域上的声源在每个水听器上各声线的传播时间以及相应的幅度从而产生相应的拷贝信道脉冲响应如式(16)所示，

${\hat{h}}_j\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，j表示接收机，j＝1和2，1为在某网格点处，声源传播到达接收机j的多途数；表示u在处的单位冲击函数。海洋环境参数包括声速梯度分布、海水深度、海水密度、海底密度、海底的衰减系数等。

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

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

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

步骤2032)由双阵元水听器接收的水声脉冲信号和步骤2031)得到的拷贝信道卷积矩阵，采用最小二乘算法，测算出每个接收阵元上的拷贝场信号其中，表示矩阵的逆阵，表示的转置，y_j表示双阵元水听器所接收到的水声脉冲信号。

步骤2033)根据式(18)测算双阵元水听器接收的水声脉冲信号的互相关函数ρ₁₂，根据式(19)测算拷贝场信号的互相关函数

式(18)

式(19)

其中，□表示互相关运算，y₁表示双阵元水听器中的一个水听器接收的水声脉冲信号，y₂表示双阵元水听器中的另一个水听器接收的水声脉冲信号，表示双阵元水听器中的一个拷贝场信号，表示双阵元水听器中的另一个拷贝场信号，v表示接收信号的离散时间索引，m＝-N+1，…，N-1。

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

$L(R,Z)=1/{\left|\right|{\mathrm{\&rho;}}_{12}-{\widehat{\mathrm{\&rho;}}}_{12}\left|\right|}^2$ 式(20)

其中：表示ρ₁₂与之间误差的范数平方和，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围。

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

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

其中：表示定位所得到的距离估计值；表示定位所得到的深度估计值。

本发明技术方案能够实现利用两垂直双阵元完成对发射信号未知的水声脉冲信号的定位。

下面例举一实施例。

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

图2是本发明的双阵元水听器的接收信号的时频表示。图2采用matlab软件绘制。图2的横坐标表示时间，单位为秒(s)，纵坐标表示频率，单位赫兹(Hz)。按照本发明的方法对水声脉冲信号双阵元定位，具体包括如下步骤：

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

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

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

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

(5)对所得的时频分布结果进行峰值搜索，获得关于目标信号的中心频率。本实施例所给出的结果为149.6Hz，接近目标信号的真实频率150Hz。

由此可知，所估计得到的信号频率较为精确。另外，从图2也可以看出：采用自适应径向高斯核函数的时频分布可以较好表征水声脉冲信号。

图3是本发明实施例的定位模糊表面。图3采用matlab软件绘制。图3中，横坐标表示距离，单位千米(km)，纵坐标表示深度，单位米(m)。完成定位过程所需的步骤包括：

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

(2)选择bellhop射线模型，将配置好的海洋环境参数以及估计所得到的频率值代入声传播模型中，在所划分的网格区域上进行前向声场的计算，获得每个网格点上的测试声源在双阵元处的拷贝信道脉冲响应，并产生相应的拷贝信道卷积矩阵。

(3)通过双阵元上的接收的信号和计算出的拷贝信道卷积矩阵，采用时域最小二乘的方法获得搜索网格域上上发射信号的估计值，并由此产生双阵元上的拷贝场信号。

(4)测算双阵元接收信号的互相关函数与拷贝场信号的互相关函数、

(5)由接收信号的互相关和拷贝场信号的互相关构造误差代价函数，确定相应的定位模糊表面，进行网格区域搜索，获得模糊表面的峰值，从而确定目标源的位置。

由估计的结果可以看出，本发明的定位模糊表面的峰值较为清晰且旁瓣较低。从图3可以看出：最终的定位结果，也即声源的位置。

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

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}_j^{*}(k{T}_s-n{T}_s)y_j(k{T}_s+n{T}_s)e^{i2\mathrm{\&pi;mk}/N}$   式(1)

其中，N表示接收信号的点长；m表示直角坐标系下的离散化的频偏索引；n表示直角坐标系下的离散化的时滞索引；T_S表示采样时间间隔；y_j表示双阵元水听器所接收到的水声脉冲信号，j＝1和2；表示y_j的共轭；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)确定声源的搜索范围：在观测范围内，对观测范围进行网格点划分，获得网格区域(R,Z)，其中，R表示搜索网格区域上的距离范围，Z表示搜索网格区域上的深度范围；

步骤202)对步骤201)所划分的网格区域(R,Z)，使用bellhop声场传播模型，将海洋环境参数和步骤1046)得到的中心频率，作为声场传播模型的输入值，通过声场传播模型测算，得到各网格区域上的声源在每个水听器上各声线的传播时间以及相应的幅度从而产生相应的拷贝信道脉冲响应如式(16)所示，

${\hat{h}}_j\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，j表示接收机，j＝1和2，l为在某网格点处，声源传播到达接收机j的多途数；表示u在处的单位冲击函数；

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

2.按照权利要求1所述的基于互相关的水声脉冲信号双阵元定位的方法，其特征在于，所述的步骤104)包括以下步骤：

步骤1041)设置常量参数及初值：收敛因子ε、最大迭代次数λ、初始步长μ(0)、初始扩展向量σ(0)、步长控制因子Δ、常数ρ₁和ρ₂，0＜ρ₁＜ρ₂＜1，Δ>0，其中，α表示核函数体积控制参数，Q表示离散时滞的总点数，Δ_ψ表示径向角分辨率，[1,…,1]表示1构成的向量，T表示转置；

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

$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)所示

σ(k+1)＝σ(k)+μ(k)▽f(k)   式(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次的迭代步长，▽f(k)表示第k次迭代的梯度向量，←表示更新，||σ(k+1)||表示第k+1次迭代的扩展向量σ(k+1)的范数；

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

$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(σ(k+1))-f(σ(k))≤ρ₁μ(k)▽f^T(k)(σ(k+1)-σ(k))，则步长为μ(k+1)＝μ(k)/Δ；若满足f(σ(k+1))-f(σ(k))＞ρ₂μ(k)▽f^T(k)(σ(k+1)-σ(k))，则步长为μ(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。

3.按照权利要求1所述的基于互相关的水声脉冲信号双阵元定位的方法，其特征在于，

所述的步骤203)包括如下步骤：

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

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

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

步骤2032)由双阵元水听器接收的水声脉冲信号和步骤2031)得到的拷贝信道卷积矩阵，采用最小二乘算法，测算出每个接收阵元上的拷贝场信号 ${\hat{y}}_j={\hat{H}}_j{\left({\hat{H}}_j^T{\hat{H}}_j\right)}^{-1}{\hat{H}}_j^T{y}_j,$ 其中，表示矩阵的逆阵，表示的转置，y_j表示双阵元水听器所接收到的水声脉冲信号；

步骤2033)根据式(18)测算双阵元水听器接收的水声脉冲信号的互相关函数ρ₁₂，根据式(19)测算拷贝场信号的互相关函数

  式(18)

  式(19)

其中，⊙表示互相关运算，y₁表示双阵元水听器中的一个水听器接收的水声脉冲信号，y₂表示双阵元水听器中的另一个水听器接收的水声脉冲信号，表示双阵元水听器中的一个拷贝场信号，表示双阵元水听器中的另一个拷贝场信号，v表示接收信号的离散时间索引，m＝-N+1,…,N-1；

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

$L(R,Z)=1/{\left|\right|{\mathrm{\&rho;}}_{12}-{\widehat{\mathrm{\&rho;}}}_{12}\left|\right|}^2$   式(20)

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

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

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

其中：表示定位所得到的距离估计值；表示定位所得到的深度估计值。

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

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