1.   引言

随着海洋探索需求增强，水下传感网络技术正在得到越来越多的关注，而为水下传感节点提供定位与导航服务正是水下传感网络的主要任务之一^[1]。对于水下传感网络定位技术，近年来国内外学者的关注点主要集中在各节点间是否需要时钟同步^[2]、是否考虑声速分布不均匀因素^[3-6]、定位协议^[7-9]以及定位算法的附加功能^[3,4,10,11]等方面。文献[2]提出了一种各节点时钟同步情况下基于双程时间测量的水声定位算法，该算法在常声速梯度假设下，采用高斯-牛顿法对目标位置进行估计。文献[3]提出了联合解决水下传感网络中定位和时钟同步的算法，待定位节点与信标节点间采用应答方式获取双程传播时延，并引入一个积分常数来修正声速分层的影响。文献[4]提出了基于双程时延测量的位置和时钟联合估计方法，采用等效常梯度声速剖面来解决声速分布不均带来的定位误差，但定位过程较复杂，周期较长。文献[5]提出了各节点间时钟已同步条件下考虑声线弯曲的水下传感网络定位算法，通过卡尔曼滤波优化接收到的时延信息，再采用凸优化方法对待定位节点位置进行估计。文献[6]提出了基于双程时延测量并结合声线跟踪的定位算法，采用高斯-牛顿法估计最优的待定位节点位置，但在信标节点布阵不佳时，算法发散。文献[7]提出了基于移动预测的水下传感网络定位算法(Scalable Localization scheme with Mobility Prediction, SLMP)，待定位节点根据自身的运动模型对自身位置进行预测，并通过信标节点的信息进行位置和运动模型的更新。文献[8,9]提出了通过单信标测距对水下潜器进行定位的方法，潜器通过自身配备的惯性导航设备获取自身轨迹，再结合行进时获取的与信标间到达时间(Time Of Arrival, TOA)或到达角(Angle Of Arrival, AOA)信息，采用虚拟长基线方法估计出自身绝对位置。文献[10]提出了通过最小二乘联合估计所有节点的时钟信息和待定位节点与各信标节点间距离，进而获取待定位节点位置的算法。文献[11]提出了一种分布式逐级网络节点定位的方法，通过对已定位节点的定位误差进行评估以将其升级为新的信标节点，从而扩大定位范围。

上述方法考虑了水下传感网络定位中的大部分情形，但均需要待定位节点反复与信标节点通信。在需要对多个待定位节点同时进行定位服务，或待定位节点因任务需求需要静默的情况下，上述方法均不适用。针对这些情况，文献[12]首次提出静默定位(Underwater Positioning Scheme, UPS)算法，该算法不需节点间时钟同步且待定位节点可全程处于静默接收状态，但该算法未考虑声速不均匀分布问题，且算法通过降阶法对待定位节点位置进行求解，存在定位盲区的问题。文献[13]采用增加信标节点的方式解决UPS算法的定位盲区的问题，但增加信标节点无疑会大大提高算法的应用成本。文献[14]对UPS算法进行了复现，并通过实验加以验证。虽然后续的研究人员在静默定位方面做出了很大的努力，但并没有解决最初的问题，即声速分布不均匀对UPS算法的影响，且没有从本质上解决定位盲区的问题。尽管在无线传感网络定位领域，文献[15]已经解决了定位盲区的问题，然而对于水下传感网络，此方法所需的通信量和计算量依然很大，这也导致其无法直接在水下传感网络中使用。

本文针对现有UPS算法的不足，提出了一种新的静默式水下传感网络定位算法。通过声线跟踪技术，解决了声速不均匀分布的问题。采用高斯-牛顿法估计最优的待定位节点位置，解决了UPS算法存在定位盲区的问题。同时针对水下传感网络各参考节点布放时容易出现的定位阵型不佳问题，引入了改进的Tikhonov正则化方法，根据迭代效果反馈控制正则化参数，提高了算法的稳健性。最后通过仿真分析，讨论了各个步骤的必要性和各项误差对算法性能的影响，验证了算法的有效性。

2.   算法原理

假定待定位节点T的通信范围内布设有4个信标节点，分别标号为A, B, C和D，其3维坐标分别表示为${{X}}={\left[ {{x_{\rm{t}}},{y_{\rm{t}}},{z_{\rm{t}}}} \right]^{\rm{T}}}$和${{{P}}_i} = {\left[ {{x_i},{y_i},{z_i}} \right]^{\rm{T}}}, i = {\rm{a,b,c,d}}$。其位置关系如图1所示。

图 1  水下传感网络静默定位流程

下载: 全尺寸图片 幻灯片

假设信标节点A为定位发起节点，该网络内信息传递的顺序为A-B-C-D。在检测到待定位节点T处于该区域内时，信标节点A周期性广播定位数据包，其中包括其当前的时钟、经纬度和深度等信息，待定位节点T于自身时钟时刻${t_{\rm{a}}}$接收到来自信标节点A的定位发起信号。后面的信标节点依次在接收到前一信标节点发出的信号的一段时间${\tau _i},i = {\rm{b,c,d}}$后，广播其自身的定位数据包和处理延迟${\tau _i}$，待定位节点T于自身时钟时刻${t_i},i = {\rm{b,c,d}}$接收到来自后面各信标节点的定位信号。若通信范围内存在更多的信标节点，该流程则依次进行下去

同理，各相邻信标节点之间的距离也可以按照式(1)的格式表示。

假设相邻信标节点之间水平距离为5 km，1个定位周期内，各节点自身时钟歪斜所引起的时间估计误差可以忽略不计，故本文定位算法不考虑时钟歪斜问题。

由于不同深度处的温度、盐度和压力等参数是不同的，水下各深度处声速分布不均匀。这导致了声波在水下不会以直线形式传播，也就是说传播时间可以被表示成声源与接收端位置的函数，即

式中，${t_{{\rm{T}}i}},i = {\rm{a,b,c,d}}$表示声信号从信标节点i直接传播到待定位节点T所需的时间。此时，求解待定位节点位置${{X}}$的问题可以转化为求解最符合观测数据的一组$\left\{ {{t_{{\rm{T}}i}},i = {\rm{a,b,c,d}}} \right\}$的问题。

以测量的时间量为核心建立观测方程，有

式中，${{F}}\left( {{X}} \right)$和${{R}}$分别称为估计向量和观测向量，可以写成

式中，${\left[ \cdot \right]^{\rm{T}}}$表示矩阵的转置；$\Delta {t_j},j = 2,3,4$表示声信号从信标节点A发出，最后经由信标节点i发送给待定位节点T与信标节点A直接发送给待定位节点T的声信号传播时间差，可以写成

式(5)中，${{\rm{tt}}_i},i = {\rm{b,c,d}}$，表示声信号从信标节点A发出，最后经由信标节点i发送给待定位节点T期间，声信号在信标节点间传播的时间，可以写成

此时，求解待定位节点位置${{X}}$的问题就变成了如下的最小化问题

可以看到，式(8)是典型的最小二乘问题，可以采用高斯-牛顿迭代方法得出最优解。求解过程中，第k次迭代中状态向量${{X}}$的变化过程为

式中，${\alpha _k}$表示第k次迭代的步长，可以通过Armijo准则选取^[16]；$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)$表示${{X}} = {{{X}}_{k - 1}}$时目标函数${{G}}\left( {{X}} \right)$的雅克比矩阵。由射线声学理论，第k次迭代的雅克比矩阵$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)$可以写成如式(10)的形式^[17]

式中，${n_{{\rm{T}}i}}$为这条声线的snell常数，即$n = \dfrac{{\cos {\vartheta _{{z_i}}}}}{{c\left( {{z_i}} \right)}} =$$\dfrac{{\cos {\vartheta _z}}}{{c\left( z \right)}}$, ${\vartheta _z}$为这条声线在深度z处的掠射角，即声线切线方向与水平方向夹角。

按照式(9)的形式进行迭代求解，目标函数${\left\| {{{G}}\left( {{X}} \right)} \right\|^2} = 0$对应的${{{X}}_k}$即为待定位节点的真实位置。但为提高计算效率和算法稳健性，可设定目标函数小于门限$\varepsilon$时，即终止迭代，当前的${{{X}}_k}$即作为待定位节点坐标。其中，$\varepsilon$表征算法最终得出的待定位节点坐标精度，可根据应用需求由${\left\| {{{G}}\left( {{X}} \right)} \right\|^2} =$${\left\| {{{F}}\left( {{X}} \right) - {{R}}} \right\|^2}$计算得出。

值得注意的是，若信标节点组成的定位阵型不佳，比如信标节点都布置在同一深度附近时，式(9)中的$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}}$矩阵可能奇异或病态，实际表现为其条件数过大，使得算法发散^[18]。可以采用Tikhonov正则化来解决这一问题，式(9)需要改写成

式中，$\eta$为正则化参数，通过设置合适的$\eta$可以使得$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}} + \eta {{I}}$的条件数降低，使得算法收敛；但Tikhonov正则化为有偏估计，这就导致了最终解算的位置与真实位置存在偏差。正则化参数$\eta$越大，矩阵$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}} +$$\eta {{I}}$的条件数越小，但Tikhonov正则化带来的偏差会越大。因此，每一次迭代的正则化参数$\eta$需要根据情况随时变化。

该算法流程如下：

(1) 在不考虑声线弯曲的情况下，得到待定位节点位置粗估计值，将其设为本算法的迭代初值${{{X}}_0}$。计算过程中，若$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}}$的条件数较大，则对其进行Tikhonov正则化处理后，再继续计算；

(2) 根据当前估计的待定位节点位置进行声线跟踪，得出声信号从各信标节点${{{P}}_i}$传播到当前估计位置所需的时间${t_{{\rm{T}}i}}$；

(3) 通过式(4)和式(5)计算当前的目标函数${\left\| {{{G}}\left( {{X}} \right)} \right\|^2}$；

(4) 将得到的目标函数与事先设定的门限相比较。若${\left\| {{{G}}\left( {{X}} \right)} \right\|^2} \le \varepsilon$，则认为当前估计位置记为待定位节点真实位置；若${\left\| {{{G}}\left( {{X}} \right)} \right\|^2} > \varepsilon$，则通过式(10)计算当前迭代周期的雅克比矩阵$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)$；

(5) 计算$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}}$的条件数。若条件数小于门限$\xi$，则认为雅可比矩阵$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)$正常，将其代入式(9)更新待定位节点位置，并重复第(2)～(4)步；若条件数大于门限$\xi$，将其代入式(11)更新待定位节点位置，并执行第(2), (3)步。若目标函数${\left\| {{{G}}\left( {{X}} \right)} \right\|^2}$降低，则令$\eta = \eta /{\kappa _1}$，并执行第(4)步；若目标函数${\left\| {{{G}}\left( {{X}} \right)} \right\|^2}$未降低，则令$\eta = \eta \times {\kappa _2}$同时减小步长${\alpha _k}$，重新代入式(11)更新待定位节点位置，并执行第(2), (3)步，直到更新的状态${{{X}}_k}$使得目标函数${\left\| {{{G}}\left( {{X}} \right)} \right\|^2}$降低为止。其中，$\xi$为矩阵病态程度判决门限，对其取值方法的讨论，请见本文仿真分析部分；${\kappa _1}$和${\kappa _2}$为大于1的常数，根据经验选取即可。

该算法流程如图2所示。

图 2  基于射线声学的水下传感网络静默定位算法流程图

下载: 全尺寸图片 幻灯片

3.   仿真分析

下面对所提算法进行仿真，并与相关算法进行对比。

3.1   定位区域及解的唯一性

设置4个信标节点的位置分别为(–2500, 2500, 3000) m, (–2500, –2500, 2010) m, (2500, –2500, 4000) m和(2500, 2500, 1050) m。为不失一般性，设置待定位节点实际位置坐标为(1000, 500, 100) m。设置按照Munk典型声速剖面设置水体声速，水深5000 m。在100 m深度平面上的目标函数${\left\| {{{G}}\left( {{X}} \right)} \right\|^2}$分布如图3(a)所示；在$y = x - 500$的坐标平面上的目标函数分布如图3(b)所示。

图 3  布阵区域内目标函数分布情况

下载: 全尺寸图片 幻灯片

由图3可见，理想布阵情况下目标函数在全空间内存在唯一最小值点，即为待定位节点真实位置，不存在多解情况。相比于UPS算法降阶方程组所带来的局限性导致其在信标节点附近存在定位盲区，且需要人工排除多解的情况，本文算法既可获得唯一解，又可在4信标通信范围内任意位置对待定位节点进行定位。

3.2   信标定位阵型不佳

当各信标节点组成的定位阵型不佳时，定位算法的性能会变得极差甚至失效。以各信标节点布置在同一深度附近为例，此时定位阵在深度维的孔径会变得很小，雅可比矩阵变得不满秩，进而导致$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}}$矩阵病态甚至奇异，从而使得算法发散。

下面针对定位阵型不佳这一情形，对本文算法进行仿真。设置4个信标节点的位置分别为(–2500, 2500, 3050) m, (–2500, –2500, 3060) m, (2500, –2500, 3070) m和(2500, 2500, 3080) m，正则化参数$\eta$的初始值设置为${10^{ - 7}}$, ${\kappa _1}$和${\kappa _2}$分别设置为10和5。其余仿真条件与前文一致，本文算法的迭代情况如图4(a)所示。值得注意的是，在第1次高斯-牛顿迭代中，$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}}$的条件数为$5.84 \times {10^4}$，矩阵明显是奇异的；经过本文改进的Tikhonov正则化处理，$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}} +$$\eta {{I}}$的条件数降至15.6。可以看出，对于定位阵型严重不佳的情况，本文算法仍能较好地对待定位节点进行定位。

图 4  定位阵型不佳时本文算法迭代情况

下载: 全尺寸图片 幻灯片

下面以本例第1次迭代为例，讨论步骤(5)中矩阵病态程度判决门限$\xi$的取值方法。在设置不同$\xi$值时，一次迭代后目标函数情况如图4(b)所示。门限$\xi$值较小时，每次迭代会采用较大的正则化参数$\eta$，虽然可降低矩阵$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right){\left( {\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}} + \eta {{I}}$的病态程度，但由于Tikhonov正则化的有偏性，迭代后的目标函数仍会较大，甚至大于迭代前的目标函数；而门限$\xi$值过大时，又会导致$\nabla {{G}}\left( {{{{X}}_{k - 1}}} \right)\cdot$${\left( {\nabla {{G}} \left( {{{{X}}_{k - 1}}} \right)} \right)^{\rm{T}}}+ \eta {{I}}$矩阵仍然病态甚至奇异，Tikhonov正则化失效。因此，矩阵病态程度判决门限$\xi$应选取算法性能稳健的部分，以本例为例，矩阵病态程度判决门限$\xi$的取值范围应为$\left[ {{\rm{10,1}}{{\rm{0}}^{\rm{4}}}} \right]$。

3.3   误差分析

由本文第2节算法原理可以看出，本文算法定位精度的主要影响因素包括信标节点位置误差、测时误差和水下环境失配等。下面采用上述仿真条件，对算法涉及的各输入参数引入一定的高斯随机误差，并进行200次蒙特卡洛仿真，以分析各参数的误差对本文算法性能的影响，并与其他算法进行对比。

3.3.1   信标节点位置误差

对水下待定位节点定位前，需对各信标节点的位置进行标校，以得到信标节点的精确位置，包含水平位置和深度信息两方面。因此，信标节点位置误差对本文算法的影响，需要分信标节点的水平位置误差和深度误差两个方面进行讨论。

信标节点深度一般由配备的深度传感器测得，因此其深度误差来源于配备的深度传感器误差。本算法和UPS算法受信标节点深度误差影响如图5(a)所示。

图 5  各误差来源对定位结果影响

下载: 全尺寸图片 幻灯片

水下信标节点的水平位置需要由测量船事先标校获得，一般地，测量船绕信标节点布放位置以圆形等对称航迹航行，期间多次对信标节点进行应答式测距，再结合测量船自身GPS信息，采用反转长基线原理等完成对信标节点水平位置的标校^[17]。信标节点的水平位置误差来源于多方面，包括信标节点深度误差、测量船测时误差、声速误差等。本文算法和UPS算法受信标节点的水平位置误差影响如图5(b)所示。

可以看出，随着信标水平位置误差和深度误差的增大，本文算法的定位误差也随之增大。由于UPS算法定位误差远大于本文算法，相比之下，其对于信标的水平位置误差和深度误差均不十分敏感。

可见，信标节点位置的准确性将直接影响本文算法的定位精度。因此，在应用本文算法前，获取信标节点的锚定准确位置，并考虑信标节点在水下位置的随机摆动，对于提高本文算法的定位精度是十分必要的。

3.3.2   测时误差

本文算法属于基于时间测量的定位算法，信号到达时间的测量误差会直接影响算法的定位误差。

本文算法中，采用双曲调频(Hyperbolic Frequency Modulated, HFM)信号通过匹配滤波实现对信号到达时间的测量。匹配滤波器的时间估计精度为$1/B$，其中B为信号的带宽，并且，各节点间的多普勒和水声信道的多途效应也会对到达时间测量精度造成一定的影响。信号到达时间测量误差对本文算法和UPS算法的影响如图5(c)所示。可以看出，随着测时误差增大，本文算法定位误差也随之增大，但误差仍显著小于UPS算法。

3.3.3   水下环境参数失配

本文算法中声线跟踪部分依赖于相对准确的水下环境参数输入，考虑到水下传感网络的应用场景主要在深海中，水底地形较为平坦，因此声速剖面的误差成为影响本文算法精确度的主要因素。

声速剖面的误差分为与深度无关的测量误差和与深度相关的时变偏差。此处引入剖面等效误差概念来表征声速剖面的失配程度^[19]，令

其中，$Z$为水深，${c_{{\rm{used}}}}\left( z \right)$为算法中使用的声速剖面，${c_{{\rm{true}}}}\left( z \right)$为水下真实的声速剖面。

对算法中使用的1000 m以上水深的声速剖面加入一定的测量误差和时变误差，其最终的定位误差随剖面等效误差的变化如图6所示。可见，本文算法对声速剖面误差较为敏感，后续研究将针对此问题进一步展开。

图 6  定位误差随剖面等效误差变化关系

下载: 全尺寸图片 幻灯片

4.   结论

针对已有静默定位算法存在定位盲区、多解以及未考虑声速不均匀带来的影响等问题，本文提出了一种改进的静默式水下传感网络定位方法。算法通过融合声线跟踪技术和高斯-牛顿迭代方法使得算法可以在信标阵列通信范围内任意位置获得唯一的待定位节点最优位置，且修正了声速分层分布对定位结果带来的影响。同时，采用引入反馈控制机制的Tikhonov正则化方法，使得算法在信标阵型不佳的情况下，仍然保持较好的定位性能。仿真结果表明：该算法较已有的UPS算法具有更好的性能，且输入参数存在一定误差的情况下，仍能保持较好的定位效果，同时算法对信标布放阵型具有一定的稳健性。