最小二乘法和加权最小二乘法的非视距传播下的3个基站的二维TDOA定位性能均方误差对比，给出预测坐标并给出完整matlab代码

本问涉及到较多的知识点，因此我会逐步解释每一个概念，以便您更好地理解。首先，我们需要了解TDOA定位和最小二乘法（LS）以及加权最小二乘法（WLS）。

TDOA定位是通过计算接收信号的到达时间差来确定接收器的位置。在非视距传播下，我们通常使用三个基站来确定接收器的位置。在二维空间中，我们需要至少三个基站来解决定位问题。

最小二乘法（LS）是一种常见的数学优化方法，用于解决线性回归问题。在定位问题中，我们可以使用LS来计算估计位置。但是，LS假设所有的观测数据都是同等重要的，这在实际应用中并不一定成立。因此，为了更好地处理观测数据中的噪声，我们可以使用加权最小二乘法（WLS）。

WLS是LS的一种扩展形式，它在计算估计位置时，给每个观测数据赋予不同的权重，这意味着一些数据点会更重要，而其他数据点则会被忽略。这些权重通常基于观测数据的信噪比（SNR）或其他因素来计算。

现在，让我们来看一下预测坐标的代码实现：

% 非视距传播下的3个基站的二维TDOA定位性能均方误差对比
clc,clear,close all

%% 仿真参数设置
N = 1000; % 仿真次数
SNR = 10; % 信噪比
sigma = 0.1; % 噪声标准差
d = 500; % 基站间距离
c = 3e8; % 光速
fs = 5e6; % 采样率
T = 1/fs; % 采样时间
L = 1024; % 信号长度
t = (0:L-1)*T; % 时间序列

%% 基站和接收器坐标设置
B1 = [0,0];
B2 = [d,0];
B3 = [0,d];
B = [B1;B2;B3];
R = [1000,2000];

%% 仿真
for i = 1:N
    % 生成随机信号
    f0 = 500e3;
    f = f0 + (rand-0.5)*1e3;
    phi = 2*pi*rand;
    s = cos(2*pi*f*t + phi);
    % 信号在基站之间的传播时间
    t1 = norm(B1-R)/c;
    t2 = norm(B2-R)/c;
    t3 = norm(B3-R)/c;
    % 信号在基站之间的传播时间差
    TDOA1 = t2-t1;
    TDOA2 = t3-t1;
    % 加入高斯噪声
    n = sigma*randn(size(s));
    x1 = s + n;
    n = sigma*randn(size(s));
    x2 = s + n;
    n = sigma*randn(size(s));
    x3 = s + n;
    % 计算互相关函数
    R12 = xcorr(x1,x2);
    R23 = xcorr(x2,x3);
    R31 = xcorr(x3,x1);
    % 计算延迟
    [~,idx] = max(abs(R12));
    tau12 = (idx-L)/fs;
    [~,idx] = max(abs(R23));
    tau23 = (idx-L)/fs;
    [~,idx] = max(abs(R31));
    tau31 = (idx-L)/fs;
    % 计算TDOA
    TDOA_LS(i,:) = [tau12-tau31,tau23-tau31];
    % 计算权重
    w1 = SNR/(SNR+10*log10(var(x1)));
    w2 = SNR/(SNR+10*log10(var(x2)));
    w3 = SNR/(SNR+10*log10(var(x3)));
    % 计算加权TDOA
    TDOA_WLS(i,:) = [w1*(tau12-tau31),w2*(tau23-tau31)];
end

%% 定位
% 最小二乘法
A = [2*(B2-B1);2*(B3-B1)];
b = c^2*(TDOA_LS(:,1).^2-TDOA_LS(:,2).^2)-norm(B1-R)^2+norm(B2-R)^2;
X_LS = (A'*A)\(A'*b);
% 加权最小二乘法
A = [(B2-B1)/norm(B2-B1);(B3-B1)/norm(B3-B1)];
b = c^2*(TDOA_WLS(:,1).^2-TDOA_WLS(:,2).^2)-norm(B1-R)^2+norm(B2-R)^2;
W = diag([SNR/(SNR+10*log10(var(x1))),SNR/(SNR+10*log10(var(x2))),SNR/(SNR+10*log10(var(x3)))]);
X_WLS = (A'*W*A)\(A'*W*b);

%% 结果展示
figure
plot(B(:,1),B(:,2),'bo','LineWidth',2)
hold on
plot(R(1),R(2),'rs','LineWidth',2)
plot(X_LS(1),X_LS(2),'g^','LineWidth',2)
plot(X_WLS(1),X_WLS(2),'c*','LineWidth',2)
legend('基站','真实位置','LS预测位置','WLS预测位置')
xlabel('X/m')
ylabel('Y/m')
axis equal

%% 结果分析
MSE_LS = mean((X_LS-R').^2);
MSE_WLS = mean((X_WLS-R').^2);
fprintf('LS预测结果坐标: (%.2f,%.2f)\n',X_LS(1),X_LS(2))
fprintf('WLS预测结果坐标: (%.2f,%.2f)\n',X_WLS(1),X_WLS(2))
fprintf('LS预测结果MSE: %.4f\n',MSE_LS)
fprintf('WLS预测结果MSE: %.4f\n',MSE_WLS)

接下来，我们将分别解释代码中的各个部分：

1. 仿真参数设置

N = 1000; % 仿真次数
SNR = 10; % 信噪比
sigma = 0.1; % 噪声标准差
d = 500; % 基站间距离
c = 3e8; % 光速
fs = 5e6; % 采样率
T = 1/fs; % 采样时间
L = 1024; % 信号长度
t = (0:L-1)*T; % 时间序列

在这里，我们设置了仿真次数N，信噪比SNR，噪声标准差sigma，基站间距离d，光速c，采样率fs，信号长度L和时间序列t。

2. 基站和接收器坐标设置

B1 = [0,0];
B2 = [d,0];
B3 = [0,d];
B = [B1;B2;B3];
R = [1000,2000];

我们设置了三个基站的坐标B1、B2和B3，以及接收器的坐标R。

3. 仿真

for i = 1:N
    % 生成随机信号
    f0 = 500e3;
    f = f0 + (rand-0.5)*1e3;
    phi = 2*pi*rand;
    s = cos(2*pi*f*t + phi);
    % 信号在基站之间的传播时间
    t1 = norm(B1-R)/c;
    t2 = norm(B2-R)/c;
    t3 = norm(B3-R)/c;
    % 信号在基站之间的传播时间差
    TDOA1 = t2-t1;
    TDOA2 = t3-t1;
    % 加入高斯噪声
    n = sigma*randn(size(s));
    x1 = s + n;
    n = sigma*randn(size(s));
    x2 = s + n;
    n = sigma*randn(size(s));
    x3 = s + n;
    % 计算互相关函数
    R12 = xcorr(x1,x2);
    R23 = xcorr(x2,x3);
    R31 = xcorr(x3,x1);
    % 计算延迟
    [~,idx] = max(abs(R12));
    tau12 = (idx-L)/fs;
    [~,idx] = max(abs(R23));
    tau23 = (idx-L)/fs;
    [~,idx] = max(abs(R31));
    tau31 = (idx-L)/fs;
    % 计算TDOA
    TDOA_LS(i,:) = [tau12-tau31,tau23-tau31];
    % 计算权重
    w1 = SNR/(SNR+10*log10(var(x1)));
    w2 = SNR/(SNR+10*log10(var(x2)));
    w3 = SNR/(SNR+10*log10(var(x3)));
    % 计算加权TDOA
    TDOA_WLS(i,:) = [w1*(tau12-tau31),w2*(tau23-tau31)];
end

在这里，我们进行了N次仿真，每次生成频率为f的正弦信号s，并将其从接收器发送到三个基站。然后，我们计算了信号在基站之间的传播时间t1、t2和t3，以及信号在基站之间的传播时间差TDOA1和TDOA2。我们还为每个信号加入了高斯噪声，并计算了互相关函数R12、R23和R31。接下来，我们计算了延迟tau12、tau23和tau31，并使用它们计算了TDOA_LS和TDOA_WLS。最后，我们将所有TDOA值存储在TDOA_LS和TDOA_WLS中。

4. 定位

% 最小二乘法
A = [2*(B2-B1);2*(B3-B1)];
b = c^2*(TDOA_LS(:,1).^2-TDOA_LS(:,2).^2)-norm(B1-R)^2+norm(B2-R)^2;
X_LS = (A'*A)\(A'*b);
% 加权最小二乘法
A = [(B2-B1)/norm(B2-B1);(B3-B1)/norm(B3-B1)];
b = c^2*(TDOA_WLS(:,1).^2-TDOA_WLS(:,2).^2)-norm(B1-R)^2+norm(B2-R)^2;
W = diag([SNR/(SNR+10*log10(var(x1))),SNR/(SNR+10*log10(var(x2))),SNR/(SNR+10*log10(var(x3)))]);
X_WLS = (A'*W*A)\(A'*W*b);

在这里，我们使用LS和WLS分别计算估计位置。对于LS，我们首先计算了矩阵A和向量b，然后使用最小二乘法求解线性方程组。对于WLS，我们使用基站之间的距离和向量b计算了权重矩阵W，并使用加权最小二乘法求解线性方程组。最终，我们将估计位置存储在X_LS和X_WLS中。

5. 结果展示

figure
plot(B(:,1),B(:,2),'bo','LineWidth',2)
hold on
plot(R(1),R(2),'rs','LineWidth',2)
plot(X_LS(1),X_LS(2),'g^','LineWidth',2)
plot(X_WLS(1),X_WLS(2),'c*','LineWidth',2)
legend('基站','真实位置','LS预测位置','WLS预测位置')
xlabel('X/m')
ylabel('Y/m')
axis equal

在这里，我们用蓝色圆圈表示三个基站的位置，用红色正方形表示真实接收器位置，用绿色三角形表示LS预测位置，用青色星形表示WLS预测位置。

6. 结果分析

MSE_LS = mean((X_LS-R').^2);
MSE_WLS = mean((X_WLS-R').^2);
fprintf('LS预测结果坐标: (%.2f,%.2f)\n',X_LS(1),X_LS(2))
fprintf('WLS预测结果坐标: (%.2f,%.2f)\n',X_WLS(1),X_WLS(2))
fprintf('LS预测结果MSE: %.4f\n',MSE_LS)
fprintf('WLS预测结果MSE: %.4f\n',MSE_WLS)

在这里，我们计算了LS和WLS预测位置与真实位置之间的均方误差（MSE），并输出预测结果的坐标和MSE。

希望这个代码实现能够对您有所帮助！