0.   引　言

随着海洋探测领域的快速发展，水下无线传感器网络(underwater wireless sensor networks, UWSN)越来越受到学术和工业研究人员的关注，已广泛应用于海洋监测、资源开发、灾害预警、海底测绘、环境保护等领域^[1]. 其中，准确的水下传感器位置信息对于完成上述任务至关重要. 在UWSN中，一般使用光波和声波信号来进行节点定位^[2]. 虽然光信号具有高宽带、低延迟、高可靠的特性，但其在水下的传输距离却极为有限. 由于声信号相对于光信号在水中传播时能量损失较小，从而可以实现远距离传输^[3]. 在水声定位技术中，长基线定位系统具有定位精度高和应用场景广的优势.

在水下定位方法中，根据接收的信号观测量，可以分为基于接收信号强度(received signal strength, RSS)、时延(time delay, TD)、到达时间差(time difference of arrival, TDOA)、到达角度(angle of arrival, AOA)、以及多普勒频移(Doppler shift, DS)的定位方法^[4]. 近年来，基于时延和多普勒频移测量的水下定位方法得到了广泛应用. 在文献[5-6]中，针对测量模型的非线性特性，提出了一种基于迭代搜索的最大似然估计(maximum likelihood estimation, MLE)方法. 然而，此方法中初始解是随机选择的，估计结果容易收敛到局部最优甚至发散. 文献[7-8]中提出了基于定位问题的闭合解(closed form solutions, CFS)方法，不需要初始解且计算复杂度较低. 然而，该方法不能很好地工作于高噪声的环境. 为了提高算法在高噪声环境下的适应性，文献[9-10]中提出了半定规划(semi-definite programming, SDP)算法，该方法可以将目标定位的非凸问题重构为凸问题进行求解. 然而，该算法的解是次优的且计算复杂度较高.

在实际的水下环境中，声传感器在信号观测期间内的运动效应相对于声波信号速度比较显著. 因此，在建立测量模型时，该运动效应不可忽略. 文献[11]考虑了传感器在信号观测期间内的运动效应，提出了一种新的TD和DS的观测模型，并分别给出了该模型下的CFS和SDP算法. 基于此工作，文献[12]在递归的测量模型基础上，提出了一种约束加权最小二乘法的定位方法. 此外，文献[13]考虑了二阶噪声对于位置估计的影响. 文献[14]理论分析了在信号观测期间内忽略传感器运动效应引起的定位模型偏差.

由于水下环境的特殊性，上述文献提出的时延和多普勒频移测量模型不能直接用于水下声传感器网络中. 一是水下声速与深度、温度、密度等许多因素相关，出现声线弯曲现象，即分层效应^[15]. 二是由于水流和潮汐等许多因素影响会导致传感器出现缓慢运动，造成位置和速度偏差^[16]. 因此，有必要研究一种能够结合声线和运动补偿的水下目标定位方法，并考虑传感器位置和速度偏差，使其更加适应于实际水下环境，以提高目标定位精度.

鉴于此，本文基于水下声传感器网络，采用长基线来实现目标定位.

1.   测量模型

水下声传感器网络主要包括M个传感器和一个待定位的目标. 其中，传感器作为移动的锚节点用于目标定位，目标位置为${{\boldsymbol{u}}^{\mathrm{o}}} = {[{x^{\mathrm{o}}},{y^{\mathrm{o}}},{z^{\mathrm{o}}}]^{\text{T}}}$.传感器i的初始位置为${\boldsymbol{t}}_i^{\mathrm{o}} = {[x_{\text{A}}^i,y_{\text{A}}^i,z_{\text{A}}^i]^{\text{T}}}$，$ i = 1,2, \cdots ,{{M}} $. 假设其在信号观测期间内以速度${\boldsymbol{v}}_i^{\mathrm{o}} = {[{v_{i,x}},{v_{i,y}},{v_{i,z}}]^{\text{T}}}$进行线性移动. 由于传感器的运动效应，传感器发射信号时的位置(${\boldsymbol{t}}_i^{\mathrm{o}}$)和信号反射回来被接收时的位置(${\boldsymbol{s}}_i^{\mathrm{o}}$)是不同的，其运动示意图如图1所示.

图  1  测量模型示意图

下载: 全尺寸图片 幻灯片

传感器i发出的探测信号到达目标后并反射回来的理想时延$ \tau _i^{\mathrm{o}} $可以表示为

式中：$\tau ({{\boldsymbol{u}}^{\mathrm{o}}},{\boldsymbol{t}}_i^{\mathrm{o}})$为探测信号从传感器i到目标的时延；$ \tau ({{\boldsymbol{u}}^{\mathrm{o}}},{\boldsymbol{s}}_i^{\mathrm{o}}) $为从信号反射回来到传感器i的时延.

根据传感器在信号观测期间内线性运动的假设，可以得到${\boldsymbol{s}}_i^{\mathrm{o}} = {\boldsymbol{t}}_i^{\mathrm{o}} + {\boldsymbol{v}}_i^{\mathrm{o}}\tau _i^{\mathrm{o}}$，表示为${\boldsymbol{s}}_i^{\mathrm{o}} = {[x_{\text{B}}^i,y_{\text{B}}^i,z_{\text{B}}^i]^{\text{T}}}$. 由于${\boldsymbol{s}}_i^{\mathrm{o}}$与$ \tau _i^{\mathrm{o}} $有关，因此$ \tau _i^{\mathrm{o}} $的表达式是递归的. 将测量得到的时延信息组合成向量形式为

式中：$ {\tau _i} $表示传感器i测量得到的时延信息；${{\boldsymbol{\tau }}^{\mathrm{o}}}$表示真实的时延信息，${{\boldsymbol{\tau }}^{\mathrm{o}}} = {[\tau _1^{\mathrm{o}}, \cdots ,\tau _M^{\mathrm{o}}]^{\text{T}}}$；$\Delta {\boldsymbol{\tau }}$为时延信息的测量噪声，$ \Delta {{\boldsymbol{\tau }}^{\mathrm{o}}} = {[\Delta {\tau _1}, \cdots ,\Delta {\tau _M}]^{\text{T}}} $，$ \Delta {\tau _i} $服从均值为0，方差为$ {\sigma _\tau } $的高斯分布.

通过对式(1)两边时间求导数，并乘以探测信号载频得到多普勒频移的真实值$ f_i^{\mathrm{o}} $，如下所示：

式中：$ f_i^{\text{c}} $为传感器i发出的探测信号载频.

同理，将测量的多普勒频移信息组合成向量形式为

式中：$ {f_i} $表示传感器i测量得到的多普勒频率信息；${{\boldsymbol{f}}^{\mathrm{o}}}$为真实的多普勒频率信息，${{\boldsymbol{f}}^{\mathrm{o}}} = {[f_1^{\mathrm{o}}, \cdots ,f_M^{\mathrm{o}}]^{\text{T}}}$；$ \Delta {\boldsymbol{f}} $为多普勒频率的测量噪声，$ \Delta {{\boldsymbol{f}}^{\mathrm{o}}} = {[\Delta {f_1}, \cdots ,\Delta {f_M}]^{\text{T}}} $，$ \Delta {f_i} $服从均值为0，方差为$ {\sigma _f} $的高斯分布.

进一步地，将时延和多普勒频率信息组合成向量的形式，可以表示为

式中：${{\boldsymbol{\alpha }}^{\mathrm{o}}} = {[{({{\boldsymbol{\tau }}^{\mathrm{o}}})^{\text{T}}},{({{\boldsymbol{f}}^{\mathrm{o}}})^{\text{T}}}]^{\text{T}}}$，$\Delta {\boldsymbol{\alpha }} = {[\Delta {{\boldsymbol{\tau }}^{\text{T}}},\Delta {{\boldsymbol{f}}^{\text{T}}}]^{\text{T}}}$服从零均值的高斯分布，其协方差矩阵为${{\boldsymbol{Q}}_{\boldsymbol{\alpha }}} = \boldsymbol{E}[\Delta {\boldsymbol{\alpha }}\Delta {{\boldsymbol{\alpha }}^{\text{T}}}]$.

考虑到传感器节点位置和速度的偏差，所有传感器节点的测量位置和速度如下所示：

式中：${\boldsymbol{t}}$和${\boldsymbol{v}}$分别表示测量得到的传感器位置和速度，具体地，${\boldsymbol{t}} = {[{{\boldsymbol{t}}_1}, \cdots ,{{\boldsymbol{t}}_M}]^{\text{T}}}$，${\boldsymbol{v}} = {[{{\boldsymbol{v}}_1}, \ldots ,{{\boldsymbol{v}}_M}]^{\text{T}}}$；${{\boldsymbol{\;\beta }}^{\mathrm{o}}} = [{({{\boldsymbol{t}}^{\mathrm{o}}})^{\text{T}}}, {({{\boldsymbol{v}}^{\mathrm{o}}})^{\text{T}}}]$表示真实传感器位置和速度，具体地，${{\boldsymbol{t}}^{\mathrm{o}}} = {[{\boldsymbol{t}}_1^{\mathrm{o}}, \cdots ,{\boldsymbol{t}}_M^{\mathrm{o}}]^{\text{T}}}$，${{\boldsymbol{v}}^{\mathrm{o}}} = {[{\boldsymbol{v}}_1^{\mathrm{o}}, \ldots ,{\boldsymbol{v}}_M^{\mathrm{o}}]^{\text{T}}}$；$\Delta {\boldsymbol{\beta }} = {[\Delta {{\boldsymbol{t}}^{\text{T}}},\Delta {{\boldsymbol{v}}^{\text{T}}}]^{\text{T}}}$表示传感器位置和速度偏差；协方差矩阵为${{\boldsymbol{Q}}_{\boldsymbol{\beta }}} = \boldsymbol{E}[\Delta {\boldsymbol{\beta }}\Delta {{\boldsymbol{\beta }}^{\text{T}}}]$，$\Delta {\boldsymbol{t}} = {[\Delta {{\boldsymbol{t}}_1}, \ldots ,\Delta {{\boldsymbol{t}}_M}]^{\text{T}}}$，$\Delta {\boldsymbol{v}} = {[\Delta {{\boldsymbol{v}}_1}, \ldots ,\Delta {{\boldsymbol{v}}_M}]^{\text{T}}}$. $\Delta {\boldsymbol{t}}$和$\Delta {\boldsymbol{v}}$分别服从均值为0，方差分别为$ {\sigma _t} $和$ {\sigma _v} $的高斯分布，协方差矩阵分别为${{\boldsymbol{Q}}_{\boldsymbol{t}}}$和${{\boldsymbol{Q}}_{\boldsymbol{v}}}$.

考虑到分层效应^[17]，在表达式(1)中$ \tau(\boldsymbol{u}^{\mathrm{o}},\boldsymbol{t}_i^{\mathrm{o}}) $和${g_i}\left\| {{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}} \right\| = {\boldsymbol{v}}_i^{\text{T}}{{\boldsymbol{t}}_i} - {\boldsymbol{v}}_i^{\text{T}}{{\boldsymbol{u}}^{\mathrm{o}}}$分别表示为

式中：a为声速梯度；当传感器i发射探测信号后，$ \theta _i^{\text{S}} $和$ \theta _i^{\text{R}} $分别表示传输声线在传感器i(${\boldsymbol{t}}_i^{\mathrm{o}}$)和目标处的射线角度；当信号从目标反射回来时，$ \tilde \theta _i^{\text{S}} $和$ \tilde \theta _i^{\text{R}} $分别表示反射的声线在传感器i(${\boldsymbol{s}}_i^{\mathrm{o}}$)和目标处的射线角度.

式(3)中$\dot \tau ({\boldsymbol{u}}_i^{\mathrm{o}},{\boldsymbol{t}}_i^{\mathrm{o}})$和$\dot \tau ({{\boldsymbol{u}}^{\mathrm{o}}},{\boldsymbol{s}}_i^{\mathrm{o}})$的具体表达式为

其中：$ \gamma _i^\beta = \displaystyle\frac{{ - {v_{i,z}}r_i^2{\text{ + }}({z^{\mathrm{o}}} - z_{\text{A}}^i){{\dot r}_i}}}{{d_i^2{r_i}}} $；$ \gamma _i^\alpha = - \displaystyle\frac{{0.5a{q_i}{{\dot r}_i} + {{(0.5a)}^2}r_i^2{v_{i,z}}}}{{(p_i^2 + q_i^2){r_i}}} $；$\displaystyle\tilde \gamma _i^\beta = \frac{ - {{\tilde v}_{i,z}}\tilde r_i^2{{ + }}({z^{\mathrm{o}}} - z_{\text{B}}^i)\tilde{\dot{r}}}{\tilde {d}_i^2{\tilde {r}_i}} $；

$ \displaystyle\tilde \gamma _i^\alpha = - \frac{{0.5a{{\tilde q}_i}{{\tilde \dot r}_i} + {{(0.5a)}^2}\tilde r_i^2{{\tilde v}_{i,z}}}}{{(\tilde p_i^2 + \tilde q_i^2){r_i}}} $；

$ {r_i} = \sqrt {{{({x^{\mathrm{o}}} - x_{\text{A}}^i)}^2} + {{({y^{\mathrm{o}}} - y_{\text{A}}^i)}^2}} $；$ {\tilde r_i} = \sqrt {{{({x_i} - x_{\text{B}}^i)}^2} + {{({y_i} - y_{\text{B}}^i)}^2}} $；$ {\dot r_i} = ({x^{\mathrm{o}}} - x_{\text{A}}^i){v_{i,x}} + ({y^{\mathrm{o}}} - y_{\text{A}}^i){v_{i,y}} $；$ \tilde{\dot{r}} = ({x^{\mathrm{o}}} - x_{\text{B}}^i){\tilde v_{i,x}} + ({y^{\mathrm{o}}} - y_{\text{B}}^i){\tilde v_{i,y}} $；$ {d_i} = \sqrt {r_i^2 + {{({z^{\mathrm{o}}} - z_{\text{A}}^i)}^2}} $；$ {\tilde d_i} = \sqrt {\tilde r_i^2 + {{({z^{\mathrm{o}}} - z_{\text{B}}^i)}^2}} $；$ {p_i} = 0.5a{r_i} $；$ {q_i} = b + 0.5a({z^{\mathrm{o}}} + z_{\text{A}}^i) $；$ {\widetilde p_i} = 0.5a{\widetilde r_i} $；$ {\tilde q_i} = b + 0.5a({z^{\mathrm{o}}} + z_{\text{B}}^i) $；b为水面声速；$[{\tilde v_{i,x}},{\tilde v_{i,y}},{\tilde v_{i,z}}] = (1 + \dot \tau ({{\boldsymbol{u}}^{\mathrm{o}}},{\boldsymbol{t}}_i^{\mathrm{o}}){{ + }}\dot \tau ({{\boldsymbol{u}}^{\mathrm{o}}},{\boldsymbol{s}}_i^{\mathrm{o}}))[{v_{i,x}},{v_{i,y}},{v_{i,z}}]$.

2.   水下定位方案

本文所提的水下定位方案主要包含两大部分：首先是建立声线直线传输假设下的测量模型，并通过CFS方法得到初始解，其次在考虑声线弯曲和锚节点运动效应的条件下，建立最大似然代价函数，并运用高斯-牛顿迭代法求解精确解. 算法流程如图2所示.

图  2  水下定位流程图

下载: 全尺寸图片 幻灯片

2.1   高斯-牛顿迭代法

时延和多普勒频移的测量模型具有递归和非线性特性从而导致该模型下的最大似然函数极其复杂. 高斯-牛顿法具有计算效率高、局部收敛快的优点. 因此，将使用高斯-牛顿迭代法来实现目标的位置估计. 根据时延和多普勒频移的观测量模型，建立关于未知量$ {\boldsymbol{u}} $的MLE，如下所示：

式中：${\boldsymbol{h}}({\boldsymbol{u}}) = {[{\boldsymbol{\tau }}{({\boldsymbol{u}})^{\text{T}}},{\boldsymbol{f}}{({\boldsymbol{u}})^{\text{T}}}]^{\text{T}}}$，${\boldsymbol{\tau }}({\boldsymbol{u}})$和${\boldsymbol{f}}({\boldsymbol{u}})$表示关于未知向量${\boldsymbol{u}}$的时延和多普勒频率的参数形式，如式(7)和(8)所示. ${\boldsymbol{\alpha }}$为测量值，如式(5)所示.

使用高斯-牛顿迭代法解算出目标位置，根据初始解${{\boldsymbol{u}}^{(0)}}$，${\boldsymbol{h}}({\boldsymbol{u}})$的线性表达式可以表示为

式中：${\boldsymbol{G}}_{\tau ,f}^{(0)}$表示迭代梯度矩阵，其具体表达式为

将式(13)带入式(12)中，并使其最小化，从而可以得到$ {\boldsymbol{u}} $的解. 因此，其迭代求解过程如下：

式中：$ k = 0,1,2, \cdots $，为迭代次数；$ {{\boldsymbol{u}}^{(k)}} $为第k次迭代后得到的解；$ {\boldsymbol{h}}({{\boldsymbol{u}}^{(k)}}) = {[{\mathbf{\tau }}{({{\boldsymbol{u}}^{(k)}})^{\text{T}}},{\boldsymbol{f}}{({{\boldsymbol{u}}^{(k)}})^{\text{T}}}]^{\text{T}}} $为$ {\boldsymbol{u}} = {{\boldsymbol{u}}^{(k)}} $时时延和多普勒频率的值；$ {\boldsymbol{G}}_{\tau ,f}^{(k)} $为在$ {\boldsymbol{u}} = {{\boldsymbol{u}}^{(k)}} $时的梯度矩阵. 迭代的终止条件设置为$ \left\| {{{\boldsymbol{u}}^{(k + 1)}} - {{\boldsymbol{u}}^{(k)}}} \right\| < \xi $或$ k > {N_{{\text{max}}}} $.

2.2   基于CFS的初始解

为了保证上述迭代实现全局收敛，需要提供一个良好的初始解${{\boldsymbol{u}}^{(0)}}$. 建立声线直线传输假设下时延和多普勒频率测量模型，并利用最小二乘法解算目标的初始位置，从而获得基于CFS的初始解${{\boldsymbol{u}}^{(0)}}$. 在声线直线传输的假设下，时延信息模型如下：

式中：c为平均声速. 对式(15)的两边同时平方，并整理可得

将$ {\tau _i} = \tau _i^{\mathrm{o}} + \Delta {\tau _i} $，${{\boldsymbol{t}}_i} = {\boldsymbol{t}}_i^{\mathrm{o}} + \Delta {{\boldsymbol{t}}_i}$，${{\boldsymbol{v}}_i} = {\boldsymbol{v}}_i^{\mathrm{o}} + \Delta {{\boldsymbol{v}}_i}$带入式(16)，忽略二阶误差，可以得到

式中：$ c_{i}=a_{i} \tau_{i}-b_{i} $；${a_i} = {{(c - {{{{\left\| {{{\boldsymbol{v}}_i}} \right\|}^2}} \mathord{\left/ {\vphantom {{{{\left\| {{{\boldsymbol{v}}_i}} \right\|}^2}} c}} \right. } c})} \mathord{\left/ {\vphantom {{(c - {{{{\left\| {{{\boldsymbol{v}}_i}} \right\|}^2}} \mathord{\left/ {\vphantom {{{{\left\| {{{\boldsymbol{v}}_i}} \right\|}^2}} c}} \right. } c})} 2}} \right. } 2}$；$b_i = {\boldsymbol{v}_i^{{\mathrm{T}}} \boldsymbol{t}_i}/{c}$；$\xi _i^\tau = {a_i}\Delta {\tau _i} - {\boldsymbol{c}}_{{\boldsymbol{t}},i}^{\text{T}}\Delta {\boldsymbol{t}} - {\boldsymbol{c}}_{{\boldsymbol{v}},i}^{\text{T}}\Delta {\boldsymbol{v}}$；${{\boldsymbol{c}}_{{\boldsymbol{t}},i}} = {{{{\boldsymbol{v}}_i}} \mathord{\left/ {\vphantom {{{{\boldsymbol{v}}_i}} c}} \right. } c} - {{\boldsymbol{\rho }}_{{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}}}$；$c_{v, I}=(\boldsymbol{v}_{I} \tau_{I}+t_{i}- \boldsymbol{u}^{0}) / c$；${\boldsymbol{\rho}} _{\boldsymbol{u f}^{\rho}-\boldsymbol{t}_{i}}=\left(\boldsymbol{u}^{0}-\boldsymbol{t}_{i}^{0}\right) /\left\|\boldsymbol{u}^{0}-\boldsymbol{t}_{i}^{0}\right\|$.

将式(17)的两边同时平方，整理可得

式中：${{\boldsymbol{\gamma }}^{\mathrm{o}}} = {[{{\boldsymbol{u}}^{{\mathrm{o}}{\text{T}}}} \odot {{\boldsymbol{u}}^{{\mathrm{o}}{\text{T}}}},{x^{\mathrm{o}}}{y^{\mathrm{o}}},{x^{\mathrm{o}}}{z^{\mathrm{o}}},{y^{\mathrm{o}}}{z^{\mathrm{o}}}]^{\text{T}}}$；${{\boldsymbol{q}}_{i,\tau }} =$$ \odot {\boldsymbol{v}}_i^{\text{T}},2{v_{i,x}}{v_{i,y}}, 2{v_{i,x}} {v_{i,z}},2{v_{i,y}}{v_{i,z}}]^{\text{T}} / {{c^2}} - {[1,1,1,0,0,0]^{\text{T}}}$.

根据$ {{\mathbf{\gamma }}^{\mathrm{o}}} $，定义新的未知向量，${{\textit{φ}}}_1^{\mathrm{o}} = {[{{\boldsymbol{u}}^{{\mathrm{o}}{\text{T}}}},{{\boldsymbol{\gamma }}^{{\mathrm{o}}{\text{T}}}}]^{\text{T}}}$，因此，式(18)可以表示为

式中：${{\boldsymbol{h}}_\tau }(i) = c_i^2 - {\boldsymbol{t}}_i^{\text{T}}{{\boldsymbol{t}}_i}$；${{\boldsymbol{A}}_\tau }(i,:) = [ - 2{[{{\boldsymbol{t}}_i} + {{{c_i}{{\boldsymbol{v}}_i}} \mathord{\left/ {\vphantom {{{c_i}{{\boldsymbol{v}}_i}} c}} \right. } c}]^{\text{T}}},$$ - {\boldsymbol{q}}_{i,\tau }^{\text{T}}] $；${{\boldsymbol{B}}_\tau }$为对角矩阵，${{\boldsymbol{B}}_\tau }(i,i) = 2\left\| {{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}} \right\|$； ${{\boldsymbol{\xi }}^\tau } = {{\boldsymbol{C}}_\tau }\Delta {\boldsymbol{\tau }} + {{\boldsymbol{C}}_{\boldsymbol{t}}}\Delta {\boldsymbol{t}} + {{\boldsymbol{C}}_{\boldsymbol{v}}}\Delta {\boldsymbol{v}}$，${{\boldsymbol{C}}_\tau }$为对角矩阵，${{\boldsymbol{C}}_\tau }(i,i) = {c_{\tau ,i}}$；${{\boldsymbol{C}}_{\boldsymbol{t}}}$和${{\boldsymbol{C}}_{\boldsymbol{v}}}$均为块对角矩阵，${{\boldsymbol{C}}_{\boldsymbol{t}}}(3i - 2:3i,i) = {\boldsymbol{c}}_{{\boldsymbol{t}},i}^{\text{T}}$；${{\boldsymbol{C}}_{\boldsymbol{v}}}(3i - 2: 3i,i) = {\boldsymbol{c}}_{{\boldsymbol{v}},i}^{\text{T}}$.

对式(16)在时间上求导，得到多普勒频移信息的观测方程，如下所示：

式中：${g_i} = {{{a_i}f_i^{\mathrm{o}}}/ {f_i^{\mathrm{c}}}} - {d_i}$；${d_i} = {{{{\left\| {{{\boldsymbol{v}}_i}} \right\|}^2}} / c}$.

将$ {f_i} = f_i^{\mathrm{o}} + \Delta {f_i} $，${{\boldsymbol{t}}_i} = {\boldsymbol{t}}_i^{\mathrm{o}} + \Delta {{\boldsymbol{t}}_i}$，${{\boldsymbol{v}}_i} = {\boldsymbol{v}}_i^{\mathrm{o}} + \Delta {{\boldsymbol{v}}_i}$带入，可以得到

式中：$\xi _i^f = {d_{f,i}}\Delta {f_i} - {\boldsymbol{d}}_{{\boldsymbol{t}},i}^{\text{T}}\Delta {\boldsymbol{t}} - {\boldsymbol{d}}_{{\boldsymbol{v}},i}^{\text{T}}\Delta {\boldsymbol{v}}$；${d_{f,i}} = {{{a_i}\left\| {{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}} \right\|} / {f_i^{\text{c}}}}$；${{\boldsymbol{d}}_{{\boldsymbol{t}},i}} = {g_i}{{\boldsymbol{\rho }}_{{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}}} + {{\boldsymbol{v}}_i}$；${{\boldsymbol{d}}_{{\boldsymbol{v}},i}} = [({{{f_i}} / {cf_i^{\mathrm{c}}}} + {2 / c}){{\boldsymbol{v}}_i} - {{\boldsymbol{\rho }}_{{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}}}]\left\| {{{\boldsymbol{u}}^{\mathrm{o}}} - {{\boldsymbol{t}}_i}} \right\|$；${{\boldsymbol{q}}_{i,f}} = {[{\boldsymbol{v}}_i^{\text{T}} \odot {\boldsymbol{v}}_i^{\text{T}},2{v_{i,x}}{v_{i,y}},2{v_{i,x}}{v_{i,z}},2{v_{i,y}}{v_{i,z}}]^{\text{T}}} - g_i^2{\left[ {1,1,1,0,0,0} \right]^{\text{T}}}$.

式(21)的向量表达式如下所示：

式中：${{\boldsymbol{h}}_f}(i) = g_i^2{\boldsymbol{t}}_i^{\text{T}}{{\boldsymbol{t}}_i} - {({\boldsymbol{t}}_i^{\text{T}}{{\boldsymbol{v}}_i})^2}$；${{\boldsymbol{A}}_f}(i,:) = [2(g_i^2{{\boldsymbol{t}}_i} - {\boldsymbol{v}}_i^{\text{T}}{{\boldsymbol{t}}_i}{{\boldsymbol{v}}_i}), {\boldsymbol{q}}_{i,f}^{\text{T}}]$；${{\boldsymbol{B}}_f}(i,i) = 2({\boldsymbol{v}}_i^{\text{T}}{{\boldsymbol{u}}^o} - {\boldsymbol{v}}_i^{\text{T}}{{\boldsymbol{t}}_i})$.${{\boldsymbol{\xi }}^f} = {{\boldsymbol{D}}_f}\Delta {\boldsymbol{f}} + {{\boldsymbol{D}}_{\boldsymbol{t}}}\Delta {\boldsymbol{t}} + {{\boldsymbol{D}}_{\boldsymbol{v}}}\Delta {\boldsymbol{v}}$，${{\boldsymbol{D}}_f}$为对角矩阵，具体地，${{\boldsymbol{D}}_f}(i,i) = {d_{f,i}}$；${{\boldsymbol{D}}_{\boldsymbol{t}}}$和${{\boldsymbol{D}}_{\boldsymbol{v}}}$均为块对角矩阵，${{\boldsymbol{D}}_{\boldsymbol{t}}}(3i - 2:3i,i) = {\boldsymbol{d}}_{{\boldsymbol{t}},i}^{\text{T}}$；${{\boldsymbol{D}}_{\boldsymbol{v}}}(3i - 2:3i,i) = {\boldsymbol{d}}_{{\boldsymbol{v}},i}^{\text{T}}$.

利用时延和多普勒频率的测量信息联合求解未知向量$ {\boldsymbol{u}} $，将式(19)和式(22)组合得到

式中：${{\boldsymbol{A}}_1} = [{{\boldsymbol{A}}_\tau },{{\boldsymbol{A}}_f}]$；${{\boldsymbol{h}}_1} = [{{\boldsymbol{h}}_\tau },{{\boldsymbol{h}}_f}]$；${{\boldsymbol{B}}_1} = {{\mathrm{blkdiag}}} \{ {{\boldsymbol{B}}_\tau },{{\boldsymbol{B}}_f}\}$；$ {\boldsymbol{\xi }} = [{{\boldsymbol{\xi }}^\tau },{{\boldsymbol{\xi }}^f}] $. 因此，可以得到未知向量${{{\textit{φ}}}_1}$的解为

式中：权重向量为${{\boldsymbol{W}}_1} = {({{\boldsymbol{B}}_1}{{\boldsymbol{Q}}_1}{\boldsymbol{B}}_1^{\text{T}})^{ - 1}}$，${{\boldsymbol{Q}}_1} = \boldsymbol{E}[{\boldsymbol{\xi }}{{\boldsymbol{\xi }}^{\text{T}}}]$. 因此，可以得到基于CFS的初始解为${{\boldsymbol{u}}^{(0)}} = {{{\textit{φ}}}_1}(1:3)$.

3.   水下定位方案

仿真环境为1.5 km×1.5 km×0.5 km的三维水域，其中存在1个待定位的水下目标和24个水下传感器. 选定的目标节点位置为${{\boldsymbol{u}}^{\mathrm{o}}} = {[500,800,80]^{\text{T}}}$，水下传感器的速度设定为$\left\| {{{\boldsymbol{v}}_i}} \right\| = 10{\text{ m/s}}$，其位置坐标随机分布于目标水域中. 为了防止传感器距离太近导致定位精度退化，其位置必须满足$\left\| {{\boldsymbol{s}}_i^{\mathrm{o}} - {\boldsymbol{s}}_j^{\mathrm{o}}} \right\| \geqslant 400{\text{ m}}$，$ i,j = 1,2, \cdots ,M $，且$ i < j $. 为了防止信号的频率干扰，传感器发出的探测信号的载频设置为$f_{i + 1}^{\text{c}} - f_i^{\text{c}} = 0.25{\text{ kHz}}$，$ i = 1,2, \cdots ,M - 1 $.水下声速模型中声速梯度a=0.1，水面声速为b=1 450 m/s.

时延信息和多普勒频率信息叠加的独立高斯白噪声的方差分别为$ \sigma _\tau ^2 $和$ \sigma _f^2 $，设置为$ \sigma _f^2{\text{ = }}100\sigma _\tau ^2 $.根据协方差矩阵${{\boldsymbol{Q}}_{\boldsymbol{\beta }}} = \sigma _{\boldsymbol{\beta }}^2 \cdot {\text{blkdiag}}\left\{ {{{\boldsymbol{I}}_{3M}},0.1{{\boldsymbol{I}}_{3M}}} \right\}$来仿真具有偏差的水下传感器位置和速度，如公式(6)所示. 通过1000次蒙特卡洛仿真，统计得到不同$ \sigma _\tau ^2 $情况下水下目标位置坐标的均方根误差(root mean square error, RMSE). 在仿真中，为了验证所提算法的有效性，与现有方法进行了对比，主要包括：

1) CFS-MLE表示本文所提出的水下定位算法，即建立MLE，然后使用高斯-牛顿迭代法寻找最优解. 然后，基于CFS的目标位置初始解，来保证上述迭代算法收敛于最优值.

2) SDP-MLE表示建立MLE，然后使用高斯-牛顿迭代法寻找最优解. 然后，基于SDP的目标位置初始解，来保证上述迭代算法收敛于最优值.

3) CFS表示在不补偿分层效应的情况下用于目标定位的CFS解决方案，如文献[11]所述.

4) SDP表示在不补偿分层效应的情况下用于目标定位的SDP解决方案，如文献[13]所述.

5) MLE-随机表示建立MLE，然后使用高斯-牛顿迭代法寻找最优解，其初始解随机取值.

6) CFS-LM表示在构建MLE模型的基础上，采用Levenberg-Marquardt(LM)算法求解最优解，该方法的初始值来自本文所提出的CFS方法.

图3给出了不同噪声值下的定位性能. 当具有小噪声水平时，CFS-MLE和SDP-MLE方法的RMSE都达到了CRLB性能. 虽然当噪声水平增加到中等或者较大水平时，CFS-MLE的性能会逐渐远离CRLB精度限，但它可以提供比其他方法更高的精度. 另外，无论测量噪声水平如何，由于CFS和SDP方法忽略了水下的分层效应，造成了测量模型偏差，均无法达到CRLB性能. MLE-随机方法的估计精度较差，因为随机的初始解可能会导致定位结果收敛到局部最优值. CFS-LM方法的性能与本文所提算法较为接近，但该方法需要通过多次迭代来调整参数，其计算复杂度相较于本文方法更高.

图  3  不同定位方法下的性能对比

下载: 全尺寸图片 幻灯片

表1显示了不同定位方法在MATLAB中的平均处理时间. CFS-MLE方法由于考虑了分层效应，因此其处理时间约是CFS方法的2.5倍，但整体时间开销相对于其他方法较低. SDP-MLE方法略高于SDP，总体来说两者的处理时间都比其他方法高很多，在耗时方面不具有优势. 此外，测量噪声水平的不同对这些方法的平均处理时间几乎没有影响.

表  1  不同定位方法下的平均处理时间 ms

定位方法	噪声水平
	0 dB	5 dB
CFS-MLE	0.296	0.297
SDP-MLE	0.679	0.675
CFS	0.083	0.085
SDP	0.634	0.647
MLE-随机	0.213	0.216
CFS-LM	0.352	0.353

下载: 导出CSV 

| 显示表格

图4给出了忽略信号观测期间内传感器运动效应下的定位性能，其中${\sigma _{\boldsymbol{\beta }}}{{ = }}0.5{\text{ m}}$. 可以看出，由于忽略传感器的运动效应会导致时延和多普勒频移观测方程的偏差，从而导致定位性能下降. 因此，在水下定位中，考虑传感器运动效应非常有必要.

图5给出了忽略传感器位置和速度偏差时的定位性能，其中$\left\| {{{\boldsymbol{v}}_i}} \right\| = 10{\text{ m/s}}$和${\sigma _\tau }{{ = }}1{\text{ ms}}$. 从结果可以看出，忽略传感器的位置和速度误差会导致定位性能的下降，且随着$ {\sigma _{\boldsymbol{\beta }}} $的增大，其偏差会增大.

图  4  忽略传感器运动效应时的定位性能

下载: 全尺寸图片 幻灯片

图  5  忽略传感器位置和速度误差时的定位性能

下载: 全尺寸图片 幻灯片

4.   结束语

本文提出了一种基于时延-多普勒频移的水下目标定位方法，该方法中联合考虑了水下环境的分层效应，传感器在信号观测期间的运动效应，以及传感器的位置和速度偏差. 在此条件下，建立时延和多普勒频移的测量模型，随后构建了基于时延-多普勒频移的MLE，并利用高斯-牛顿迭代法来估计目标的位置. 为了保证迭代方法能够实现全局收敛，提出了在声线直线传输条件下的数学CFS，并作为高斯-牛顿迭代的初始解. 仿真结果显示，该方法具有良好的定位性能，并且在水下定位系统中忽视传感器的运动效应、位置和速度误差将会导致测量模型的偏差，从而降低定位精度. 此外，解决锚节点的非线性运动特性以及多元传感器融合等问题将作为下一步的主要研究方向.