电子科技大学航空航天学院导航原理基础课程利用广播星历计算GPS卫星位置坐标。

rain_crying

于 2024-04-03 11:25:04 发布

阅读量323

点赞数 8

文章标签： matlab 大作业学习

本文链接：https://blog.csdn.net/rain_crying/article/details/137337311

版权

本文详细介绍了如何根据给定的GPS星历参数，通过数据预处理和一系列计算公式，求解卫星的位置坐标。作者使用Python代码实现了从卫星发射时刻到计算卫星位置的过程，并提到与老师答案的误差在分米级别，符合课程要求。

摘要由CSDN通过智能技术生成

一.根据所给的初始化数据进行数据预处理

PS：每一届老师所给的参数可能不一致，建议大家自行修改

%% 数据预处理
clear;
clc;
 
% 卫星发射时刻
t = 535460.931374331;
 
% GPS星历参数
toc = 3.528000000000e+05;
sqrtA = 5.15378158e+03;
toe = 532800;
deltan = 4.07552691e-09;
M0 = -0.95904256173;
e = 0.001125632785;
i0 = 0.96091568;
IDOT = -1.3786289e-10;
omega = -3.007907646;
Cuc = -2.84984708e-06;
Cus = 1.289881766e-05;
Crc = 1.3146875e+02;
Crs = -53.59375;
Cic = 3.725290299e-08;
Cis = -6.70552254e-08;
omega0 = 0.056568024125;
omegaDot = -7.7674664e-09;
%常用数据
omegadote = 7.2921151467e-5;%自传角速度
miu= 3.986005e14;%μ

二.根据所给公式进行计算

公式具体内容博主记不太清了，详细历程请见CSDN博主文章，基本历程几乎相同，课程仅仅要求到计算位置坐标，无需进行下属步骤。使用广播星历进行 GPS 卫星位置的计算_gps 广播星历 toe-CSDN博客、

代码如下：

%% 计算GPS卫星位置

%计算归化时间                          
tk = t - toe;

%计算平近点角
A = sqrtA^2;
n = sqrt(miu/(A^3)) + deltan;
M = M0 + n*tk;


%计算偏近点角[迭代到1e-5时候改变不大，故采用最大值]
E0 = 0;
while 1
    E_k = M+e*sin(E0);
    delt_E = E_k - E0;
    delt_E = rem(delt_E,2*pi);
    if abs(delt_E) < 1e-5
        break;
    else
        E0 = E_k;
    end
end
 

%计算真近点角和维度幅角
u_k = atan2(sqrt(1-e^2)*sin(E_k),cos(E_k)-e);
phi = u_k + omega;

%计算摄动修正后的升交点角距、轨道半径、轨道倾角
delt_u = Cuc*cos(2*phi)+Cus*sin(2*phi);
delt_r = Crc*cos(2*phi)+Crs*sin(2*phi);
delt_i = Cic*cos(2*phi)+Cis*sin(2*phi);
u_k = phi + delt_u;
r = A*(1-e*cos(E_k))+delt_r;
i = i0 + IDOT*tk + delt_i;

%计算卫星轨道平面的坐标
x_k=r*cos(u_k);
y_k=r*sin(u_k);


%计算升交点赤经
OMEGA = omega0 + (omegaDot-omegadote)*tk - omegadote*toe;

%卫星坐标
positionX = cos(u_k)*r*cos(OMEGA) - sin(u_k)*r*cos(i)*sin(OMEGA);
positionY = cos(u_k)*r*sin(OMEGA) + sin(u_k)*r*cos(i)*cos(OMEGA);
positionZ = sin(u_k)*r*sin(i);



% 命令行打印结果
fprintf('卫星位置：[%f, %f, %f]\n',positionX,positionY,positionZ);

根据课程教师所给参数以及相关公式计算，博主算得的结果为：

卫星位置：[-438453.887492, 24873592.589880, 9235528.251340]

与老师的答案误差为分米级，达到计算要求。

三.总结

综上所述，综合代码如下：

%% 数据预处理
clear;
clc;
 
% 卫星发射时刻
t = 535460.931374331;
 
% GPS星历参数
toc = 3.528000000000e+05;
sqrtA = 5.15378158e+03;
toe = 532800;
deltan = 4.07552691e-09;
M0 = -0.95904256173;
e = 0.001125632785;
i0 = 0.96091568;
IDOT = -1.3786289e-10;
omega = -3.007907646;
Cuc = -2.84984708e-06;
Cus = 1.289881766e-05;
Crc = 1.3146875e+02;
Crs = -53.59375;
Cic = 3.725290299e-08;
Cis = -6.70552254e-08;
omega0 = 0.056568024125;
omegaDot = -7.7674664e-09;
%常用数据
omegadote = 7.2921151467e-5;%自传角速度
miu= 3.986005e14;%μ

%% 计算GPS卫星位置

%计算归化时间                          
tk = t - toe;

%计算平近点角
A = sqrtA^2;
n = sqrt(miu/(A^3)) + deltan;
M = M0 + n*tk;


%计算偏近点角[迭代到1e-5时候改变不大，故采用最大值]
E0 = 0;
while 1
    E_k = M+e*sin(E0);
    delt_E = E_k - E0;
    delt_E = rem(delt_E,2*pi);
    if abs(delt_E) < 1e-5
        break;
    else
        E0 = E_k;
    end
end
 

%计算真近点角和维度幅角
u_k = atan2(sqrt(1-e^2)*sin(E_k),cos(E_k)-e);
phi = u_k + omega;

%计算摄动修正后的升交点角距、轨道半径、轨道倾角
delt_u = Cuc*cos(2*phi)+Cus*sin(2*phi);
delt_r = Crc*cos(2*phi)+Crs*sin(2*phi);
delt_i = Cic*cos(2*phi)+Cis*sin(2*phi);
u_k = phi + delt_u;
r = A*(1-e*cos(E_k))+delt_r;
i = i0 + IDOT*tk + delt_i;

%计算卫星轨道平面的坐标
x_k=r*cos(u_k); 
y_k=r*sin(u_k);


%计算升交点赤经
OMEGA = omega0 + (omegaDot-omegadote)*tk - omegadote*toe;

%卫星坐标
positionX = cos(u_k)*r*cos(OMEGA) - sin(u_k)*r*cos(i)*sin(OMEGA);
positionY = cos(u_k)*r*sin(OMEGA) + sin(u_k)*r*cos(i)*cos(OMEGA);
positionZ = sin(u_k)*r*sin(i);



% 命令行打印结果
fprintf('卫星位置：[%f, %f, %f]\n',positionX,positionY,positionZ);