一种基于水下多基地声纳的TDOA偏置削减定位方法

一种基于水下多基地声纳的tdoa偏置削减定位方法
技术领域
1.本发明涉及水下目标源定位技术领域，尤其涉及一种基于水下多基地声纳的tdoa偏置削减定位方法。


背景技术：

2.声纳定位是水下目标源定位的重要手段。传统的声纳按照工作模式不同可以分为被动式声纳和主动式声纳。被动式声纳靠直接接收水下目标机械工作所发出的噪声发现目标，隐蔽性较好。然而随着潜艇隐身技术的深入研究，潜艇发出的噪声越来越小，被动式声纳的定位性能显著下降。主动式声纳通过自主发射声波信号，再接收目标回波对目标进行定位。作用距离相对较远，但是由于需要主动发射信号，隐蔽性不强。相比于上述两种声纳，多基地声纳收发装置是分开放置的，由于发射站可以主动发射信号，具有主动式声纳的优势，并且它的接收站是被动工作的，隐蔽性较好。由于多基地声纳隐蔽性好、抗干扰能力强、机动性能高以及作用距离远的优点，已成为国内外学者的研究热点。
3.多基地声纳定位的原理是由单个或多个发射站发射声波信号，由多个接收站接收目标回波，并根据从信号中获得的时域、频域、空域或者能量域等参数信息对目标进行定位。这些信息包括到达时间、到达时间差(tdoa)、到达频率、到达频率差、到达方位角、到达仰角、接收信号强度和信号到达增益比等。基于上述观测信息，越来越多的定位算法被提了出来。目前的算法一般可以分为闭式解类的算法以及迭代类的算法两种。闭式解类的算法是通过公式推导可以得出求解目标位置的具体公式，它的计算过程更加明确，且计算简洁，计算量较小，因而有很多相关的算法被提了出来；迭代类的算法是通过迭代的过程来求解目标位置，相比于闭式解算法而言，它的计算过程更加复杂，且需要考虑初值选取的问题。
4.闭式解类的算法虽然计算过程简单，且计算的复杂度较低，但一般都需要多次根据观测方程中的各种物理量间的关系进行公式的转换，在此过程中往往会产生很多被忽略的二阶误差项，这样可能造成定位偏置的增大。对于多基地声纳定位来说，其对应的观测方程方程中涉及的物理量更多(有发射站、接收站、tdoa观测量和声速)，进行公式转换时所产生的二阶误差项更多，对系统中的各种误差量更加敏感，因此更需要对其算法的偏置进行削减。针对此问题，本专利将设计一种基于水下多基地声纳的tdoa偏置削减定位方法，以对其定位估计偏置进行削减。


技术实现要素：

5.本发明针对现有的水下多基地声纳tdoa定位方法中偏置较大的问题，提出一种基于水下多基地声纳的tdoa偏置削减定位方法，可以对目前的定位方法的偏置进行削减。
6.为了实现上述目的，本发明采用以下技术方案：
7.本发明首先利用多基地声纳中关于发射站和接收站的到达时间差的观测量，建立相应的观测方程，将所有的观测方程整合为一个矩阵函数的形式并对其进行变形，使其变为一个未知量包含目标位置的伪线性矩阵方程，确定求解该矩阵方程的代价函数，并通过
引入增广矩阵和扩维向量对代价函数进行变形。接着确定一个常数约束条件，并与基于最小二乘准则的求解模型构成新模型，通过求解一个矩阵束的最小广义特征值对应的广义特征向量求解此新模型，得出一个包含目标位置估计的初始解。然后根据一阶泰勒级数展开法确定新模型对于目标源位置估计的误差所满足的方程，并利用最小二乘相关原理得出其解。最后用初始解中对于目标位置的估计减去求出的估计误差的结果，就是最终的目标位置的估计结果。具体如下：
8.本发明的一种基于水下多基地声纳的tdoa偏置削减定位方法，包括：
9.步骤1：利用关于m个发射站和n个接收站的相应的tdoa观测量，建立相应的tdoa观测方程。
10.步骤2：将所有的tdoa观测方程进行整合，并将其变为一个未知量包含目标源位置u的伪线性矩阵方程。
11.步骤3：确定求解步骤2的伪线性方程的代价函数j，并引入增广矩阵a和扩维向量v对代价函数j进行变形。
12.步骤4：基于变形后的代价函数j求解一个常数约束条件，并与基于最小二乘准则的求解模型构成一个未知量包含水下目标源位置u的优化求解模型。
13.步骤5：通过求解矩阵束(a
t
w1a,ω)的最小广义特征值对应的广义特征向量求解优化求解模型，得出一个包含目标位置估计的初始解。
14.步骤6：根据一阶泰勒级数展开法确定优化求解模型对于目标源位置估计的误差所满足的方程，并利用最小二乘相关原理得出其解。
15.步骤7：用步骤5得出的初始解中对于目标位置的估计减去步骤6得出的结果，得出最终的目标位置估计。
16.进一步地，所述步骤1中，第i个发射站的位置为ti，其中1≤i≤m，m表示发射站的个数；第j个接收站的位置为sj，其中1≤j≤n，j表示接收站的个数；u表示一水下目标源；声速为c；那么相关的到达时间差(tdoa)τ
ij
的观测方程为
[0017][0018]
式中中表示声速的先验值，σc表示声速对应的先验误差；
[0019][0019]
表示第i个发射站、第j个接收站对应的带有误差的tdoa观测量，τ
ij
表示第i个发射站、第j个接收站对应的tdoa观测量，
△
τ
ij
表示第i个发射站、第j个接收站对应的tdoa观测量误差；ε
ij
＝c
△
τ
ij
τ
ij
σc，表示对应的误差值。
[0020]
进一步地，所述步骤2中，tdoa观测方程进行变形后的伪线性形式为
[0021][0022]
其中b1表示方程的误差矩阵，ε表示所有ε
ij
构成的向量形式，h1表示观测向量，g1表示观测矩阵，表示方程中的未知量。
[0023]
式中各个元素的具体表达式为
[0024][0025]
ε＝[(ε1)
t (ε2)
t
ꢀ…ꢀ
(εm)
t
]
t
[0026]
[0027][0028][0029]
其中im表示维度为m的单位矩阵；表示矩阵的kronecker积；b
11
表示有关b1的子对角矩阵；εi表示ε的子向量；h
1i
表示h1的子向量；g
1i
表示g1的子矩阵。
[0030]
其中
[0031]b11
＝2diag([||u-s1|| ||u-s2||
ꢀ…ꢀ
||u-sn||])
[0032]
εi＝[ε
i1 ε
i2
ꢀ…ꢀ
ε
in
]
t
[0033][0034][0035]
其中diag(*)表示以向量*中的各项元素为对角的矩阵；οi×j表示一个行数为i，列数为j的零矩阵。
[0036]
进一步地，所述步骤3中，需要求解步骤2的代价函数为
[0037][0038]
式中w1表示相应的加权矩阵，具体表示为
[0039][0040]
e(εε
t
)是有关ε的协方差矩阵，具体表示为
[0041][0042]
其中q
τ
表示tdoa观测量的误差协方差矩阵；τ表示tdoa观测量的矩阵形式，式中
[0043]
τ＝[(τ1)
t (τ2)
t
ꢀ…ꢀ
(τm)
t
]
t
[0044]
τi＝[τ
i1 τ
i2
ꢀ…ꢀ
τ
in
]
t
[0045]
引入增广矩阵a＝[-g
1 h1]和扩维向量则可以将代价函数j变形为
[0046]
j＝v
tat
w1av
[0047]
进一步地，所述步骤4中，需要确定一个常数约束条件，具体推导如下：
[0048]
增广矩阵a可以分解为a＝a
ο

△
a，其中a
ο
表示增广矩阵a的精确已知部分，
△
a表示a的误差部分，分别表示为
[0049][0050]
△
a＝[
‑△
g1ꢀ△
h1]
[0051]
式中表示g1的不含误差项部分，表示h1的不含误差项部分，可表示为
[0052]
[0053][0054]
△
g1表示g1的误差项部分，
△
h1表示h1的误差项部分，可表示为
[0055][0056][0057]
式中
[0058][0059][0060]
△h1i
＝d
i1
△
τi d
i2
σc[0061]
式中d
i1
表示有关
△
τi的误差相关矩阵，d
i2
表示有关σc的误差相关向量，具体表示为
[0062]di1
＝2diag([c||t
i-s1|| c2τ
i1 c||t
i-s2|| c2τ
i2
ꢀ…ꢀ
c||t
i-sn|| c2τ
in
])
[0063][0064]
代价函数j可以表示为
[0065]
j＝v
t
(a
ο

△
a)
t
w1(a
ο

△
a)v
[0066] ＝v
taοt
w1a
ο
v 2v
taοt
w1△
av v
t
△at
w1△
av
[0067]
然后对代价函数的两边取期望可得
[0068]
e(j)＝v
taοt
w1a
ο
v v
t
e(
△at
w1△
a)v
[0069]
将e(j)中的第二项作为一个常数约束条件，则可以构建如下新的关于的优化求解模型
[0070][0071]
其中ω表示二阶误差相关项，k表示任意常数。
[0072]
ω可以具体表示为
[0073][0074]
式中ω1、ω2和ω3分别表示ω的不同子矩阵，可分别表示为
[0075][0076][0077]
其中ω1(m1:m2,:)表示矩阵ω1的第m1～m2行构成的子矩阵；ω1(m1:m2,n1:n2)表示矩阵的第m1～m2行以及第n1～n2列构成的子矩阵；ω2(m:n)表示向量的第m～n列构成的子向量；ω2(m)表示向量ω2的第m个元素；trace(*)表示矩阵*的迹。
[0078]
进一步地，所述步骤5中，若将优化求解模型关于v的估计结果记为那么的估计值可以表示为
[0079][0080]
进一步地，所述步骤6中，根据一阶泰勒级数展开法，中的各个元素可以表示为
[0081][0082]
表示的估计误差，那么根据则可以得出
[0083][0084]
若将优化求解模型对于目标源位置估计的误差记为根据上一组等式，可得如下矩阵方程
[0085][0086]
式中h2表示h1关于的转换形式，g2表示g1关于的转换形式，分别表示为
[0087]
[0088][0089]
根据最小二乘的相关理论，的估计值可以表示为
[0090][0091]
式中w2表示相应的加权矩阵
[0092][0093]
表示的估计协方差矩阵。
[0094]
进一步地，所述步骤7中，目标源位置的最终估计结果可以表示为
[0095][0096]
与现有技术相比，本发明具有的有益效果：
[0097]
本发明利用定位系统中的各种已知量以及先验信息，通过求解最小广义特征值对应的广义特征向量，并且通过最小二乘的基本思想，对水下目标源进行定位，可以有效对水下多基地声纳的tdoa定位方法进行偏置削减。
附图说明
[0098]
图1为本发明实施例一种基于水下多基地声纳的tdoa偏置削减定位方法基本流程图；
[0099]
图2为多基地声纳定位几何示意图；
[0100]
图3为tdoa观测误差为0.001s的定位结果图与误差椭圆曲线；
[0101]
图4为目标位置的估计均方根随tdoa观测误差的变化曲线；
[0102]
图5为目标位置的估计偏置随tdoa观测误差的变化曲线；
[0103]
图6为目标位置为(90,90,90)
t
m的定位结果图与误差椭圆曲线；
[0104]
图7为目标位置的估计均方根随不同目标位置的变化曲线；
[0105]
图8为目标位置的估计偏置随不同目标位置的变化曲线。
具体实施方式
[0106]
下面结合附图和具体的实施例对本发明做进一步的解释说明：
[0107]
如图1所示，本发明的一种基于水下多基地声纳的tdoa偏置削减定位方法的具体实施步骤如下：
[0108]
步骤1：利用关于m个发射站和n个接收站的相应的tdoa观测量，建立相应的tdoa观测方程。
[0109]
步骤2：将所有的tdoa观测方程进行整合，并将其变为一个未知量包含目标源位置u的伪线性矩阵方程。
[0110]
步骤3：确定求解步骤2的伪线性方程的代价函数j，并引入增广矩阵a和扩维向量v对代价函数j进行变形。
[0111]
步骤4：基于变形后的代价函数j求解一个常数约束条件，并与基于最小二乘准则的求解模型构成一个未知量包含水下目标源位置u的优化求解模型。
[0112]
步骤5：通过求解矩阵束(a
t
w1a,ω)的最小广义特征值对应的广义特征向量求解优化求解模型，得出一个包含目标位置估计的初始解。
[0113]
步骤6：根据一阶泰勒级数展开法确定优化求解模型对于目标源位置估计的误差所满足的方程，并利用最小二乘相关原理得出其解。
[0114]
步骤7：用步骤5得出的初始解中对于目标位置的估计减去步骤6得出的结果，得出最终的目标位置估计。
[0115]
进一步地，所述步骤1中，第i个发射站的位置为ti，其中1≤i≤m，m表示发射站的个数；第j个接收站的位置为sj，其中1≤j≤n，j表示接收站的个数；u表示一水下目标源；声速为c；那么相关的到达时间差(tdoa)τ
ij
的观测方程为
[0116][0117]
式中式中表示声速的先验值，σc表示声速对应的先验误差；
[0118][0118]
表示第i个发射站、第j个接收站对应的带有误差的tdoa观测量，τ
ij
表示第i个发射站、第j个接收站对应的tdoa观测量，
△
τ
ij
表示第i个发射站、第j个接收站对应的tdoa观测量误差；ε
ij
＝c
△
τ
ij
τ
ij
σc，表示对应的误差值。
[0119]
多基地声纳定位几何示意图如图2所示。
[0120]
进一步地，所述步骤2中，tdoa观测方程进行变形后的伪线性形式为
[0121][0122]
其中b1表示方程的误差矩阵，ε表示所有ε
ij
构成的向量形式，h1表示观测向量，g1表示观测矩阵，表示方程中的未知量。
[0123]
式中各个元素的具体表达式为
[0124][0125]
ε＝[(ε1)
t (ε2)
t
ꢀ…ꢀ
(εm)
t
]
t
[0126][0127][0128][0129]
其中im表示维度为m的单位矩阵；表示矩阵的kronecker积；b
11
表示有关b1的子对角矩阵；εi表示ε的子向量；h
1i
表示h1的子向量；g
1i
表示g1的子矩阵。
[0130]
其中
[0131]b11
＝2diag([||u-s1|| ||u-s2||
ꢀ…ꢀ
||u-sn||])
[0132]
εi＝[ε
i1 ε
i2
ꢀ…ꢀ
ε
in
]
t
[0133][0134][0135]
其中diag(*)表示以向量*中的各项元素为对角的矩阵；οi×j表示一个行数为i，列数为j的零矩阵。
[0136]
进一步地，所述步骤3中，需要求解步骤2的代价函数为
[0137][0138]
式中w1表示相应的加权矩阵，具体表示为
[0139][0140]
e(εε
t
)是有关ε的协方差矩阵，具体表示为
[0141][0142]
其中q
τ
表示tdoa观测量的误差协方差矩阵；τ表示tdoa观测量的矩阵形式，式中
[0143]
τ＝[(τ1)
t (τ2)
t
ꢀ…ꢀ
(τm)
t
]
t
[0144]
τi＝[τ
i1 τ
i2
ꢀ…ꢀ
τ
in
]
t
[0145]
引入增广矩阵a＝[-g
1 h1]和扩维向量则可以将代价函数j变形为
[0146]
j＝v
tat
w1av
[0147]
进一步地，所述步骤4中，需要确定一个常数约束条件，具体推导如下：
[0148]
增广矩阵a可以分解为a＝a
ο

△
a，其中a
ο
表示增广矩阵a的精确已知部分，
△
a表示a的误差部分，分别表示为
[0149][0150]
△
a＝[
‑△
g1ꢀ△
h1]
[0151]
式中表示g1的不含误差项部分，表示h1的不含误差项部分，可表示为
[0152][0153][0154]
△
g1表示g1的误差项部分，
△
h1表示h1的误差项部分，可表示为
[0155][0156]
[0157]
式中
[0158][0159][0160]
△h1i
＝d
i1
△
τi d
i2
σc[0161]
式中d
i1
表示有关
△
τi的误差相关矩阵，d
i2
表示有关σc的误差相关向量，具体表示为
[0162]di1
＝2diag([c||t
i-s1|| c2τ
i1 c||t
i-s2|| c2τ
i2
ꢀ…ꢀ
c||t
i-sn|| c2τ
in
])
[0163][0164]
代价函数j可以表示为
[0165]
j＝v
t
(a
ο

△
a)
t
w1(a
ο

△
a)v
[0166] ＝v
taοt
w1a
ο
v 2v
taοt
w1△
av v
t
△at
w1△
av
[0167]
然后对代价函数的两边取期望可得
[0168]
e(j)＝v
taοt
w1a
ο
v v
t
e(
△at
w1△
a)v
[0169]
将e(j)中的第二项作为一个常数约束条件，则可以构建如下新的关于的优化求解模型
[0170][0171]
其中ω表示二阶误差相关项，k表示任意常数。
[0172]
ω可以具体表示为
[0173][0174]
式中ω1、ω2和ω3分别表示ω的不同子矩阵，可分别表示为
[0175][0176]
[0177][0178]
其中ω1(m1:m2,:)表示矩阵ω1的第m1～m2行构成的子矩阵；ω1(m1:m2,n1:n2)表示矩阵的第m1～m2行以及第n1～n2列构成的子矩阵；ω2(m:n)表示向量的第m～n列构成的子向量；ω2(m)表示向量ω2的第m个元素；trace(*)表示矩阵*的迹。
[0179]
进一步地，所述步骤5中，若将优化求解模型关于v的估计结果记为那么的估计值可以表示为
[0180][0181]
进一步地，所述步骤6中，根据一阶泰勒级数展开法，中的各个元素可以表示为
[0182][0183]
表示的估计误差，那么根据则可以得出
[0184][0185]
若将优化求解模型对于目标源位置估计的误差记为根据上一组等式，可得如下矩阵方程
[0186][0187]
式中h2表示h1关于的转换形式，g2表示g1关于的转换形式，分别表示为
[0188][0189][0190]
根据最小二乘的相关理论，的估计值可以表示为
[0191]
[0192]
式中w2表示相应的加权矩阵
[0193][0194]
表示的估计协方差矩阵。
[0195]
进一步地，所述步骤7中，目标源位置的最终估计结果可以表示为
[0196][0197]
为验证本发明效果，进行如下具体示例：
[0198]
假设多基地声纳中有2个发射站和5个接收站，2个发射站的位置分比为(600,900,600)
t
m，(600,700,800)
t
m；5个接收站的位置分比为(-500,500,600)
t
m，(600,-600,600)
t
m，(700,700,-600)
t
m，(-600,-600,600)
t
m，(750,-600,-700)
t
m。声速在水下的值为1500m/s。
[0199]
(1)假设声速对应的先验误差为0.5m/s，tdoa观测误差假设为一个随自变量σ1变化的值0.001σ1s，图3给出了本专利公开的偏置削减定位方法当tdoa观测误差为0.001s的定位结果图与误差椭圆曲线，可以直观看出本专利公开的方法的定位结果的基本情况；图4给出了目标位置的估计均方根随tdoa观测误差的变化曲线，可以看出本专利公开的偏置削减定位方法与现有的多步加权算法以及泰勒级数法相比，估计均方误差不会改变，且不存在泰勒级数法初值选取不合适估计均方误差显著增大的情况；图5给出了目标位置的估计偏置随tdoa观测误差的变化曲线，可以看出本专利公开的偏置削减定位方法与现有的多步加权算法相比，偏置明显进行了削减，且与初值为真实值的泰勒级数法的偏置接近，同时没有初值选取不合适时估计偏置增大的问题。
[0200]
(2)假设声速对应的先验误差为0.5m/s，目标位置可以表示为(50 40σ2,50 40σ2,50 40σ2)
t
m，图6给出了本专利公开的偏置削减定位方法当目标位置为(50,50,50)
t
m时的定位结果图与误差椭圆曲线，可以直观看出本专利公开方法的定位结果的基本情况；图7给出了目标位置的估计均方根随不同目标位置的变化曲线，可以看出本专利公开的偏置削减定位方法与现有的多步加权算法以及泰勒级数法相比，估计均方误差不会改变，且不存在泰勒级数法初值选取不合适估计均方误差显著增大的情况；图8给出了目标位置的估计偏置随不同目标位置的变化曲线，可以看出本专利公开的偏置削减定位方法与现有的多步加权算法相比，偏置明显进行了削减，且与初值为真实值的泰勒级数法的偏置接近，同时没有初值选取不合适时估计偏置增大的问题。
[0201]
综上，本发明利用定位系统中的各种已知量以及先验信息，通过求解最小广义特征值对应的广义特征向量，并且通过最小二乘的基本思想，对水下目标源进行定位，可以有效对水下多基地声纳的tdoa定位方法进行偏置削减。
[0202]
以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。