一种结合声速反演的水声传感网络节点定位误差修正方法与流程

本发明属于水下定位技术领域，具体是涉及水下分布式定位和声速反演技术领域，涉及一种结合声速反演的水声传感网络节点定位误差修正方法。

背景技术：

在大多数水声传感网络(underwateracousticsensornetwork，uasn)的应用中，定位都是一个基础且重要的任务。例如,在位置相关的数据采集任务中，如果缺失了位置信息，数据将变得无用。此外，节点的位置信息有助于增强其他水下任务的完成，如用于辅助路由算法的设计，帮助提高网络的通信效率，以及降低节点间通信能量消耗等。

水声传感网络中一种最基本的定位方法是三边测量定位技术。为了确定目标节点的位置，我们需要一些位置信息已知的网络节点作为参考节点，并测量目标节点到这些参考节点的距离。在一个三维的坐标系中，如果一个点到其他三个参考点的距离已知，那么这个点的位置坐标可以通过求三个球面的相交位置得出。在水下环境中，节点的深度信息可以比较方便的通过压力传感器测得，因此基于时间测量的定位方法，其性能的好坏很大程度上取决于节点间距离估计。

节点间距离一般由传播时间结合声速信息转化得到。常见的两种距离估计方法：

(1)传播时间乘以等效声速，这是最简单也是最常用方法。一般地，如果网络节点都位于同一深度，则取目标位置的声速值作为等效声速，否则使用垂直剖面上的平均声速值。由于假设信号沿直线传播，忽略了水体的分层效应，因此采用这种方法的距离估计存在一定的偏差，这个偏差是距离有关的，并且受到声速剖面的影响。随着距离的增大，实际声线路径逐渐偏离直线，变得越来越弯曲，直线传播假设带来的偏差会越来越大。

(2)利用射线声学知识求解。由折射定律可以知道声线传播时间和声速剖面的关系为：

给定声速剖面c(z)与节点间信号传播时间，我们可以先求出k0，然后再把求得的k0代入下式便可以求得参考节点和目标节点间的距离：

其中，

该方法需要确知的声速剖面，并且假设起止深度zs和zd间任一深度声线只经过了一次。

由于水下介质的不均匀与海底海面对信号的反射，声线的传播路径并不是直线，会发生弯曲，具体的声线传播路径受到声速分布的影响。此外，海洋环境的动态变化，特别是声速剖面的扰动，会造成声线路径和信号传播时间的变化，影响定位结果的稳定性。

技术实现要素：

本发明的目的是针对现有技术的不足，提供一种结合声速反演的水声传感网络节点定位误差修正方法，利用水声传感节点间的信号测量对网络覆盖区域的平均声速进行估计，并基于射线跟踪模型完成传统定位方法基础上的再定位，能够对传统模型的固有偏差和环境动态变化引起的定位不确定进行修正，使水声传感网络节点定位不受声速变化的影响。

本发明的具体步骤是：

步骤(1).输入待定位节点处到达时间差观测向量δt、水声传感网络(underwateracousticsensornetwork，uasn)节点位置信息、m对水声传感网络节点间信号传播时间组成的向量tpert＝[tp1,…,tpm,…,tpm]^t、参考声速剖面c0(z)、经验正交函数(empiricalorthogonalfunction，eof)fk(z)、到达时间测量协方差矩阵r∈；

步骤(2).声速剖面反演：根据水声传感网络节点间信号传播时间tpert，反演网络覆盖区域平均声速剖面c(z)；具体方法如下：

(2-1).求解经验正交函数对应系数：根据水声传感网络节点间传播时间扰动和经验正交函数对应系数之间的线性关系y＝eα+v，求解得到经验正交函数对应系数向量α＝[α1,…,αk,…,αk]^t；αk为对应第k阶的经验正交函数对应系数，k＝1,2,…,k，k为选取的经验正交函数的数量，上标t表示转置；

y为m个传播时间扰动测量，y＝[δτ1,…,δτm,…,δτm]^t，其中δτm＝tpm-trm，t表示转置；v为模型噪声；观测矩阵e的第(i,j)项元素表示第j阶经验正交函数在第i层对应声线路径上对时延扰动的影响，其中fj(·)为第j阶经验正交函数，γi(-)为声线在第i层的路径，c0(·)为参考声速剖面；根据声学计算模得到参考声速剖面c0(z)下对应的信号传播时间tref＝[tr1,…,trm,…,trm],所述的声学计算模型包括bellhop、kraken等；

(2-2).声速剖面估计：

声速剖面估计值即为网络覆盖区域平均声速剖面c(z)：c0(z)为参考声速剖面。

步骤(3).寻找初始位置估计：

将待定位节点处的到达时间差(timedifferenceofarrival，tdoa)测量输入常规定位方法，如双曲定位方法，得到初步定位结果x⁽⁰⁾。

步骤(4).定位误差修正：根据网络覆盖区域平均声速剖面c(z)、到达时间差观测向量δt、初步定位结果x⁽⁰⁾，对定位结果进行修正；具体如下：

(4-1).对声速剖面进行分段线性近似处理：根据c(z)的变化特性，对c(z)进行分段近似处理，分段点记为zp，求出每段的声速梯度a和起始深度处的声速b，得到声速剖面的分段线性表征cp(z)＝apz+bp；其中，cp(z)为第p段声速剖面的线性表征，ap为第p段声速剖面的声速梯度，bp为第p段起始深度处的声速。

(4-2).设置声线路径解析点：根据(4-1)得到的分段线性近似处理结果，设置水声传感网络节点和当前搜索位置x^(k)之间声线路径的解析点(初始时刻k＝0)；由于第p段声速剖面的声速梯度为常数ap，根据斯涅尔定理，对应的声线路径为一圆弧，圆弧由起点和终点坐标确定。在水平距离深度二维平面内，第p声线路径的起点记为终点记为其中，字母s表示start，字母e表示end，和为第p段声线起点和终点对应的水平距离，为待求变量；和为第p段声线起点和终点对应的深度，和(4-1)中的分段点zp对应。

(4-3).建立声线路径方程：根据声线路径的连续性，可得边界条件即第p段声线的终点处的掠射角和第p+1段声线起点处的掠射角相等。对等式两边同时取tan(·)，结合第p段声线起点和终点处掠射角之间的关系其中，αp是起点和终点连线与水平方向的夹角，βp是第p段声线掠射角与水平夹角αp的差值；由三角函数的性质，整理得到声线路径方程其中xp＝tanβp，kp＝tanαp·tanβp。

(4-4).求解声线路径方程：根据关系式求出kp的值，代入(4-3)中的声线路径方程后，解方程求出xp。

(4-5).确定声线路径：由关系式求出(4-4)中每个xp对应的和并结合(4-2)路径解析点的设置，得到第p段声线的起点和终点所有的声线路径解析点一起确定了完整的声线路径。

(4-6).计算模型到达时间差观测：根据(4-4)求得的kp、xp，由βp＝arctanxp，以及得到第p段声线起点处的掠射角和终点处的掠射角代入关系式得到当前搜索位置x^(k)和水声传感网络节点之间的信号传播时间向量h(x^(k))；选择h(x^(k))的第一个元素对应的水声传感网络节点作为参考信标，h(x^(k))的其他元素与之做差，得到模型到达时间差观测h′(x^(k))；

(4-7).计算模型在当前迭代位置处的梯度：边界条件以及信号传播时间表达式关于迭代位置x^(k)＝(r^t,z^t)求偏导，其中，r^t,z^t的上标t表示target，将得到的偏导关系(4-1)求得的第p段声速剖面的声速梯度ap，(4-6)求得的第p段声线起点处掠射角和终点处掠射角代入传播时间关于搜索位置的梯度关系式和得到传播时间关于当前搜索位置x^(k)的梯度类似(4-6)，选择第一个元素对应的水声传感网络节点作为参考信标，其他元素与第一个元素做差得到迭代位置x^(k)处对应的模型梯度

(4-8).位置估计更新：将实际到达时间差测量δt、模型tdoa观测h'(x^(k))、模型梯度输入得到第k+1步的迭代估计x^(k+1)；

(4-9).输出修正后的位置估计：重复(4-2)～(4-8)，直到算法收敛，输出收敛后的位置估计x^(k)。

本发明通过对声速剖面的线性分段近似，建立了节点间声线路径的解析表达式，并将任意两点间声线路径跟踪转化为多项式求根问题。在此基础上，提出了一种利用tdoa测量的水下节点定位误差修正方法。本发明具有以下有益技术效果：

(1)本发明方法能够通过水声传感网络对所覆盖区域的平均声速进行估计，为定位模型更新声速剖面信息，修正环境声速变化造成的定位误差；

(2)本发明方法针对不同海域、不同深度的目标都能进行定位误差修正；

(3)本发明方法不需要改变已有的水声传感网络拓扑，只需要选择合适的节点为声速估计提供测量信息。

附图说明

图1反演节点拓扑与定位信标拓扑关系；

图2结合声速反演的定位误差修正原理图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，一种结合声速反演的水声传感网络节点定位误差修正方法的反演和定位节点拓扑示意图。这里以两个正六边形组成的水声传感网络为例，图中每个小正三角形区域对应一个定位覆盖区域。本发明的主要思路是用大拓扑中的几个节点(图中节点a、b、c、d围成的实线四边形)进行声速反演，然后用于修正小正三角形拓扑中(图中节点a、e、f围成的实线三角形)的定位。整个流程分为两个部分，如图2所示，第一部分是基于反演拓扑节点间的传播时间测量，对声速剖面信息反演。流程的第二部分是将估计得到的声速剖面和初始定位结果输入到基于声线跟踪的定位模型，对目标位置进行高精度搜索，算法收敛后输出定位修正结果。声速剖面反演是利用水声传感网络节点间信号传播时间扰动，基于微扰法，对信号经过水体的声速扰动进行测量，具体方法是：

(1.1)经验正交函数系数求解：

声速剖面由eof形式表征，可降低待求参数的维度，eof为经验正交函数(empiricalorthogonalfunction)，得到网络节点间传播时间扰动和经验正交函数对应系数之间的线性关系为：

y＝eα+v(s1)；

其中，y为m个传播时间扰动测量，y＝[δτ1,δτ2,…,δτm]^t，t表示转置；α为待求参数，α＝[α1,α2,…,αk]^t,α中的元素对应前k阶eof系数；v为模型噪声；观测矩阵e的第(i,j)项元素表示第j阶经验正交函数对应系数在第i层对应声线路径上对时延扰动的影响，其中fj(·)为第j阶经验正交函数，γi(-)为声线在第i层的路径，c0(·)为参考声速剖面；

(1.2)声速剖面估计；

由(s1)求得经验正交函数对应系数后，声速剖面估计为：(s2)；c0(z)为参考声速剖面，由历史声速剖面平均得到。

定位误差修正是基于反演得到的声速剖面，结合声线跟踪模型，对定位结果进行修正，具体方法是：

(2.1)声线路径求解；

水声传感网络定位修正的关键是对声线路径进行建模。本发明基于梯度分层模型求解节点间声线路径的解析表达式。

梯度分层模型将声速剖面近似为等梯度的分段线性函数，即对(s2)的声速剖面线性分段得到：cp(z)＝apz+bp(s3)；其中，ap为第p段声速剖面的声速梯度，bp为该段起始深度处的声速。

采用该模型以简化声线路径的计算。在声线路径求解中，将多层水体中的声线路径，拆分成一个个相连的单层路径，并加上边界条件确定，边界条件为：p＝1,2,…,p-1(s4)；即第p段声线的终点处的掠射角和第p+1段声线起点处的掠射角相等。

由斯涅尔定理以及声线路径的几何关系，第p段声线由其起点和终点以及水体的声速确定，即：

tanβp＝xp,tanαp＝yp(s55)，

其中，xp、yp是辅助变量；和分别是声线在起点和终点处的掠射角；αp是起点和终点连线与水平方向的夹角，βp是第p段声线掠射角与水平夹角αp的差值，tp是第p段声线对应的信号传播时间。

本发明对声线路径进行分段建模。由于第p段声速剖面的声速梯度为常数ap，根据斯涅尔定理，对应得声线路径为一圆弧，圆弧可由起点和终点坐标确定。对应二维笛卡尔坐标平面下圆弧的起点(上标s表示start),则对应圆弧的终点(上标e表示end)。

根据(s3)中声速剖面的线性分段，第p段声线对应的起点和终点深度信息和是已知的；第p段声线由起点和终点水平位置和确定。将(s56)代入式(s4)，并两边同时求tan(·)，由三角函数和差化积的性质，结合(s55)，整理后得到(s4)边界条件对应的多项式方程：

该式是关于xp或者rp的高阶多项式。

求解上述高阶多项式，由(s57)得到p段声线总的信号传播时间为：

当声线发生海底海面反射时，边界方程修改为：

而(s6)对应高阶多项式的阶数不发生改变。

(2.2)迭代算法修正定位误差；

由于深度信息均已知，(s6)的p-1个边界条件和(s7)联立，组成包含p个描述声线路径的解析点(由于深度已知zp，解析点位置由水平距离rp确定)和一个描述信号总传播时间t的包含p个方程的多项式方程组，将该关系记为函数h(·)，并取rp,p＝1,2,…,p作为函数的自变量，t作为因变量，则函数h(·)建立了总传播时间和声线路径之间的解析关系，该声线路径由rp唯一确定。

本方明的定位修正算法是基于声线路径跟踪的，定位场景为待定位节点与水声传感网络信标之间不同步的场景，也即基于待定位节点处到达时间差(timedifferenceofarrival,tdoa)测量的定位场景。

设水声传感网络中共有n个信标节点，tdoa测量模型为：δt＝h'(x)+∈；其中，δt为tdoa观测向量，h'(·)为tdoa关于位置x＝[x,y,z]的映射关系，根据h(·)函数做差得到，∈为模型噪声。假设待定位节点处水声传感网络各节点发射信号的测量噪声均为高斯白噪声，则协方差矩阵为：为水声传感网络各信标信号的测量噪声能量。

采用基于高斯牛顿迭代算法对定位结果进行修正，收敛后输出定位修正结果。其中，迭代过程中的一步迭代为：其中，x^(k+1)为第k步迭代结果，迭代起点x⁽⁰⁾为初始定位结果；h'(x^(k))为当前搜索位置x^(k)对应的模型tdoa测量；δt为实际到达时间差测量；r∈为测量误差协方差矩阵。

该方法具体操作为：

step1：水声传感网络节点测量到m对有效的uasn节点间信号传播时间，记为：

tpert＝[tp1,…,tpm,…,tpm]^t；

step2：根据tpert中对应的收发节点关系，计算参考声速剖面c0(z)下对应的信号传播时间，以及对应的声线路径，分别记为tref＝[tr1,…,trm,…,trm]，γref＝[γr1,…,γrm,…,γrm]；

step3：tpert与tref相减得到测量时间扰动，计算公式为δτi＝tpi-tri，再将声线路径γref，以及历史数据得到k阶的eof，代入y＝eα+v得到观测方程；

step4：求解y＝eα+v得到k阶eof系数估计，并代入得到估计声速c(z)；

step5：将待定位节点处的tdoa测量输入常规定位方法，如双曲定位方法，得到初步定位结果x⁽⁰⁾；

step6：参数初始化：对step4求得的声速剖面进行线性分段，求出每段的声速梯度ap以及起始深度处的声速bp；计算测量误差协方差矩阵r∈；设置初始迭代位置x⁽⁰⁾为step5中求得的结果；

step7：根据(s51)～(s57)、(s6)～(s7)求出uasn信标关于迭代位置的h'(·)函数；

step8：求出第k步迭代位置x^(k)对应的模型tdoa观测h'(x^(k))，以及计算x^(k)对应的梯度

step9：将实际tdoa测量δt，模型tdoa观测h'(x^(k))，模型梯度输入下式，得到第k+1步的迭代估计，

step10：重复step8和step9，直到算法收敛；

step11：输出收敛后的位置估计x^(k)。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。