定位问题基于matlab TDOA+taylor算法移动基站无源定位含Matlab源码 2098期

Posted 2022-09-07 海神之光

一、chan+taylor算法移动基站无源定位简介

1 引言
随着无人机的普及，低空空域的安全问题受到人们的极大关注．针对该问题，本研究对“非合作型”无人机采用一种基于时差法的无源定位算法对其进行实时定位．基于时差法的无源定位方法是根据求解无人机信号到达主站和各辅站的距离差，并联合各基站坐标所构成的双曲线方程组来实现．该方法定位精度高，且不对外发射信号，可在机场等区域安全使用．

目前，Chan算法和Taylor算法是2种经典的时差定位算法．其中，Chan算法在时差值精确的情况下，可以实现较高精度的定位，但如果时差值精度不够，其定位精度会大幅降低．Taylor算法则是在已有的定位坐标基础上，进行迭代递归，使定位出的坐标接近于目标的真实坐标．虽然Taylor算法定位精度较高，但需要提供初始估计坐标，否则就无法实现准确定位．基于2种算法的特点，本研究提出Chan-Taylor联合算法，其思路是，将Chan算法解算出的目标坐标作为初始估计坐标值赋给Taylor算法进行迭代运算，即使获取的时差值存在一定误差，使初始估计坐标的精度不高，但可以通过迭代来提高定位坐标的精度．通过算法对比和仿真分析表明，Chan-Taylor联合算法较Chan算法具有更高的定位精度和稳定性，较Taylor算法更具有实用性．

2 算法描述
2.1 Chan算法
基于Chan算法的无源定位是通过求解目标源信号到达辅站与主站之间的时差并联立各基站坐标所得的双曲线方程组来实现的．该算法是一种非迭代算法，不需要初始值，在时差精确、视距传输的情况下，其定位效果良好，但在工程上，很难获得满足要求的时差初值．因此，Chan算法可作为其他算法的前置条件．

本研究以4站三维定位系统为例建立3组方程，该方程组为超定方程组．通常情况下，由于该方程组导出的矩阵不存在逆矩阵，方程组无法正常求解．所以，本研究利用伪逆法联合最小二乘法对方程组进行解算，即Chan算法．4站定位系统的定位原理如图1所示．

图1 4站定位系统示意图
图1中，主站坐标联立3个辅站坐标，通过分别计算出的时差可构建3条双曲线，其交点就是无人机的位置．

假设无人机位置为P=[px,py,pz]，各基站坐标为Pi=[pix,piy,piz],i∈[0,3]，则无人机到基站i的距离ri2可表示为，

本研究若不特别指明，均默认i∈[1,3]．对式(1)整理可得，

式中，Ri表示基站i到坐标原点的距离;R0为主站到坐标原点的距离;ri0为无人机到辅站与主站间的距离差．

4站三维定位系统存在一个由3组式(2)的关系式结合而成的方程组，如式(3)所示．当A≠0时，线性方程组(3)有解．

式中，A是方程组的系数矩阵，b是方程组的输出向量．

利用伪逆法可求得无人机坐标为，

2.2 Chan-Taylor联合算法
因为Chan算法是非递归算法，对时差精度要求高，因此，本研究对该算法的定位结果进行二次处理．Taylor算法是利用局部最小二乘解进行迭代[5]的递归算法，其定位精度高，但需要初始估计坐标，否则无法进行定位．

基于低空无人机时差定位的实际需求，本研究结合Chan算法和Taylor算法提出了一种改进的算法，即Chan-Taylor联合算法．Chan-Taylor联合算法是将Chan算法的解算结果作为初始估计坐标送入Taylor算法，以达到对无人机坐标进行误差计算和定位修正的作用．算法在迭代时，将误差与设定的阈值进行比较，若误差值大于阈值，则继续迭代;若误差值小于设定阈值，则终止迭代并输出结果．

2.3 Chan-Taylor联合算法流程及计算原理
Chan-Taylor联合算法流程如图2所示，具体为:首先，算法获得无人机信号到主站与各辅站之间的时差;然后，将时差用于Chan算法部分进行初始估计坐标值的计算，并利用该坐标值在Taylor算法部分做误差向量的计算，用以定位修正;同时，对误差进行阈值比较．如不满足条件，则继续迭代，如满足条件，则结束迭代，并输出最终结果．

图2 Chan-Taylor算法流程

式中，α为目标差值向量，b为时差的差值向量，e为时差估计误差向量，H为时差估计的梯度矩阵．它们可分别表示为，

由式ri0=ri-r0=cτi0与站址坐标，可得，

将式(10)与式(11)联立，化简可得，

二、部分源代码

clear all
clc
close all
tic
X = 150;
Y = 150; %空间范围大小
N = 200; %仿真点数
c = 3e8;
num = 0;
distribute = []; %用来统计不同误差范围内的点数
num0 = 0;
num1 = 0;
num2 = 0;
num3 = 0;
num4 = 0;
num5 = 0;

% *直角三角形布局
o1x = 0;o1y = 0; %坐标原点
o2x = X;o2y = 0; %x轴方向
o3x = 0;o3y = Y; %y轴方向
o4x = X;o4y = Y; %辅助站
MS(1,1) =25;
MS(1,2) = 25;
% 算法开始：
BS = [0, X, 0,X
0, 0, Y,Y]; % 参考站坐标
% TDOA协方差矩阵Q：
BSN=4;
Noise =rand(1); %测距误差1.5厘米
figure(1)

plot(o2x,o2y,‘rs’,‘LineWidth’,2,…
‘MarkerEdgeColor’,‘k’,…
‘MarkerFaceColor’,‘r’,…
‘MarkerSize’,10);hold on
plot(o3x,o3y,‘rs’,‘LineWidth’,2,…
‘MarkerEdgeColor’,‘k’,…
‘MarkerFaceColor’,‘r’,…
‘MarkerSize’,10);hold on
plot(o4x,o4y,‘rs’,‘LineWidth’,2,…
‘MarkerEdgeColor’,‘k’,…
‘MarkerFaceColor’,‘r’,…
‘MarkerSize’,10);hold on
grid on
axis ([-500 500 -500 500]);

for j=1:N
mse = 0;

XX = 25;
YY =25;
o1x=o1x-0.01*N;
o1y=o1y-0.01*N;
user = [XX YY];
plot(o1x,o1y,'rs','LineWidth',2,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor','r',...
    'MarkerSize',10);hold on
%****************************解方程组1(1 2 3)********************************
%*****************************************************************
beta21 = o2x^2+o2y^2-o1x^2-o1y^2;
beta31 = o3x^2+o3y^2-o1x^2-o1y^2;

D1 = sqrt((user(1)-o1x)^2+(user(2)-o1y)^2);
D2 = sqrt((user(1)-o2x)^2+(user(2)-o2y)^2);
D3 = sqrt((user(1)-o3x)^2+(user(2)-o3y)^2);

D12 = D1-D2;
D13 = D1-D3;
D32 = D3-D2;

%   tdoa_bias = 5e-11;

tdoa_bias=(5e-11)*(2*round(rand)-1)
%**************TDOA偏差0.05ns******************************
t12 = D12/c + tdoa_bias;
t13 = D13/c + tdoa_bias;
t32 = D32/c + tdoa_bias;
a1 = t12*(o3x-o1x)-t13*(o2x-o1x);
b1 = t12*(o3y-o1y)-t13*(o2y-o1y);
g1 = 1/2*(c^2*t12*t13*t32+t12*beta31-t13*beta21);
A = -b1/a1;
B = g1/a1;

E = (1/(c*t12))*((o2x-o1x)*A+(o2y-o1y));
F = (c*t12)/2+(1/(2*c*t12))*(2*(o2x-o1x)*B-beta21);

G = A^2-E^2+1;
H = 2*(A*(B-o1x)-o1y-E*F);
I = (B-o1x)^2+o1y^2-F^2;

y11 = -H/(2*G)+sqrt((H/(2*G))^2-I/G);
y12 = -H/(2*G)-sqrt((H/(2*G))^2-I/G);

x11 = A*y11+B;
x12 = A*y12+B;

mse11 = sqrt((x11-user(1))^2+(y11-user(2))^2);
mse12 = sqrt((x12-user(1))^2+(y12-user(2))^2);
position11 = [x11 y11];
position12 = [x12 y12];
    %**********************判断合理解****************************
%方程组1得到的用户坐标到参考站1与到参考站2的距离差
distance_diff11 = sqrt((position11(1))^2+(position11(2))^2)-sqrt((position11(1)-50)^2+(position11(2))^2);
distance_diff12 = sqrt((position12(1))^2+(position12(2))^2)-sqrt((position12(1)-50)^2+(position12(2))^2);

if(t12==0)
    if(t13==0)
        position = [25 25];
        mse = 0;
        k = 0;
    else    %用户在竖中垂线上，只能用t42作判断
        k = 1; %这种情况下第二组方程组有解，直接去判断第二组中的合理解。
    end
elseif(t13==0)
    if(t12==0)
        position = [25 25];
        mse = 0;
    else    %用户在横中垂线上，只能用t43作判断
        k = 2;  %这种情况下两组方程均无解。。。。。。。。。。另想办法，再加一组解？
    end               
else%必须要有下面的先通过坐标判断，因为有distance_diff11和distance_diff12相等的情况，如点（5，4），偏差3时。
    if((x11>50)||(x11<-0.015)||(y11>50.015)||(y11<-0.015))
        position = position12;
        mse = mse12;
    elseif((x12>50)||(x12<-0.015)||(y12>50.015)||(y12<-0.015))
        position = position11;
        mse = mse11;
    elseif(t12>0)
        if(distance_diff11>distance_diff12)
            position = position11;
            mse = mse11;
        elseif(distance_diff12>distance_diff11)
            position = position12;
            mse = mse12;
        end
    elseif(t12<0)
        if(distance_diff11<distance_diff12)
            position = position11;
            mse = mse11;
        elseif(distance_diff12<distance_diff11)
            position = position12;
            mse = mse12;
        end
    end
    k = 3;%得到定位结果后。为了判断是否在大误差范围内而设的标志位。
end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 罗正华,雷林,周方均,李霞.基于Chan-Taylor联合算法的低空无人机时差定位研究[J].成都大学学报(自然科学版). 2019,38(04)

3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除