|
TDOA算法原理
TDOA(Time Difference of Arrival)——时间差到达算法，利用了几何数学中双曲线的特点—— 双曲线上的任意点到达两焦点的距离差是固定值。 这个距离差它天然可以抹去用户设备(UE)和基站的之间时钟误差。 
P
1
C
1
=
c
⋅
(
t
11
+
Δ
t
)
P_1C_1 = c·(t_{11}+\Delta t)
P1C1=c⋅(t11+Δt)其中
Δ
t
\Delta t
Δt是UE和基站之间的钟差(在UE与基站不完全同步的情况下)，这个钟差我们没法直接获得。
P
1
C
2
=
c
⋅
(
t
12
+
Δ
t
)
P_1C_2 = c·(t_{12}+\Delta t)
P1C2=c⋅(t12+Δt)
则
∣
P
1
C
1
−
P
1
C
2
∣
=
c
⋅
(
t
11
−
t
12
)
则|P_1C_1-P_1C_2 |=c·(t_{11}-t_{12})
则∣P1C1−P1C2∣=c⋅(t11−t12)可见这里的钟差
Δ
t
\Delta t
Δt被消除了，之后使用数学方法求出两个双曲线的焦点。但这同时也暗示着 基站的时钟需要同步才能被消除。 所以TDOA算法特性：UE和基站无需同步，基站之间需要同步，最少三个基站能测得焦点。
Chan算法介绍
在TDOA的解算方法上，有直接求解析解的Chan算法、Fang算法。也有迭代算法如Taylor算法(它是通过不断计算当前误差来调整参数，这个误差需要真实的位置标签来对比，但我们有真实标签后为什么还需要估计呢？这个是我对Taylor算法的疑惑，欢迎大家一起探讨👏)。基于解析解的方差理论上是没有误差的，它只受限于计算机的计算精度。
Chan算法公式推导
  上述算法中，
A
x
=
r
0
C
+
D
Ax=r_0C+D
Ax=r0C+D，先求
A
x
=
C
Ax=C
Ax=C的解
x
a
x_a
xa，再求
A
x
=
D
Ax=D
Ax=D的解
x
b
x_b
xb，再求
r
0
r_0
r0，最后按照
x
=
r
0
x
a
+
x
b
x=r_0x_a+x_b
x=r0xa+xb组合起来。
Chan算法实现2D
%%
clc;
clear;
close all;
format long;
figure;
%设置UE位置
for i=1:100
ue_x = randi(100);
ue_y = randi(100);
scatter(ue_x, ue_y, '*');
hold on;
%注意注意！！！基站的y不能全部相同，否则在第57行的矩阵A第二列元素全为0，Ax=C或Ax=Ds时求不出唯一解
stations = [-40 0; 20 0; 40 50; 10 10]; %第四个基站是为了提出伪解
hold on;
r0_real = distance(ue_x, ue_y,stations(1,1), stations(1,2));
r1_real = distance(ue_x, ue_y,stations(2,1), stations(2,2));
r2_real = distance(ue_x, ue_y,stations(3,1), stations(3,2));
r3_real = distance(ue_x, ue_y,stations(4,1), stations(4,2));
%ri_real只是用来计算tds，实际上它会带有时钟误差，而这个误差我们不能直接得到
tds = [r1_real-r0_real r2_real-r0_real, r3_real-r0_real];
position = TDOA(stations, tds);
scatter(position(1), position(2), 'o');
hold on;
end
%%
function [position] = TDOA(stations, tds)
x0 = stations(1,1);y0 = stations(1,2);
x1 = stations(2,1);y1 = stations(2,2);
x2 = stations(3,1);y2 = stations(3,2);
x3 = stations(4,1);y3 = stations(4,2);
r10 = tds(1);
r20 = tds(2);
r30 = tds(3); %ue对3号基站和0号基站的距离差，真实的
scatter(x0,y0,120,'d', 'filled'); text(x0,y0,'Anchor1');
scatter(x1,y1,120,'d', 'filled'); text(x1,y1,'Anchor2');
scatter(x2,y2,120,'d', 'filled'); text(x2,y2,'Anchor3');
hold on;
x10 = x1 - x0;
x20 = x2 - x0;
y10 = y1 - y0;
y20 = y2 - y0;
k0 = x0^2 + y0^2;
k1 = x1^2 + y1^2;
k2 = x2^2 + y2^2;
A = [x10 y10; x20 y20];
C = -[r10; r20];
D = [(k1-k0-r10^2)/2; (k2-k0-r20^2)/2];
%求解Ax = r0 * C + D
a = A\C;
b = A\D;
%求解r0
A_ = a(1)^2 + a(2)^2-1;
B_ = a(1) * (b(1) - x0) + a(2) * (b(2) - y0);
C_ = (x0 - b(1))^2 + (y0 - b(2))^2;
if B_^2-A_*C_ |