一、理论基础

无人机(UAV)因其体积小,灵活性高,成本低等优势得到快速发展并被广泛应用于军事战争,城市管理,民用,地质,抢险救灾等各个领域,与此同时,无人机定位技术也得到了深入研究,其中无线电探测与定位技术备受众多学者关注.基于电磁波参数的无线电探测与定位技术主要包括角度估计,信号到达时间差估计(TDOA),接收信号强度估计(RSS),三种方法均可实现独立定位.信号到达时间差估计要求时钟同步,接收信号强度估计对信道环境比较敏感,而针对角度估计的联合定位算法普遍需要发射阵列和接收阵列角度配对.本文旨在提高无人机定位的估计精确度和降低联合估计过程中的运算复杂度,围绕多输入多输出(MIMO)双基地雷达下波离方向(DOD)和波达方向(DOA)联合估计技术开展研究,提出消除配对过程的联合估计算法,并结合角度和信号到达时间差实现无人机的位置定位.

首先介绍下TDOA的概念，如图所示，假设我们在空间中有一个声源(记为s(t)，其在空间的位置为S)两个麦克风(记为m1和m2，它们在空间的位置分别为M1和M2，接收到的信号为x1(t)和x2(t)

那么麦克风m1和m2收到的信号分别为：

其中τ1和τ2分别是声源到达两个麦克风的延迟时间，n1(t)和n2(t)为加性噪声。那么声源信号到达两个麦克风的TDOA为

其中c是声速。一般情况下，我们选择一个麦克风的信号作为参考信号，例如我们把M2作为参考信号，那么τ2=0。在麦克风阵列几何形状已知的情况下，声源定位问题变为对时延的估计问题。

DOA定位技术原理是利用接收机处的阵列天线和波达方向（DOA）估计技术来确定一个从接收机到信源的波达方向线，即为方向线（LOB），最后利用多个接收机估计的DOA进行三角测量，方向线的交点就是信源的估计位置。

波达方向（Direction Of Arrival）估计，又称为角谱估计（Angle spectral estimation）、波达角（Angle Of Arrival）估计。一个信源有很多可能的传播路径和到达角。如果几个发射机同时工作，每个信源在接收机处形成潜在的多径分量。因此，接收天线能估计出这些到达角就显得很重要，目的是估计出哪个发射机在工作以及发射机所处的方向，简单的说就是利用己方雷达接收来自目标发射机的来波方向进行估计；其物理原理是利用电磁波的直线传播原理。
通过测量辐射信号的波达方向（DOA全称Direction Of Arrival）或波达角（AOA）来估测辐射源的位置。理论上这种估计只需要两个接收阵元就可以确定辐射源的位置，但在实际中，由于受到角度分辨率、多径和噪声限制，所需阵元通常要多于两个。

二、核心程序

%-----------------------------------------------------------------------
%  --- 外辐射源基于DOA联合TDOA时间积累下二维平面GDOP分析 ---一发一收体制----
clc;
clear;
close all;
warning off;
addpath(genpath(pwd));

c=3e8;                         % 传播速度

 

% % 市区
xo=-0.25e3;yo=0;                 % 接收站1的位置   %市区 
xa=0.25e3;ya=0;                  % 发射站1的位置   %市区 
% 目标位置
xt=-2000:53:3000;
yt=-2000:53:2000;

da=pi/180;                     % 方位测量误差标准差
dtao=1e-7;                     % 时差测量误差标准差

R=[(da)^2 0;0 (dtao)^2];       % 测量误差协方差矩阵



for l=1:length(xt)
    for j=1:length(yt)
        
        a=(sqrt((xt(l)-xo)*(xt(l)-xo)+(yt(j)-yo)*(yt(j)-yo))+sqrt((xt(l)-xa)*(xt(l)-xa)+(yt(j)-ya)*(yt(j)-ya)))/2;
        tao=2*(a-xo)/c;
        
        afa_diff_x =-(yt(j)-yo)/(xt(l)-xo)^2/(1+(yt(j)-yo)^2/(xt(l)-xo)^2);
        afa_diff_y =1/(xt(l)-xo)/(1+(yt(j)-yo)^2/(xt(l)-xo)^2);
        
        f_diff_x =2*xt(l)/(xo+c*tao)^2;
        f_diff_y =2*yt(j)/(1/4*c^2*tao^2+c*tao*xo);
        f_diff_tao =-2*xt(l)^2/(xo+c*tao)^3*c-yt(j)^2/(1/4*c^2*tao^2+c*tao*xo)^2*(1/2*c^2*tao+c*xo);
        
        H=[afa_diff_x afa_diff_y;-f_diff_x/f_diff_tao -f_diff_y/f_diff_tao];
        
        Px=pinv(H)*R*pinv(H)';
        
        GDOP(l,j)=sqrt(Px(1,1)+Px(2,2));
        
    end
end

 
V=[0.005 0.007 0.01 0.015 0.02 0.025  0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.08 0.10 0.15 0.2 0.3 0.5 1];%市区 %1e-7
figure(1);
[pic1]=contour(xt./1000,yt./1000,GDOP.'/1000,V);
clabel(pic1);xlabel('X(Km) ');ylabel('Y(Km) ');title( 'CONTOUR of valleys of GDOP in T_1R_1 (Km)' );
hold on;
plot(xo/1000,yo/1000,'gP',xa./1000,ya./1000,'rP');
axis([-1.25,1.75,-1.5,1.5]);%市区
% axis([-30,50,-40,40]);%郊区
grid on