基于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);