0 引　言

目标被动定位在水声领域具有重要的研究意义。常见的定位方法主要有三子阵测距^[1]、目标运动分析^[2]、匹配场^[3]、匹配模^[4]、多普勒分析^[5]等方法。三子阵测距利用目标到各子阵之间时延差的不同实现，当距离较远时，由于曲率减小导致方法性能退化。目标运动分析利用水平阵估计目标来波方向解算出目标参数，相关研究层出不穷^[6-9], 不过需要接收平台保持机动。匹配类的方法需要环境参数先验信息，当先验信息与实际环境参数存在失配，估计结果会产生明显偏差。多普勒分析适用频率较高且运动速度较快目标，水下运动目标往往线谱频率偏低且速度较慢。而基于声场特性定位方法从海洋声传播特性出发可有效避免上述方法面临的问题。已有方法包括阵不变量^[10]、群速度时延差^[11]、互谱估计^[12-13]、波导不变量^[14-15]等。阵不变量和群速度时延差需假设目标发射宽带脉冲信号，基于简正波到达结构实现目标距离估计。互谱估计通常需要对相对速度随时间的变换关系进行长时段拟合实现目标距离估计。

波导不变量利用运动目标辐射噪声的LOFAR谱中连续谱条纹特征进行距离估计。仅观测一段时长内的谱图中条纹特征，根据条纹特征实现距离估计。根据是否存在CPA可分为2类：无CPA时即目标沿水听器径向匀速运动条件下，具有代表性的有Cockrell等提出的划分子窗稳健波导不变量测距方法^[16]；存在CPA情况的代表有Turgut等^[17]提出的针对时间方位历程图谱和LOFAR谱进行Hough变换处理获得目标航向角和距速比（最近通过距离与目标速度之比）参数。通过对阵元上LOFAR谱条纹图进行Hough变换可获得波导不变量值 $ \beta $ ，但无法估计出目标的距离信息。余赟等在Turgut的基础上，提出利用双阵元测距方法，实现目标最近通过距离和最近通过时间的估计及被动测距^[18-19]。然而，由于该测距方法是利用波导不变量关系式 $ r\left( t \right) = \beta \omega {{{\text{d}}r\left( t \right)} \mathord{\left/ {\vphantom {{{\text{d}}r\left( t \right)} {{\text{d}}\omega }}} \right. } {{\text{d}}\omega }} $ ，存在条纹斜率求导项。当声源临近最近通过点时，距离与时间呈非线性变化，该式的适用性受限。

针对上述问题，本文在双阵元测距算法的基础上，提出基于简正波水平波数差与波导不变量关系式的被动测距方法。利用运动目标声源通过CPA位置在水平阵两端水听器产生的宽带连续谱干涉条纹，获得最近通过时间 $ {t_{{\text{cpa}}1}} $ 和 $ {t_{{\text{cpa2}}}} $ ，再利用几何关系和简正波水平波数差的频散表达式，估计出阵元1和阵元2最近通过距离 $ {r_{{\text{cpa}}1}} $ 和 $ {r_{{\text{cpa2}}}} $ 以及目标运动速度 $ {v_0} $ ，最终实现运动目标相对接收阵几何中心的距离估计。本文方法与文献[19]相比不同之处在于推导过程中使用水平波数差频散关系式替代了条纹斜率求导项，可应用于临近CPA位置处的时间-频率域具有一定曲率干涉条纹。

1 基础理论 1.1 目标运动方程及干涉条纹方程

匀速直线运动目标相对水听器阵中心的距离 $ r $ 随时间 $ t $ 的变化关系可表示为：

$ r\left( t \right) = \sqrt {v_0^2{{\left( {t - {t_{{\text{cpa}}}}} \right)}^2} + r_{{\text{cpa}}}^2}。$

(1)

式中： $ {r_{{\text{cpa}}}} $ 为阵中心与最近通过点的距离； $ {t_{{\text{cpa}}}} $ 为通过CPA位置的时刻； $ {v_0} $ 表示目标运动速度。目标运动轨迹满足如下关系：

$ \tan \left( {\theta - \varphi } \right) = \frac{{{r_{{\text{cpa}}}}}}{{{v_0}\left( {t - {t_{{\text{cpa}}}}} \right)}} 。$

(2)

式中： $ \theta $ 为相对于端射方向的方位角， $ \varphi $ 为目标航向角。式（2）可进一步写为：

$ \theta \left( t \right) = {\tan ^{ - 1}}\left( {\frac{{\left( {t - {t_{{\text{cpa}}}}} \right)\sin \varphi + \frac{{{r_{{\text{cpa}}}}}}{{{v_0}}}\cos \varphi }}{{\left( {t - {t_{{\text{cpa}}}}} \right)\cos \varphi - \frac{{{r_{{\text{cpa}}}}}}{{{v_0}}}\sin \varphi }}} \right)。$

(3)

目标与水平阵布放示意图如图1所示。在航向角 $ \varphi $ 为45°， $ {t_{{\text{cpa}}}} $ 为250 s， $ {r_{{\text{cpa}}}} $ 为1250 m， $ {v_0} $ 为5 m/s的理想条件下，运动目标时间方位谱图如图2所示。

可以看出，目标方位角在0 s时刻对应 $ \theta $ 为180°，随着目标运动，方位角度减小，在500 s时刻目标处于水平阵正横方向为90°，随后持续减小。可看出，方位角度变化率随着时间 $ t $ 先增加随后减小，即从数学角度看变化曲线先是凸函数而后是凹函数，其驻点处时刻恰好对应最近通过时间。

对目标方位时间谱图做Hough变换后的能量分布伪彩图，如图3所示。可以看出，能量峰值对应的航向角为45°，对应最近通过距离与速度比值为250，与实际值一致。

另外，宽带运动目标声源的距离-频率干涉声强可以表示为^[20]：

$ \begin{split} & I\left( {r,\omega } \right) \propto p\left( {r,\omega } \right) \cdot {p^ * }\left( {r,\omega } \right)= \\ & \left( {\mathop \sum \limits_{m{\text{ = 1}}}^M B_m^2 + 2\mathop \sum \limits_{m = 1}^M \mathop \sum \limits_{n = 1,n \ne m}^M {B_m}B_n^ * \cos \left( {{{\Delta }}{k_{mn}}r} \right)} \right) 。\end{split} $

(4)

根据波导不变量定义，干涉条纹的斜率与距离、频率满足如下关系：

$ \frac{{{\text{d}}\omega }}{{{\text{d}}r}} = \frac{\omega }{r}\beta ，$

(5)

由于 $ \omega \left( r \right) $ 和 $ r\left( t \right) $ 都可导，时间-频率域干涉条纹的斜率如下式：

$ \frac{{{\text{d}}\omega }}{{{\text{d}}t}} = \frac{{{\text{d}}\omega }}{{{\text{d}}r}}\frac{{{\text{d}}r}}{{{\text{d}}t}}。$

(6)

由式（1）可知：

$ \frac{{{\text{d}}r}}{{{\text{d}}t}} = \frac{{v_0^2\left( {t - {t_{{\text{cpa}}}}} \right)}}{{\sqrt {v_0^2{{\left( {t - {t_{{\text{cpa}}}}} \right)}^2} + r_{{\text{cpa}}}^2} }} ，$

(7)

将式（5）、式（7）代入式（6），可得：

$ \frac{{{\text{d}}\omega }}{\omega } = \frac{{\beta v_0^2\left( {t - {t_{{\text{cpa}}}}} \right)}}{{v_0^2{{\left( {t - {t_{{\text{cpa}}}}} \right)}^2} + r_{{\text{cpa}}}^2}}{\text{d}}t 。$

(8)

对上式两侧同时进行积分，整理可得频率关于时间的变化关系：

$ \omega = {\omega _{{\text{cpa}}}}{\left[ {1 + {{\left( {\frac{{{v_0}}}{{{r_{{\text{cpa}}}}}}} \right)}^2}{{\left( {t - {t_{{\text{cpa}}}}} \right)}^2}} \right]^{{\beta \mathord{\left/ {\vphantom {\beta 2}} \right. } 2}}}。$

(9)

式中， $ {\omega _{{\text{cpa}}}} $ 为时间-频率域条纹的最近通过点时刻对应频率，对应时间频率变化关系曲线的极值点。不同 $ {\omega _{{\text{cpa}}}} $ 值下的时间频率关系曲线对应LOFAR图上不同的干涉条纹。

1.2 水平阵测距方法

水平阵两端处的2个阵元之间的距离为 $ d $ ，目标声源航向角为 $ \varphi $ ，目标相对接收阵元1的最近通过时刻为 $ {t_{{\text{cpa}}1}} $ ，目标相对接收阵元2的最近通过时刻为 $ {t_{{\text{cpa}}2}} $ 。根据几何关系，目标在时刻 $ {t_{{\text{cpa}}2}} $ 相对接收阵元1的距离可表示为：

$ {r_1}\left( {{t_{{\text{cpa}}2}}} \right) = \sqrt {r_{{\text{cpa}}1}^2 + {{\left( {d\cos \varphi } \right)}^2}}。$

(10)

式中， $ {r_{{\text{cpa}}1}} $ 表示目标相对阵元1的最近通过距离。

简正波的水平波数差随频率的变化可由波导不变量近似表示为^[20]：

$ \Delta {k_{mn}} = {\gamma _{mn}}{\omega ^{ - {1 \mathord{\left/ {\vphantom {1 \beta }} \right. } \beta }}}，$

(11)

式中， $ {\gamma _{mn}} $ 表示第 $ m $ 阶与第 $ n $ 阶简正波的频散参数。声强相干项可看作由频率和距离共同影响，它们之间的关系满足

$ \omega _{1j}^{ - {1 \mathord{\left/ {\vphantom {1 \beta }} \right. } \beta }}\left( {{t_{{\text{cpa}}1}}} \right){r_{{\text{cpa}}1}} = \omega _{1j}^{ - {1 \mathord{\left/ {\vphantom {1 \beta }} \right. } \beta }}\left( {{t_{{\text{cpa}}2}}} \right){r_1}\left( {{t_{{\text{cpa}}2}}} \right)，$

(12)

式中， $ {\omega _{1j}}\left( {{t_{{\text{cpa}}1}}} \right) $ 为阵元1接收信号第 $ j $ 条干涉条纹在最近通过时刻 $ {t_{{\text{cpa}}1}} $ 处的频率。 $ {\omega _{1j}}\left( {{t_{{\text{cpa}}2}}} \right) $ 为阵元1接收信号第 $ j $ 条干涉条纹在阵元2最近通过时刻 $ {t_{{\text{cpa2}}}} $ 处的频率。式（12）避免了条纹斜率求导，仅利用 $ {t_{{\text{cpa}}1}} $ 和 $ {t_{{\text{cpa2}}}} $ 两时刻的频率和距离关系，它可应用于临近CPA位置处的时间-频率域具有一定曲率干涉条纹。

将式（10）和式（12）联立有如下表达：

$ {r_{{\text{cpa}}1}} = \sqrt {\frac{{{{\left( {d\cos \varphi } \right)}^2}}}{{\left[ {{{\left( {{{{\omega _{1j}}\left( {{t_{{\text{cpa}}1}}} \right)} \mathord{\left/ {\vphantom {{{\omega _{1j}}\left( {{t_{{\text{cpa}}1}}} \right)} {{\omega _{1j}}\left( {{t_{{\text{cpa}}2}}} \right)}}} \right. } {{\omega _{1j}}\left( {{t_{{\text{cpa}}2}}} \right)}}} \right)}^{ - {2 \mathord{\left/ {\vphantom {2 \beta }} \right. } \beta }}} - 1} \right]}}}。$

(13)

由此可获得目标相对阵元1的最近通过距离。同理，目标在时刻 $ {t_{{\text{cpa}}1}} $ 相对阵元2的距离为：

$ {r_2}\left( {{t_{{\text{cpa}}1}}} \right) = \sqrt {r_{{\text{cpa}}2}^2 + {{\left( {d\cos \varphi } \right)}^2}} 。$

(14)

阵元2接收信号干涉条纹上的频率和距离之间关系如下式：

$ \omega _{2i}^{ - {1 \mathord{\left/ {\vphantom {1 \beta }} \right. } \beta }}\left( {{t_{{\text{cpa}}2}}} \right){r_{{\text{cpa}}2}} = \omega _{2i}^{ - {1 \mathord{\left/ {\vphantom {1 \beta }} \right. } \beta }}\left( {{t_{{\text{cpa}}1}}} \right){r_2}\left( {{t_{{\text{cpa}}1}}} \right) 。$

(15)

式中： $ {\omega _{2i}}\left( {{t_{{\text{cpa}}2}}} \right) $ 为阵元2接收信号第 $ i $ 条干涉条纹在最近通过时刻 $ {t_{{\text{cpa}}2}} $ 处的频率； $ {\omega _{2i}}\left( {{t_{{\text{cpa}}1}}} \right) $ 为阵元2接收信号第 $ i $ 条干涉条纹在目标相对阵元1最近通过时刻 $ {t_{{\text{cpa}}1}} $ 处的频率。联立式（14）和式（15）有：

$ {r_{{\text{cpa}}2}} = \sqrt {\frac{{{{\left( {d\cos \varphi } \right)}^2}}}{{\left[ {{{\left( {{{{\omega _{2i}}\left( {{t_{{\text{cpa}}2}}} \right)} \mathord{\left/ {\vphantom {{{\omega _{2i}}\left( {{t_{{\text{cpa}}2}}} \right)} {{\omega _{2i}}\left( {{t_{{\text{cpa}}1}}} \right)}}} \right. } {{\omega _{2i}}\left( {{t_{{\text{cpa}}1}}} \right)}}} \right)}^{ - {2 \mathord{\left/ {\vphantom {2 \beta }} \right. } \beta }}} - 1} \right]}}} 。$

(16)

那么水平阵几何中心处最近通过时间和最近通过距离可表示为：

$ {t_{{\text{cpa}}}} = {{\left( {{t_{{\text{cpa}}1}} + {t_{{\text{cpa}}2}}} \right)} \mathord{\left/ {\vphantom {{\left( {{t_{{\text{cpa}}1}} + {t_{{\text{cpa}}2}}} \right)} 2}} \right. } 2}，$

(17)

$ {r_{{\text{cpa}}}} = {{\left( {{r_{{\text{cpa}}1}} + {r_{{\text{cpa}}2}}} \right)} \mathord{\left/ {\vphantom {{\left( {{r_{{\text{cpa}}1}} + {r_{{\text{cpa}}2}}} \right)} 2}} \right. } 2} 。$

(18)

根据BTR图做Hough变换得到的阵中心的距速比，可获得目标声源运动速度 $ {v_0} $ ，于是目标与阵中心的水平距离随时间的变化关系如下：

$ r\left( t \right) = \sqrt {v_0^2{{\left( {t - {t_{{\text{cpa}}}}} \right)}^2} + r_{{\text{cpa}}}^2}。$

(19)

通过式（19）可计算出各时刻目标与阵中心之间的距离。

2 数值仿真

浅海水深为100 m，水体声速等梯度为1500 m/s，海底为半空间，声速为1650 m/s，密度为1850 kg/m³，衰减系数为0.8 dB/λ。声源深度 $ {z_s} $ 为50 m，水听器深度 $ {z_r} $ 为70 m。运动目标相对水平阵的航向角 $ \varphi $ 为45°，行进速度 $ {v_0} $ 为5 m/s。阵元1最近通过距离 $ {r_{{\text{cpa}}1}} $ 为1000 m，最近通过时间 $ {t_{{\text{cpa}}1}} $ 为200 s，阵元2最近通过距离 $ {r_{{\text{cpa}}2}} $ 为1500 m，最近通过时间 $ {t_{{\text{cpa}}2}} $ 为300 s，阵元1与阵元2的间距 $ d $ 为500 $ \sqrt 2 $ m。目标与接收阵元的相对位置示意图见图1，运动时间区间0～400 s，宽带声源频率区间100～350 Hz。阵元1和阵元2的干涉条纹如图4所示。

可知，阵元1的最近通过时间估计值 $ {\hat t_{{\text{cpa}}1}} $ 为200 s，阵元2的最近通过时间 $ {\hat t_{{\text{cpa2}}}} $ 为300 s，则阵中心的最近通过时间 $ {\hat t_{{\text{cpa}}}} $ 为250 s。选取阵元1时频域中的第 $ j $ 条条纹，在200 s和300 s对应的能量峰值处频率分别为 $ {\hat \omega _{1j}}\left( {{{\hat t}_{{\text{cpa}}1}}} \right) = 293 $ Hz和 $ {\hat \omega _{1j}}\left( {{{\hat t}_{{\text{cpa2}}}}} \right) = 328 $ Hz。选取阵元2时频域中的第 $ i $ 条条纹，在200 s和300 s对应的能量峰值处频率分别为 $ {\omega _{2i}}\left( {{{\hat t}_{{\text{cpa}}1}}} \right) = 314 $ Hz和 ${\omega _{2i}}\left( {{{\hat t}_{{\text{cpa2}}}}} \right) = 299$ Hz。根据式（13）和式（16）可计算获得 $ {\hat r_{{\text{cpa}}1}} $ 为993.7057 m和 $ {\hat r_{{\text{cpa2}}}} $ 为1559.1 m，相对误差分别为0.63%和3.9%。阵元1和阵元2的最近通过距离估计误差产生的主要原因是峰值频率和波导不变量假设值偏差造成的。阵中心的最近通过距离 $ {\hat r_{{\text{cpa}}}} $ 为1276.4 m，由距速比 ${{r_{{\text{cpa}}}}} / {{v_0}}$ 为250可获得目标运动速度 $ {\hat v_0} = 5.1 $ m/s。根据式（19）可得运动目标与接收阵中心的估计距离随时间的变化关系，如图5所示。

可知，运动目标与接收阵中心的估计距离与真实值存在一定偏差，各时间点的距离估计相对误差为2%，不随时间的变化而变化。由于航向角和距速比是真实值，由式（19）可知，距离估计相对误差与目标运动速度估计相对误差大小一致。

2.1 距速比对估计结果的影响

考虑使用Hough变换获得距速比的估计误差对距离估计结果的影响。运动目标的航向角为真实值45°，讨论距速比估计值分别为230、250、270时，运动目标与接收阵中心的距离估计结果如图6所示。

可知，距速比会影响距离估计结果的偏离程度。随着距速比的增加，曲线开口逐渐增大。当距速比大于真实值250时，曲线位于实际距速比曲线的下方，反之，小于真实值的曲线位于实际距速比曲线的上方。各距速比下的距离估计值与距离真实值的相对误差随着时间的变化曲线，如图7所示。

可知：1）当时间对应最近通过时间250 s时，此时估计结果只受最近通过距离估计值的影响，无论距速比估计值为何值，相对误差都为2%左右。2）当距速比为真实值250时，距离估计相对误差结果不随时间的变化而发生改变，保持在2%左右。3）当距速比小于真实值时，距离估计相对误差随着时间远离 $ {t_{{\text{cpa}}}} $ 逐渐增大，0时刻误差达到最大为7%左右。由式（19）可知，最近通过距离和速度一定，当时间偏离 $ {t_{{\text{cpa}}}} $ ，误差会随之增加。4）当距速比大于真实值时，随着时间远离 $ {t_{{\text{cpa}}}} $ ，相对误差先减小后增大。由于最近通过距离估计值偏小，当时间偏离 $ {t_{{\text{cpa}}}} $ ，距离估计值趋近真实值，之后逐渐出现偏差。以上现象分析是针对目标速度估计值大于真实值的情况，若其小于真实值，则情况恰好相反。

2.2 航向角对估计结果的影响

考虑不同航向角估计值对距离估计结果的影响。假设距速比为真实值，航向角估计值分别为40°、45°、50°情况下，距离估计结果如图8所示。

可知，随着航向角估计值相较真实值45°偏大（偏小），式（13）和式（15）的最近通过距离估计结果偏小（偏大），相对接收阵中心的目标距离估计结果较实际值偏小（偏大），因此估计距离随时间变化曲线下侧（上侧）。不同航向角的距离估计相对误差，如图9所示。可知，距离估计相对误差会保持在一个定值，与时间变化无关。这是由于距速比是真实值，距离估计相对误差仅与目标运动速度估计误差有关。

由上述分析可知，距速比和航向角的估计偏差会在一定程度上影响目标距离估计结果，在实际中应尽量保证距速比和航向角先验信息的准确性。

3 结　语

本文提出一种基于简正波水平波数差与波导不变量关系式的被动测距方法。利用运动目标声源通过CPA位置在水平阵两端水听器产生的宽带连续谱干涉条纹，获得其最近通过时间 $ {t_{{\text{cpa}}1}} $ 和 $ {t_{{\text{cpa2}}}} $ ，再利用几何关系及简正波水平波数差和波导不变量关系式，估计出阵元1和阵元2的最近通过距离 $ {r_{{\text{cpa}}1}} $ 和 $ {r_{{\text{cpa2}}}} $ 以及目标运动速度 $ {v_0} $ ，最终估计出运动目标相对接收阵几何中心的距离。数值仿真结果表明，方法可有效利用临近CPA位置处的时间-频率域具有一定曲率干涉条纹，估计出运动目标相对接收阵几何中心的最近通过距离及目标相对水听器距离随时间的变化关系。目标距速比估计值和目标航向角估计值会一定程度影响目标距离估计结果，在实际中应尽量保证距速比和航向角先验信息的准确获取。