% 交叉定位 - 最小二乘法定位算法模拟

% 参数设置
numIterations = 1000; % 模拟迭代次数
maxDistance = 1000; % 最远定位距离（设定范围）
speedOfSound = 343; % 声速（单位：m/s）

% 预警机坐标
source = [0, 0];

% 初始化结果
crbResults = zeros(numIterations, 1);
rmseResults = zeros(numIterations, 1);

% 模拟迭代
for i = 1:numIterations
% 随机生成无人机坐标
drone = generateRandomPosition(maxDistance);

% 计算距离
distance1 = calculateDistance(source, drone);
distance2 = calculateDistance(source, drone);% 添加测量误差
measurement1 = distance1 + normrnd(0, 1);
measurement2 = distance2 + normrnd(0, 1);% 定位算法 - 最小二乘法
estimatedPosition = leastSquaresLocalization(source, measurement1, measurement2, speedOfSound);% 计算CRB
crb = calculateCRB(source, drone, speedOfSound);
crbResults(i) = crb;% 计算RMSE
rmse = norm(estimatedPosition - drone);
rmseResults(i) = rmse;

end

% 寻找最远定位距离
maxDistanceIdx = find(rmseResults == max(rmseResults));
maxDistanceValue = sqrt(crbResults(maxDistanceIdx));

% 显示结果
fprintf(‘最远定位距离：%.2f m\n’, maxDistanceValue);

% 生成随机位置
function position = generateRandomPosition(maxDistance)
angle = rand * 2 * pi;
distance = rand * maxDistance;
position = distance * [cos(angle), sin(angle)];
end

% 计算距离
function distance = calculateDistance(source, target)
distance = norm(target - source);
end

% 最小二乘法定位算法
function estimatedPosition = leastSquaresLocalization(source, measurement1, measurement2, speedOfSound)
A = 2 * [source(1) - measurement1(1), source(2) - measurement1(2); …
source(1) - measurement2(1), source(2) - measurement2(2)];
b = [measurement1(1)^2 - source(1)^2 + measurement1(2)^2 - source(2)^2 - speedOfSound^2 * measurement1(3)^2; …
measurement2(1)^2 - source(1)^2 + measurement2(2)^2 - source(2)^2 - speedOfSound^2 * measurement2(3)^2];
estimatedPosition = (A’ * A) \ (A’ * b);
end

% 计算CRB（Cramér-Rao下界）
function crb = calculateCRB(source, target, speedOfSound)
distance = norm(target - source);
crb = (speedOfSound^2 / (4 * pi^2)) * (1 / distance)^2;
end