%伪距定位算法(多星)
%距离单位为km
%----------------------------begin-----------------------------------------
clear;clc;
c=299792.458;%光速
%通过星历参数解算到所在地可见卫星的坐标位置
sat13=[-7134.52924416113.64883623709.205570];
sat22=[-22383.70004018533.2331685307.245613];
sat23=[-5384.90131728971.6223232079.796362];
sat14=[637.46657128016.0538419347.297933];
sat12=[-11568.199533-3328.51154326977.312423];
sat21=[-28908.916747-577.0617606051.375658];
sat5=[-1205.65118128296.890128-8397.025036];
sat4=[16456.52732412347.28249421199.173063];
sat=[sat13;sat22;sat23; sat14; sat12; sat21; sat5; sat4];%多卫星位置矩阵
%所在地实际大地坐标,用来与定位结果作比较
nanjing=[-2604.2985334743.2972173364.978513];%理想伪距测量值
r0 =zeros(8,1);for i =1:8r0(i,1)=norm(sat(i,:)-nanjing);
end
%加入零均值,方差为80的随机高斯分布,模拟含有误差的伪距r
r = r0 +sqrt(80)*randn(size(r0))*1e-3;%--------------------------------------------------------------------------%Newton迭代法
%设定迭代初值,若无法估计则全部假设为0
xyzt=[0000];for i =1:100% 最大迭代次数设为100次
f =dingwei_fun(xyzt,sat,r);
df =dingwei_dfun(xyzt,sat);%左除,具有良好的数值稳定性,在MATLAB中此已经为最小二乘意义下的解
%delta(xyzt)=G^(-1)*b
delta=-df\f;xyzt(1)=xyzt(1)+delta(1);xyzt(2)=xyzt(2)+delta(2);xyzt(3)=xyzt(3)+delta(3);xyzt(4)=xyzt(4)+delta(4);
p=norm(delta)%定位精度
if(p<1e-10)break;
end
end
disp('迭代次数为')
i
disp('用户位置为')xyzt(1:3)disp('用户钟差为')xyzt(4)%--------------------------------end---------------------------------------
二、子主函数
function f=dingwei_fun(xyzt,sat,r)%多卫星定位方程函数
%xyzt(1:3)为用户位置(km)%xyzt(4)为用户钟差(s)%r为测量得到的伪距(km)%sat为卫星数据
c=299792.458;%光速
[m,n]=size(sat);for i=1:m
f(i)=norm(sat(i,:)-xyzt(1:3))+c*xyzt(4)-r(i);
end
f=f';
%多颗卫星定位的偏导数矩阵
function df=dingwei_dfun(xyzt,sat)%xyzt为用户位置及钟差
%sat为卫星数据
c=299792.458;%光速
[m,n]=size(sat);
xyz=xyzt(1:3);for i=1:m
for j=1:3df(i,j)=(xyz(j)-sat(i,j))/norm(sat(i,:)-xyz(j));
end
end
df(:,4)=c;%线性函数ct,对t求偏导为c