本发明涉及一种tdoa定位方法,具体的说,涉及了一种接收机位置误差下多无人机目标的tdoa定位方法。
背景技术:
近几十年来,利用无人机上的辐射源信号到达不同接收机的时差(timedifferenceofarrival,tdoa)来对无人机定位的问题得到了广泛的关注。然而,由于tdoa测量关于辐射源位置参数高度非线性,利用tdoa测量对无人机目标进行定位并非易事。
得益于国内外对基于tdoa的辐射源定位问题的持续研究,目前基于tdoa的辐射源定位问题已经发展出了一些有效算法,其中,两步加权最小二乘算法作为处理tdoa定位问题的经典线性方法,一直备受关注。y.t.chan等人针对基于tdoa的辐射源定位问题,提出了经典的两步加权最小二乘优化思想,其基本思想是:在第一步加权最小二乘估计中,将辐射源目标到参考接收机的距离作为辅助参数,从而实现tdoa测量方程的伪线性化;然后在第二步加权最小二乘估计中,利用辅助参数与辐射源位置参数之间的函数关系来构造方程,从而进一步提高辐射源位置的估计精度。理论分析和仿真实验表明,在tdoa测量误差较小时,2wls算法的定位误差可以达到克拉美罗界(cramér-raolowerbound,crlb)。chan的算法建立在接收机位置准确已知的基础上,然而现代定位系统中,接收机通常被布设在卫星、飞机、空中无人机等运动平台上在动态移动平台上,这将导致接收机的位置误差增加。当接收机位置误差增大,从而无法忽略时,将导致y.t.chan的算法显著增加,从而无法达到crlb。k.c.ho等人定量分析了接收机位置误差对辐射源定位精度的影响,并针对接收机位置存在误差时,提出了适用于接收机位置误差下的2wls代数解算法。理论推导和仿真结果表明,ho的算法在接收机位置出现误差时,定位精度依然可以达到crlb。但是,上述算法均假设场景中仅存在单个辐射源目标。而实际应用中,定位场景中常常会同时出现多个辐射源目标(例如典型的无人机探测场景)。当待监视区域中同时出现多个辐射源目标时,在利用数据关联技术建立起辐射源目标和tdoa测量的关联配对后,可以利用上述针对单辐射源场景设计的定位算法分别对各辐射源进行定位。然而,由于对应于这多个辐射源目标的tdoa测量误差之间的统计相关性,以及多个辐射源目标的tdoa测量方程与接收机位置误差的函数关系,这种利用针对单辐射源场景设计的定位算法分别对各辐射源进行定位的方式,得到的辐射源位置估计理论上并非全局最优。l.yang等人在文献中提出了一种适用于多辐射源场景的tdoa定位代数解算法,该算法将经典的两步加权最小二乘思想扩展到多辐射源场景,其在第一步加权最小二乘估计中将各辐射源到参考接收机距离作为辅助参数,从而将非线性的rd测量方程转化为伪线性形式,并利用加权最小二乘估计得到各辐射源位置和辅助参数的估计,然后在第二步加权最小二乘估计中利用辅助参数与各辐射源位置参数之间的函数关系来构建一组线性方程,从而提高辐射源位置的估计精度。理论分析和仿真实验表明,yang的算法在多辐射源场景下定位精度可以达到crlb,但是,yang的算法方法具有较高的计算复杂度,因为它需要同时估计辐射源和接收机的位置。为此,m.sun等人在文献针对多辐射源目标场景,进一步改进了yang的算法,提出了一种更高效的2wls代数解算法。sun的算法无需估计接收机位置,因此计算复杂度更低。理论分析和仿真实验表明,sun的算法可以达到crlb。然而,经过分析可知,sun的算法的第二步加权最小二乘估计结果需要进行开方处理,并且开方后解的符号需要借助第一步加权最小二乘估计得到的并不准确的估计结果来确定。此外,在一些特殊场景下,例如当目标某一坐标与参考接收机的对应坐标重合和接近时,sun的算法的第二步加权最小二乘估计还会出现缺秩问题,导致算法的定位稳健性和定位精度下降。因此,有必要针对接收机位置误差下的多辐射源场景下的tdoa定位问题进行进一步研究。
为了解决以上存在的问题,人们一直在寻求一种理想的技术解决方案。
技术实现要素:
本发明的目的是针对现有技术的不足,从而提供了一种接收机位置误差下多无人机目标的tdoa定位方法。
为了实现上述目的,本发明所采用的技术方案是:一种接收机位置误差下多无人机目标的tdoa定位方法,每个无人机上预置1个辐射源,其包括以下步骤:
s1,利用m个接收机截获n个辐射源的信号,提取对应的tdoa测量信息,构建rd测量方程:
r=r°+δr(1)
其中,r为rd测量向量矩阵,且r=[r1t,r2t,…,rnt]t;ri为辐射源i到接收机j=2,3,…,m的测量距离与其到达参考接收机1的测量距离之差的向量矩阵,且ri=[r21,i,r31,i,…,rm1,i]t,i=(1,2,…n);rj1,i为辐射源i到接收机j=2,3,…,m的测量距离与其到达参考接收机的测量距离之差,且
s2,通过引入辅助变量以及对rd测量方程进行非线性变换,构建线性估计方程;
s201,引入辅助变量
根据式(2)计算获得辅助变量θ的估计协方差矩阵:
式中,w1为加权矩阵,
且
δh1表示误差向量矩阵,δh1,i=b1,iδri+d1,iδs(7)
其中,
且
s202,根据辅助变量
其中,w2为加权矩阵,且
s3,采用迭代加权最小二乘法获得辐射源目标的定位结果估计值;
步骤s301,令迭代次数k=0,设置加权矩阵
步骤s302,将加权矩阵w1代入式(2)中进行计算,获得辅助变量θ的初步估计结果
步骤s303,将辅助变量θ的初步估计结果
步骤s304,将b1,i,d1,i代入式(8),求得b1,d1;
步骤s305,将b1,d1代入式(4),计算更新w1,累计迭代次数k,判断k是否大于2,若不大于,则返回执行步骤s302;若大于,则执行步骤s306;
步骤s306,将更新后的加权矩阵w1代入式(3)中,得到辅助变量θ的初步估计结果
步骤s307,将协方差矩阵
步骤s308,将加权矩阵w2代入式(10),即得到最终定位结果。
基于上述,步骤s1的具体为:
s101,以接收机1作为参考接收机,计算辐射源i发射的信号到达接收机j的时间与其到达参考接收机的时间之差
其中,
s102,对公式(13)做变形,得到辐射源i和接收机j的rd测量公式:
考虑到实际中测量误差的影响,则进一步将辐射源i和接收机j的rd测量公式表示为:
对于系统中的n个辐射源和m个接收机,一共可以得到n(m-1)个rd测量公式;
s103,定义
则上述n(m-1)个rd测量公式表示为:
r=r°+δr(1)
其中,rd测量误差向量矩阵δr服从零均值的高斯分布,其协方差矩阵为qr=e(δrδrt)。
基于上述,s201的具体步骤如下:
将辅助参数与辐射源目标位置参数之间的约束关系为
对于辐射源i,引入辅助向量
其中,
δh1,i表示为误差向量矩阵,且δh1,i=b1,iδri+d1,iδs;
b1,i为(m-1)×(m-1)的矩阵,d1,i为(m-1)×3m的矩阵,二者的具体表达式为
对于n个辐射源目标,令
g1θ°=h1+δh1(16)
式中,
δh1为总的误差向量,且δh1=b1δr+d1δs;
对式(16)进行求解,获得加权最小二乘解:
其中,w1为加权矩阵,且w1选取为:
忽略矩阵g1和d1中的误差,计算获得式(2)中θ的估计误差δθ为:
将式(17)乘以其转置,并求期望,可得θ的估计协方差矩阵为:
基于上述,步骤s202的具体步骤如下:
将辅助参数与辐射源目标位置参数之间的约束关系为
令
将
对于辐射源i,构建如下方程:
其中,
δh2,i为误差向量,且δh2,i=b2,iδθi+d2,iδs1,δθi为辅助变量θi的估计误差;
对于n个辐射源,获得n个矩阵
g2u°=h2+δh2(19)
式中,
δh2为总的误差向量,且δh2=b2δθ+d2δs,δθ为辅助变量矩阵θ的估计误差;
其中,
对式(19)进行求解,获得加权最小二乘解:
其中,w2为加权矩阵,且w2选取为:
本发明相对现有技术具有突出的实质性特点和显著的进步,具体的说,本发明按照两步加权最小二乘思想框架,推导辐射源目标位置的代数解。首先在第一步加权最小二乘估计中,通过对rd测量方程进行非线性变换和引入辅助参数的方式实现对rd测量方程的伪线性化处理,并粗略估计辐射源目标位置;接着,在第二步加权最小二乘估计中,利用辐射源位置与辅助参数之间的函数关系提出了一种新的线性化方法,由于式(19)中待估参量为辐射源位置u°,因此后续对其进行求解时,不再需要进行额外的开方操作,也就避免了开方后的符号模糊问题;同时由于式(12)中,并不包含辐射源位置与参考接收机位置之差,因此,当任意辐射源的任意一维坐标与参考接收机的对应坐标相等时,也不会出现矩阵缺秩问题,能够精确估计辐射源目标位置。
附图说明
图1是本发明的接收机位置误差下多辐射源目标定位场景。
图2是仿真定位场景1。
图3是不同rd测量误差下算法的定位误差。
图4是不同接收机位置误差下算法的定位误差。
图5是仿真定位场景2。
图6是不同角度下算法的稳健性比较。
具体实施方式
下面通过具体实施方式,对本发明的技术方案做进一步的详细描述。
如图1所示,本文所针对的三维定位场景中存在n个无人机,每个无人机上设置有一个辐射源,故每个无人机作可以为一个辐射源目标,;将待估计的全部辐射源位置参数表示为一个3n×1的列向量
实际中接收机位置存在误差的非理想定位场景,因此,这里的接收机真实位置向量s°是未知的,无法直接用于辐射源定位。取而代之地,我们只能得到含有误差的接收机位置测量值
那么,数学上可以表示为
其中,δsj=[δxs,j,δys,j,δzs,j]t,j=1,2,...,m为接收机j位置测量误差向量。
对于系统中的m个接收机,其真实位置可以合并表示为一个3m×1的列向量s=s°+δs,其中,
为了估计n个辐射源目标的位置向量u°,本发明提出一种接收机位置误差下多无人机目标的tdoa定位方法,其包括以下步骤:
s1,利用m个接收机截获n个辐射源的信号,提取对应的tdoa测量信息,构建rd测量方程:
r=r°+δr(1)
其中,r为rd测量向量矩阵,且r=[r1t,r2t,…,rnt]t;ri为辐射源i到接收机j=2,3,…,m的测量距离与其到达参考接收机1的测量距离之差的向量矩阵,且ri=[r21,i,r31,i,…,rm1,i]t,i=(1,2,…n);rj1,i为辐射源i到接收机j=2,3,…,m的测量距离与其到达参考接收机的测量距离之差,且
s2,通过引入辅助变量以及对rd测量方程进行非线性变换,构建线性估计方程;
s201,引入辅助变量
根据式(2)计算获得辅助变量θ的估计协方差矩阵
其中,w1为加权矩阵,
且g1=diag([g1,1,g1,2,...g1,n]),
δh1表示误差向量矩阵,δh1,i=b1,iδr+d1,iδs(7),
其中,
s202,根据辅助变量
其中,w2为加权矩阵,且
s3,采用迭代加权最小二乘法获得辐射源目标的定位结果估计值;
步骤s301,令迭代次数k=0,设置加权矩阵
步骤s302,将加权矩阵w1代入式(2)中进行计算,获得辅助变量θ的初步估计结果
步骤s303,将辅助变量θ的初步估计结果
步骤s304,将b1,i,d1,i代入式(8),求得b1,d1;
步骤s305,将b1,d1代入式(4),计算更新w1,累计迭代次数k,判断k是否大于2,若不大于,则返回执行步骤s302;若大于,则执行步骤s306;
步骤s306,将更新后的加权矩阵w1代入式(3)中,得到辅助变量θ的初步估计结果
步骤s307,将协方差矩阵
步骤s308,将加权矩阵w2代入式(10),即得到最终定位结果
具体的,s1的具体为:
s101,以接收机1作为参考接收机,计算辐射源i发射的信号到达接收机j的时间与其到达参考接收机的时间之差
其中,
s102,对公式(13)做变形,得到辐射源i和接收机j的rd测量公式:
考虑到实际中测量误差的影响,则进一步将辐射源i和接收机j的rd测量公式表示为:
对于系统中的n个辐射源和m个接收机,一共可以得到n(m-1)个rd测量公式;
s103,定义
则上述n(m-1)个rd测量公式表示为:
r=r°+δr(1)
其中,rd测量误差向量矩阵δr服从零均值的高斯分布,其协方差矩阵为qr=e(δrδrt)。
具体的,s201的具体步骤如下:
将辅助参数与辐射源目标位置参数之间的约束关系为
对于辐射源i,引入辅助向量
其中,
δh1,i表示为误差向量矩阵,且δh1,i=b1,iδri+d1,iδs;
b1,i为(m-1)×(m-1)的矩阵,d1,i为(m-1)×3m的矩阵,二者的具体表达式为
对于n个辐射源目标,令
g1θ°=h1+δh1(16)
式中,
δh1为总的误差向量,且δh1=b1δr+d1δs;
对式(16)进行求解,获得加权最小二乘解:
其中,w1为加权矩阵,且w1选取为:
忽略矩阵g1和d1中的误差,计算获得式(2)中θ的估计误差δθ为:
将式(17)乘以其转置,并求期望,可得θ的估计协方差矩阵为:
具体的,s202的具体步骤如下:
将辅助参数与辐射源目标位置参数之间的约束关系为
令
将
对于辐射源i,构建如下方程:
其中,
δh2,i为误差向量,且δh2,i=b2,iδθi+d2,iδs1,δθi为辅助变量θi的估计误差;
对于n个辐射源,获得n个矩阵
g2u°=h2+δh2(19)
式中,
δh2为总的误差向量,且δh2=b2δθ+d2δs,δθ为辅助变量矩阵θ的估计误差;
其中,
对式(19)进行求解,获得加权最小二乘解:
其中,w2为加权矩阵,且w2选取为:
为验证本发明技术方案的有效性,下面将推导算法定位的crlb和理论估计误差协方差矩阵,并分析二者之间的等价性:
(1)crlb
对于本文针对的接收机位置误差下多辐射源目标定位场景,记待估参量
由于rd测量误差向量δr服从零均值的高斯分布,协方差矩阵为qr;接收机位置误差向量δs服从零均值的高斯分布,其协方差矩阵为qs,并且rd测量误差和接收机位置误差相互独立,因此,给定待估参量
其中,κ是与待估参数
通过对式取对数,并对待估参量
其中,子矩阵x、y、z的详细表达式为:
根据rd测量公式:
其余位置元素为0。
根据crlb的定义可知,系统的crlb等于fim矩阵的逆,故对于式(21)中所示的系统的fim矩阵,利用分块矩阵求逆引理,得到对应的crlb为:
crlb(u°)=(x-yz-1yt)-1
=x-1+x-1y(z-ytx-1y)-1ytx-1(24)
显然,式(24)中,等号右侧的x-1即为接收机位置准确已知时目标辐射源位置估计的crlb,而x-1y(z-ytx-1y)-1ytx-1即为接收机位置误差的出现引起的目标辐射源位置估计误差的增量。
(2)理论估计协方差矩阵与crlb之间的等价性
将式(11)、式(3)和式(4)依次代入式
式中的
其中:
比较式(25)中的理论估计协方差矩阵cov(u)与式(24)中的克拉美罗界crlb(u°),可以发现两者在结构上是相同的;进一步比较式(22)(23)和式(26)(27),可以发现,在接收机位置误差和rd测量误差较小时,
也就是说,在接收机位置误差和rd测量误差较小时,
即本文所提定位方法的估计误差在接收机位置误差和rd测量误差较小时,可以达到克拉美罗界。
仿真实验
通过一组蒙特卡洛仿真实验来验证本文算法的定位精度、稳健性和算法的计算复杂度;与此同时,为了客观地评估发明定位方法的定位性能,还选用集中代表性的定位方法进行比较,分别为:chan的算法(单辐射源目标场景,未考虑接收机位置误差)、ho的算法(单辐射源目标场景,考虑接收机位置误差)、yang提出的算法(多辐射源目标场景,考虑接收机位置误差)以及sun的算法(多辐射源目标场景,考虑接收机位置误差)。
表1接收机的位置
蒙特卡洛仿真实验的场景设置如下:假设仿真场景中有m=5个接收机,其真实位置如表1所示;n=2个辐射源目标,其位置分别为:
为了模拟含有误差的rd测量,在rd真实值向量r°上添加零均值的高斯白噪声δr,δr的协方差矩阵设置为
为了模拟含有误差的接收机位置,在接收机真实位置向量s°上添加零均值的高斯白噪声δs,δs的协方差矩阵设置为
在给定rd测量误差和接收机位置误差条件下,利用算法进行蒙特卡洛仿真定位实验,并通过统计5000次蒙特卡洛仿真的均方根误差(rootmeansquareserror,rmse)来评估算法的定位误差,rmse的定义如下:
其中,u(l)表示辐射源位置向量u°在第l次蒙特卡洛仿真中的估计值。
仿真1:不同测量误差条件下算法的定位性能
首先,为了分析不同定位方法在不同rd测量误差条件下的定位性能,设定接收机位置误差水平为
如图2所示,待定位的两个辐射源对应的r1=10000,φ1=30°、r2=11000,φ2=30°,即两个辐射源位置为
利用本文的定位方法、chan的算法、ho的算法、yang的算法以及sun的算法[20]对两个辐射源进行仿真定位实验,对应的定位结果如图3所示。
图3展示了不同rd测量误差条件下算法对辐射源的定位误差。从图3可以看出,chan的算法和ho的算法在不同的rd测量误差条件下,定位rmse始终大于crlb,原因在于chan的算法没有考虑多辐射源目标场景以及接收机位置误差的影响,而ho的算法虽然考虑了接收机位置误差的影响,但是其针对单辐射源目标场景设计,在多辐射源目标场景下对目标的定位rmse无法达到crlb。
而本文的定位方法、yang的算法和sun的算法这三种算法是针对多辐射源目标场景设计的,均将接收机位置考虑到定位模型中,因此在rd测量误差较小时,这三种算法的定位rmse都可以逼近crlb。但是从局部放大子图中可以看出,随着rd测量误差增大,这三种算法的定位误差开始逐渐偏离crlb,但是三种算法开始偏离crlb的阈值不同。yang的算法和sun的算法在rd测量误差为
接下来,为了评估算法在不同接收机位置误差下的定位性能,我们固定rd测量误差为
图4给出了不同接收机位置误差条件下算法对辐射源的定位误差。在不同的接收机位置误差下,chan的算法和ho的算法的定位rmse均高于crlb。在接收机位置误差较小时,yang的算法和sun的算法都可以达到crlb,但是当接收机位置误差达到40db,yang的算法和sun的算法定位rmse开始偏离crlb,然而在这种程度的接收机位置误差下,本文算法依然可以给出crlb精度的定位结果。
仿真2:稳健性评估
如图5所示,我们令辐射源对应的角度φ1和φ2在0°~360°之间变化,对应半径设置为r1=10000、r2=11000。rd测量误差设置为
图6给出了几种算法对辐射源定位的rmse曲线,从图6可以看出,在大部分角度下,几种算法的定位rmse均较为接近crlb,说明此时几种算法的定位结果较为稳定。然而,在一些特定角度(0°、90°、180°、270°、360°)下,chan的算法[、ho的算法、yang的算法[以及sun的算法会出现rmse突然增大的异常情况。这说明上述几种算法存在不稳健性,而这种不稳健性主要是由于chan的算法、ho的算法、yang的算法以及sun的算法这四种算法的第二步加权最小二乘估计中矩阵b2在0°、90°、180°、270°、360°等角度下会出现缺秩现象,从而导致后续矩阵求逆时误差增大。相比之下,本文算法在各个角度下均具有较高的定位精度,说明本文算法的稳健性优于现有几种算法。
仿真3:计算复杂度比较
我们通过5000次蒙特卡洛仿真,统计了几种算法的平均运行时间,计算平台的配置如下:intel(r)core(tm)i7-8550u@1.80ghz;结果如表2所示。
表2不同算法计算复杂度比较
表2列出了各种算法的平均运行时间,chan的算法的计算复杂度最低,在表注所示的计算平台上运行一次的平均时间约为0.293ms。ho的算法在定位模型中考虑了接收机位置误差的影响,因此计算复杂度增加了约3倍。yang的算法、sun的算法和本文算法均针对多辐射源目标场景设计,这三种算法中,yang的算法计算复杂度最高,平均运行时间约为8.985ms,这是由于yang的算法需要同时估计目标辐射源位置和接收机位置;sun的算法计算复杂度大幅简化,平均运行时间约为0.504ms;本文算法相比于sun的算法,计算复杂度略微降低,约为0.492ms,这是由于本文算法无需额外的开方运算。
上述理论分析和仿真实验均表明,本文算法在多目标场景下的定位误差可以达到crlb,并且由于避免了矩阵缺秩和开方后的符号模糊问题,本文算法的定位精度、稳健性和计算量均优于现有算法。
最后应当说明的是:以上实施例仅用以说明本发明的技术方案而非对其限制;尽管参照较佳实施例对本发明进行了详细的说明,所属领域的普通技术人员应当理解:依然可以对本发明的具体实施方式进行修改或者对部分技术特征进行等同替换;而不脱离本发明技术方案的精神,其均应涵盖在本发明请求保护的技术方案范围当中。