基于MATLAB的GPS卫星绕地运行轨迹动态模拟仿真
1.算法运行效果图预览
2.算法运行软件版本
matlab2022a
3.算法理论概述
模拟24颗GPS卫星的轨道运行是一个复杂的任务，涉及到多个卫星的轨道计算和绘制。以下是一个大致的步骤和示例代码，用于在MATLAB中模拟和绘制这些卫星的轨道运行。实际GPS卫星的轨道参数非常复杂，而且卫星之间的相对运动也需要考虑。下面的示例代码是一个简化的版本，用于演示基本的概念。要通过MATLAB模拟24个GPS卫星的轨道运行效果，你可以按照以下步骤进行操作：
获取卫星轨道参数： GPS卫星的轨道参数可以从相关文献或数据源中获取。这些参数包括每颗卫星的半长轴、偏心率、轨道倾角、升交点赤经等。
计算卫星的轨道： 使用获取的轨道参数，可以通过开普勒运动方程计算每颗卫星在每个时间点的位置和速度。
选择模拟时间范围： 选择一个适当的时间范围，以便观察卫星在地球上的运动。
生成时间序列： 生成一系列时间点，可以使用MATLAB中的时间函数，如linspace来生成等间隔的时间点。
计算卫星位置： 对于每个时间点，使用计算得到的轨道参数和时间，计算每颗卫星的位置和速度。
绘制卫星轨道： 使用MATLAB的绘图函数，如plot3，绘制每颗卫星在三维空间中的轨道。你可以在地球球面上绘制卫星的位置，也可以绘制在三维笛卡尔坐标系中的轨迹。
通过使用GPS卫星星历（Almanac data）信息，来计算模拟24个GPS卫星的轨道。每个卫星用PRN1-24来编号，假设GPS卫星轨道是圆的。
自己卫星轨道模拟，这个卫星是离地面350Km的太阳同步轨道卫星。这个轨道是椭圆的，使用轨道倾角98度
计算自己卫星和每个GPS卫星的多普勒频移。
注意：这里使用GPS卫星是发送源，自己的卫星是接收源。假设发送频率1575.42 MHz
GPS卫星轨道周期几乎是24小时，而自己的卫星在太阳同步轨道上的周期大概是1.5个小时，那么就是说太阳同步轨道已经绕几周了，GPS卫星才饶一周。所以当算多普勒频移的时候只需要算出GPS一个周期时间内的多普勒频移就好了。就是说，如果在算多普勒频移的时候，如果算多过24小时，那么多普勒频移就会重复了。我只需要24小时GPS轨道周期内的多普勒频移就好了。
这里，首先介绍一下星历文件的含义：
Prn
卫星编号
iode
电文中给出的当前参考历元的有效期
Crs
电文中给出的轨道半径角距的改正项—正弦振幅
delta_n
电文中给出的平地点角改正值
M_zero
电文中给出的参考时刻平近点角
Cuc
电文中给出的升交点赤经的改正项—余弦振幅
e1
电文中给出的轨道椭圆偏心率
Cus
电文中给出的升交点赤经的改正项—正弦振幅
sqrt_a
电文中给出的卫星轨道椭圆长半轴的平方根
toe
电文中给出的参考时刻
Cic
电文中给出的倾角角距的改正项—余弦振幅
OMEGA_zero
电文中给出的参考时刻升交点赤经
Cis
电文中给出的倾角角距的改正项—正弦振幅
i_zero
电文中给出的参考时刻轨道倾角
Crc
电文中给出的轨道半径角距的改正项—余弦振幅
omega
电文中给出的轨道近地点角距
OMEGA_dot
电文中给出的升交点赤经变化率
i_dot
电文中给出的轨道倾角变化率
这里需要注意的时候，由于GPS距离地面的高度一般为20000km，而这里的同步卫星只有350km，所以看上去会效果不明显，所以这里我们把这里的参数设置的大些，这样看上去效果稍微明显点。然后你再写论文的时候，如果用到其中的数据，只要把他改回350即可。另外，其周期为1.5小时，这样在房子的时候，速度太快，不容易观察，这里稍微设置的大些，使用周期为6小时。
4.部分核心程序
Prn = NavData(PRNS_SEL,1);%识别导航数据中的PRN
iode = NavData(PRNS_SEL,11);%企业日期序列号
Crs = NavData(PRNS_SEL,12);%轨道半径正弦谐波校正的振幅
delta_n = NavData(PRNS_SEL,13);%与计算值的平均运动差
M_zero = NavData(PRNS_SEL,14);%参考时间的平均异常
Cuc = NavData(PRNS_SEL,15);%纬度变元余弦谐波校正项的振幅
es = NavData(PRNS_SEL,16);%偏心率
Cus = NavData(PRNS_SEL,17);%纬度变元正弦谐波校正项的振幅
sqrt_a = NavData(PRNS_SEL,18);%半长轴的平方根
toe = NavData(PRNS_SEL,19);%星历重新出现时间（GPS周第二次）
Cic = NavData(PRNS_SEL,20);%倾角余弦谐波校正项的振幅
OMEGA_zero = NavData(PRNS_SEL,21);%参考时间赤经
Cis = NavData(PRNS_SEL,22);%倾斜角正弦谐波校正项的振幅
i_zero = NavData(PRNS_SEL,23);%参考时间的倾角
Crc = NavData(PRNS_SEL,24);%轨道半径余弦谐波校正的幅度
omega = NavData(PRNS_SEL,25);%近地点论点
OMEGA_dot = NavData(PRNS_SEL,26);%赤经角变化率
i_dot = NavData(PRNS_SEL,27);%倾斜度变化率
%计算真实异常
%校正平均运动
n_initial = sqrt(Gravitational_constant)/(sqrt_a*sqrt_a*sqrt_a);
n = n_initial + delta_n;
%计算自参考历元以来的时间
t_k = func_tk_limits(t,toe);
%T_K平均异常
M_k = M_zero + n * t_k + 2*pi;
E_k = M_k;
E_k1= M_k - es*sin(E_k);
%E_K的迭代求解
while (abs(E_k - E_k1) > 1e-9)
E_k1 = E_k;
E_k = M_k + es * sin(E_k1);
end;
%TRUE ANOMALY
v_k = 2*atan(sqrt((1+es)/(1-es))*tan(E_k/2));
%%纬度论
%未修正纬度的自变量
phi_k = v_k + omega;
%校正
C_u = Cus*sin(2*phi_k) + Cuc*cos(2*phi_k);
C_r = Crs*sin(2*phi_k) + Crc*cos(2*phi_k);
C_i = Cis*sin(2*phi_k) + Cic*cos(2*phi_k);
