一种基于TDoA的搜索加权定位方法

一种基于tdoa的搜索加权定位方法
技术领域
1.本发明涉及定位与导航技术领域，尤其涉及一种基于tdoa的搜索加权定位方法。


背景技术：

2.信号到达时间差(time difference of arrival，tdoa)是定位领域被广泛使用的测量之一，它是指目标接收到不同基站发射信号的时间差，或者是不同的基站接收到目标发射信号的时间差。基于tdoa的定位系统的优点是不需要基站和目标间的时间同步，从而降低了系统设计的成本和复杂度。
3.基于tdoa的定位算法已得到广泛的研究，主要分为3类：非线性优化方法、线性方法和几何方法。非线性优化方法有高斯-牛顿法、梯度下降法等，缺点是性能与选择的初值有关，在迭代求解过程中可能陷入局部最优解；线性方法主要有两步加权最小二乘算法、约束加权最小二乘算法和分离的约束加权最小二乘算法等，两步加权最小二乘算法和约束加权最小二乘算法在噪声较小的情况下达到克拉美罗下界，但是在特殊的几何条件下时定位效果很差，分离的约束加权最小二乘算法通过把位置坐标与辅助变量分离在等式的两边，可避免系统矩阵的病态。总之，由于加权最小二乘类方法忽视了噪声的高次项，它在较小的量测噪声下可以接近获得克拉美罗下届的估计性能，但是随着量测噪声的增加，估计的误差会偏离克拉美罗下届，尤其在基站呈均匀圆形阵列分布且目标靠近基站时，抗噪能力较差；几何方法主要包括缩圆法，这类方法计算计算复杂度低，但无法利用冗余的量测信息来提高定位精度。


技术实现要素：

4.本发明的目的是提供一种基于tdoa的搜索加权定位方法，该方法不仅考虑了测量误差的影响，还考虑了几何条件对位置估计的影响，充分利用了冗余信息，提高了位置估计精度。
5.本发明的目的是通过以下技术方案实现的：
6.一种基于tdoa的搜索加权定位方法，所述方法包括：
7.步骤1、当基站数目等于3时，通过建立方程直接求得一个解析解，该解析解即为最终目标估计位置；
8.步骤2、当基站数目大于3时，计算三基站的组合个数，并针对每个组合基于tdoa求解出目标位置的解析解；
9.步骤3、基于步骤2求出的目标位置的解析解，构造矩形搜索区域；
10.步骤4、然后在所构造的矩形搜索区域内选取候选位置，针对所有的候选位置计算对应候选位置的权重；
11.步骤5、将步骤4计算的权重代入构建的目标函数中，计算得到所有候选位置对应的目标函数值；其中，所述目标函数是测量值与真实值之间偏差的加权平方和；
12.步骤6、从计算得到的所有候选位置对应的目标函数值中找到最小值，该最小值对
应的候选位置就是最终目标估计位置。
13.由上述本发明提供的技术方案可以看出，上述方法不仅考虑了测量误差的影响，还考虑了几何条件对位置估计的影响，充分利用了冗余信息，提高了位置估计精度，能反映众多测量中是否存在异常，并具有屏蔽异常测量的能力。
附图说明
14.为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
15.图1为本发明实施例提供的基于tdoa的搜索加权定位方法流程示意图；
16.图2为本发明实施例中矩形搜索区域的示意图；
17.图3为本发明实施例中坐标系变换后候选位置与单支双曲线的最短距离示意图；
18.图4为本发明实施例一般情况下搜索方法求解r示意图；
19.图5为本发明实施例特殊情况下的r求解示意图。
具体实施方式
20.下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，这并不构成对本发明的限制。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
21.如图1所示为本发明实施例提供的基于tdoa的搜索加权定位方法流程示意图，所述方法包括：
22.步骤1、当基站数目等于3时，通过建立方程直接求得一个解析解，该解析解即为最终目标估计位置；
23.在该步骤中，以二维为例，同样适用于三维情形，具体过程为：
24.假设第i个基站的坐标为si(xi，yi，hi)，其中hi是基站i与目标间的垂直高度差，i＝1，2，
…
，n，n是基站总数，第k时刻目标所处的实际位置为pk(xk，yk)，对其的估计为
25.实现基于tdoa的定位至少需要3个基站，即要求n≥3，已知的tdoa量测可转化为与距离相关的量测，如式(1)所示：
[0026][0027]
其中，c为声速，单位为m/s，是以基站1作为参考点，目标与基站si间以及与参考基站间的距离之差；
[0028]
当n＝3时，用量测代替真实的rdoa
i1，k
，求得的解是对pk的估计，具体求解过程如下：
[0029][0030]
式(2)中的
②
和
③
分别减去
①
得：
[0031][0032]
其中为目标与基站1之间的距离；x
12
＝x
1-x2，y
12
＝y
1-y2，x
13
＝x
1-x3，y
13
＝y
1-y3，，
[0033]
式(3)的求解分为两种情况：
[0034]
三基站不在同一条直线，即x
12y13
≠x
13y12
；三基站在同一条直线上，即x
12y13
＝x
13y12
；在这两种情况下，三个基站均在不同的位置上；
[0035]
当x
12y13
≠x
13y12
时，从式(3)可得到与的关系，见式(4)：
[0036][0037]
其中，
[0038][0039]
把式(4)代入式(2)中的
①
式可得到一个关于的一元二次方程，如式(6)：
[0040][0041]
其中，
[0042]
[0043]
显然式(6)有两个解，但只能选择其中一个作为的解，令δ＝b
2-4ac，选择的准则如下：
[0044]
①
若δ≥0，式(6)有两个实数根，即(-b+√δ)/(2a)和(-b-√δ)/(2a)，如果两个根都是正值，则选择较小的作为的解，否则选择较大的；
[0045]
②
若δ＜0，式(6)有两个复数根，那么取其实部作为的解，即
[0046]
把代入式(4)中可求得目标的估计位置
[0047]
当x
12y13
＝x
13y12
时，此时式(4)是不成立的，因为式(5)中的分母均为零，因此需要通过如下的方法求解式(3)，分为3种情形：
[0048]
1)若x
12
＝0，一定有y
12
≠，否则基站1和基站2在同一位置，从而有x
13
＝0且y
13
≠0，令式(3)中的
①
式两边乘以y
13
，且
②
式两边乘以y
12
，那么可通过式(8)求解：
[0049][0050]
2)若y
12
＝0，同情形
①
，一定有x
12
≠0，那么有x
13
≠0且y
13
＝0，此时通过式(9)求解：
[0051][0052]
3)若x
12
≠0且y
12
≠0，那么可通过式(8)或者(9)求解；
[0053]
在求出后，把其代入式(2)中的式
①
和
②
可求得目标的估计位置
[0054]
具体实现中，当n＝3时的目标位置解算方法，不仅可以求解三基站不在一条直线上的情形，对于三基站在同一条直线上的情形仍可求解。
[0055]
步骤2、当基站数目大于3时，计算三基站的组合个数，并针对每个组合基于tdoa求解出目标位置的解析解；
[0056]
在该步骤中，当基站数目大于3，即n》3时，对n个基站，共有m个3基站组合，m值由下式确定，其中
‘
c’表示组合操作符，
‘
！’表示求阶乘操作符：
[0057][0058]
针对m个3基站组合，每个基站组合利用步骤1中n＝3时的求解方法，求出目标估计位置的解析解，共m个，形成一个集合s。
[0059]
步骤3、基于步骤2求出的目标位置的解析解，构造矩形搜索区域；
[0060]
在该步骤中，如图2所示为本发明实施例中矩形搜索区域的示意图，构造矩形搜索区域的过程具体为：
[0061]
利用所有三基站组合解的横坐标最小值和最大值以及纵坐标最小值和最大值构建，矩形搜索区域的顶点坐标分别为：
[0062]
和
[0063]
步骤4、然后在所构造的矩形搜索区域内选取候选位置，针对所有的候选位置计算对应候选位置的权重；
[0064]
在该步骤中，在所构造的矩形搜索区域内选取候选位置的过程具体为：
[0065]
如图2所示，使用m条水平线和n条竖直线等分矩形搜索区域，以水平线与竖直线的交点、水平线与矩形边的交点、竖直线与矩形边的交点、矩形顶点构成所有候选位置，共计(m+2)
·
(n+2)个。值得注意的是，m和n的取值根据需要设置，可以相等，也可以根据矩形边长的比例选取。
[0066]
具体实现中，对应候选位置的权重计算过程为：
[0067]
首先测量噪声在侯选位置与两基站间的距离差上得到
[0068]
经过坐标系的旋转和平移，针对一般情况，以两基站坐标和构建单支标准双曲线，求解所述候选位置到该单支标准双曲线的最短距离r；
[0069]
针对特殊情况，以两基站坐标和构建直线或射线，求解所述候选位置到该直线或射线的最短距离r；
[0070]
则权重α
i，j
按照如下公式计算：
[0071]
α
i，j
＝1/r2。
[0072]
举例来说，所述权重α
i，j
的具体计算过程为：
[0073]
选择两个不同的基站si和sj，其坐标分别为si(xi，yi)，sj(xj，yj)；
[0074]
选取一个候选位置，坐标为t(x
t
，y
t
)，那么目标未知待定位置与基站si和sj的距离差rdoa
i，j
的值可由式(12)计算得到：
[0075][0076]
设置当前测量噪声为δd
i，j
，令
[0077]
假设有一个位置点q，满足其中||qsi||表示点q与基站si的距离；
[0078]
在一般情况下，即当时，点q在一条双曲线上；
[0079]
在特殊情况下，当时，点q在一条直线上，当时，点q在射线上；
[0080]
①
针对一般情况，即
[0081]
为计算权重的方便，首先需要对坐标系进行变换，即进行旋转和平移操作，使得在
坐标系变换后，基站si和sj处在x轴上，且它们的中心处在原点上，坐标系旋转、平移的方法见式(13)：
[0082][0083]
其中：
[0084][0085][0086]
坐标系变换后，以为距离差构建单支标准双曲线，经坐标系的旋转和平移，构建的单支标准双曲线方程如式(14)：
[0087][0088]
其中：
[0089]
b2＝c
2-a2[0090]
若则取单支标准双曲线(14)的左支；若则取单支标准双曲线(14)的右支；
[0091]
假设候选位置t(x
t
，y
t
)经公式(13)旋转、平移之后的坐标为t
′
(x
′
t
，y
′
t
)，以t
′
(x
′
t
，y
′
t
)为圆心画圆，当该圆与单支标准双曲线(14)的单支相切时，半径r就是待定点位置至双曲线的最短距离；r越大意味着当前测量误差和几何条件下，对应测量对最终位置估计的精度贡献越小，相应的权重应越小；反之，r越小意味着对应测量对最终位置估计的精度贡献越大，权重越大。
[0092]
假设切点坐标为p
′0(x
′0，y
′0)，那么有r＝||t
′
p
′0||，如图3所示为本发明实施例中坐标系变换后候选位置与单支双曲线的最短距离示意图，然而p
′0(x
′0，y
′0)不易直接求解，故本发明为此提出了一种数值搜索求解方法，即在单支标准双曲线上搜索到一个点p
′0(x
′0，y
′0)，满足[(x
′
t-x
′0)2+(y
′
t-y
′0)2]最小；
[0093]
为加快数值求解速度，如图4所示为本发明实施例所述一般情况下搜索方法求解r示意图，本发明限定搜索范围在单支标准双曲线上的a和b两点之间，其中点a(x
′a，y
′a)是过点t
′
(x
′
t
，y
′
t
)，且与渐近线垂直的直线与单支标准双曲线的交点；点b(x
′b，y
′b)是直线y
′
＝y
′
t
与单支标准双曲线的交点；
[0094]
下面以右支标准双曲线为例，介绍点a和b坐标的确定过程：
[0095]
当y
′
t
≥0时，过点t
′
(x
′
t
，y
′
t
)且与渐近线垂直的直线方程取：
[0096][0097]
而当y
′
t
＜0时，过点t
′
(x
′
t
，y
′
t
)且与渐近线垂直的直线方程取：
[0098][0099]
令将式(15)代入式(14)有：
[0100][0101]
求解式(17)可得到两个解，由于双曲线取右支，故取其中一个正值作为x
′a；把x
′a的值代入式(15)求得y
′a；同理，把式(16)代入式(14)中，可求得当y
′
t
＜0时的a和b的坐标；
[0102]
把y
′
＝y
′
t
代入式(14)得：
[0103][0104]
求解式(18)可得到两个解，由于双曲线取右支，故取其中一个正值作为x
′b，y
′b＝y
′
t
；
[0105]
在确定搜索范围之后，使用折半查找方法快速缩小查找范围，直至||ab||小于阈值ε1；最后，r＝min(||t
′
a||，||t
′
b||)；
[0106]
仍以右支标准双曲线为例，折半查找的过程如下：
[0107]
1)将y
′
＝(y
′a+y
′b)/2代入式(14)得：
[0108][0109]
选取正值的根作为x
′c的值，y
′c＝(y
′1+y
′2)/2，即直线y
′
＝(y
′1+y
′2)/2与右支标准双曲线的交点c坐标为(x
′c，y
′c)；
[0110]
2)比较||t
′
a||、||t
′
b||和||t
′
c||；若||t
′
a||最大，搜索范围缩小至bc之间，更新c为a；若||t
′
b||最大，搜索范围缩小至ac之间，令c为b；
[0111]
3)计算||ab||，假如||ab||小于阈值ε1，则终止循环，r＝min(||t
′
a||，||t
′
b||)，否则继续；
[0112]
4)重复步骤(1)～(3)；
[0113]
同理，采用同样的方法可确定候选位置与左支标准双曲线的最短距离r；
[0114]
②
针对特殊情况，即或或
[0115]
当双曲线变成了直线方程：x
′
＝0，如图5所示为本发明实施例在特殊情况下的r求解示意图，此时r＝x
′
t
；
[0116]
当时，即双曲线变为以s
′j为端点的射线，若x
′
t
≥x
′j，则r＝y
′
t
；否则，r＝||t
′
sj||；
[0117]
当时，即双曲线变为以s
′i为端点的射线，若x
′
t
≤x
′i，则r＝y
′
t
；否则，r＝||t
′
si||；
[0118]
最后，权重α
i，j
按照如下公式(20)计算：
[0119]
α
i，j
＝1/r2(20)。
[0120]
步骤5、将步骤4计算的权重代入构建的目标函数中，计算得到所有候选位置对应的目标函数值；
[0121]
其中，所述目标函数是测量值与真实值之间偏差的加权平方和，所述目标函数f(x，y)表示为：
[0122][0123]
其中,c
n2
表示从n个基站中任意选取2个，是组合公式的表示形式；rdoa
i，j
为目标未知待定位置与基站si和sj的距离差；是目标距离两个基站距离差的测量值；α
i，j
为权重；
[0124]
由式(11)即可得到所有候选位置对应的目标函数值。
[0125]
步骤6、从计算得到的所有候选位置对应的目标函数值中找到最小值，该最小值对应的候选位置就是最终目标估计位置。
[0126]
在该步骤中，具体是通过找到一个位置使得目标函数(11)的值最小，即误差的加权平方和最小，那么该最小值对应的候选位置就是最终目标估计位置，记为pf。
[0127]
值得注意的是，本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。
[0128]
综上所述，本发明实施例所述方法具有如下优点：
[0129]
1、与常见的加权最小二乘(weightedleastsquare,wls)类算法相比，本方法不存在忽略噪声的高次项问题，所以更加鲁棒，且针对基站在同一条直线上的情形，仍可以求解；
[0130]
2、与传统的非线性优化方法相比，本方法避免了初值的选择，且防止陷入局部最优解，更重要的是，本方法不仅考虑了测量误差的影响，还考虑了几何条件对位置估计的影响，而传统的方法往往忽视了几何条件这一重要因素；
[0131]
3、与几何方法相比，本方法充分利用了冗余信息，提高了位置估计精度；
[0132]
4、在测量冗余的条件下，本方法能反映众多测量中是否存在异常的测量值，即误差极大的观测，并具有屏蔽异常测量的能力，实现更加鲁棒的位置估计，这是其他任何基于tdoa的定位方法无法比拟的。
[0133]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，
都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。本文背景技术部分公开的信息仅仅旨在加深对本发明的总体背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。