目录

一、理论基础

二、核心程序

三、测试结果

一、理论基础

1. 通过使用GPS卫星星历（Almanac data）信息，来计算模拟24个GPS卫星的轨道。每个卫星用PRN1-24来编号，假设GPS卫星轨道是圆的。
2. 自己卫星轨道模拟，这个卫星是离地面350Km的太阳同步轨道卫星。这个轨道是椭圆的，使用轨道倾角98度
3. 计算自己卫星和每个GPS卫星的多普勒频移。

注意：这里使用GPS卫星是发送源，自己的卫星是接收源。假设发送频率1575.42 MHz

        GPS卫星轨道周期几乎是24小时，而自己的卫星在太阳同步轨道上的周期大概是1.5个小时，那么就是说太阳同步轨道已经绕几周了，GPS卫星才饶一周。所以当算多普勒频移的时候只需要算出GPS一个周期时间内的多普勒频移就好了。就是说，如果在算多普勒频移的时候，如果算多过24小时，那么多普勒频移就会重复了。我只需要24小时GPS轨道周期内的多普勒频移就好了。

二、核心程序

....................................................................
        %计算GPS卫星轨道
        %计算GPS卫星轨道
        Prn        = NavData(PRNS_SEL,1);
        %电文中给出的当前参考历元的有效期
        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);
        %电文中给出的倾角角距的改正项—余弦振幅
        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);  
..........................................................................

        %计算卫星在平面坐标系中的位置    
        X_k = x_k*cos(OMEGA_k) - y_k*cos(i_k)*sin(OMEGA_k);
        Y_k = x_k*sin(OMEGA_k) + y_k*cos(i_k)*cos(OMEGA_k);
        Z_k = y_k*sin(i_k);     
        %实际空间轨道坐标
        A1=[32.8,92.8,152.8,212.8,272.8,332.8];
        for k=1:GDNUM
            A(k) = A1(k)*pi/180;
            B(k) = 55*pi/180;
            C(k) = pi/100;
            R1= [1       0          0;
                 0       cos(B(k)) -sin(B(k));
                 0       sin(B(k))  cos(B(k))]; 
            R2= [cos(C(k)) -sin(C(k))  0;
                 sin(C(k))  cos(C(k))  0; 
                 0          0          1];
            R3= [cos(A(k)) -sin(A(k))  0;
                 sin(A(k))  cos(A(k))  0;
                 0          0          1];
            R312    = R3*R1*R2;
            Ans     = R312*[x;y;z];
            x1(k,j) = Ans(1,:);
            y1(k,j) = Ans(2,:);
            z1(k,j) = Ans(3,:);        
        end
        %计算信号发射时刻的E_k_dot  
        M_k_dot = n;
        E_k_dot = M_k_dot/[1 - es*cos(E_k)];
        %计算信号发射时刻的phi_k_dot
        v_k_dot   = sqrt(1-es*es)*E_k_dot/(1-es*cos(E_k));
        phi_k_dot = v_k_dot;    
.........................................................................
        
        %计算每个时刻的频偏
        %计算每个时刻的频偏
        if j == 1
           Fx_GPS(i,j)  = 0;
           Fy_GPS(i,j)  = 0;
           Fz_GPS(i,j)  = 0;  
        else
           Fx_GPS(i,j)  = 1575.42e6/3e8*abs(Vx_GPS(i,j)-Vx_GPS(i,j-1));
           Fy_GPS(i,j)  = 1575.42e6/3e8*abs(Vy_GPS(i,j)-Vy_GPS(i,j-1));
           Fz_GPS(i,j)  = 1575.42e6/3e8*abs(Vz_GPS(i,j)-Vz_GPS(i,j-1));  
        end
    end
end
for i = 1:24
    for j = 1:144 
        Fre{i}(j,:) = [Fx_GPS(i,j), Fy_GPS(i,j), Fz_GPS(i,j)];
    end
end


%显示卫星轨迹，动态显示
%显示卫星轨迹，动态显示
L  = length(x1);
X1 = x1;
Y1 = y1;
Z1 = z1;

X2 = [x1(:,round(1*L/4)+1:end),x1(:,1:round(1*L/4))];
Y2 = [y1(:,round(1*L/4)+1:end),y1(:,1:round(1*L/4))];
Z2 = [z1(:,round(1*L/4)+1:end),z1(:,1:round(1*L/4))];

X3 = [x1(:,round(2*L/4)+1:end),x1(:,1:round(2*L/4))];
Y3 = [y1(:,round(2*L/4)+1:end),y1(:,1:round(2*L/4))];
Z3 = [z1(:,round(2*L/4)+1:end),z1(:,1:round(2*L/4))];

X4 = [x1(:,round(3*L/4)+1:end),x1(:,1:round(3*L/4))];
Y4 = [y1(:,round(3*L/4)+1:end),y1(:,1:round(3*L/4))];
Z4 = [z1(:,round(3*L/4)+1:end),z1(:,1:round(3*L/4))];
colortable;
figure;
%根据计算得到的速度，模拟动态的卫星运动效果
for time = 1:Days
    for i=1:GDNUM
        plot3(x1(i,:),y1(i,:),z1(i,:),colors{i});
        hold on
        plot3(X1(i,time),Y1(i,time),Z1(i,time),colors1{i});
        hold on
        plot3(X2(i,time),Y2(i,time),Z2(i,time),colors2{i});
        hold on        
        plot3(X3(i,time),Y3(i,time),Z3(i,time),colors3{i});
        hold on 
        plot3(X4(i,time),Y4(i,time),Z4(i,time),colors4{i});
        hold on   
    end
    plot3(Px_MY(:),Py_MY(:),Pz_MY(:),'r','linewidth',2);
    hold on   
    plot3(Px_MY(time),Py_MY(time),Pz_MY(time),'--rs','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5);     
    hold on 
    func_earth();
    pause(0.01);
    axis([-30e6,30e6,-30e6,30e6,-30e6,30e6]);
    xlabel('X');
    xlabel('Y');
    xlabel('Z');
    title('GPS卫星动态模拟效果');
    grid on;
    hold off;
end

%计算各个时刻各个卫星的频移,并以表格形式输出
%计算各个时刻各个卫星的频移,并以表格形式输出
%A01-82

三、测试结果

在matlab2021a中仿真得到如下的效果：

         这里，首先介绍一下星历文件的含义：

       这里需要注意的时候，由于GPS距离地面的高度一般为20000km，而这里的同步卫星只有350km，所以看上去会效果不明显，所以这里我们把这里的参数设置的大些，这样看上去效果稍微明显点。然后你再写论文的时候，如果用到其中的数据，只要把他改回350即可。另外，其周期为1.5小时，这样在房子的时候，速度太快，不容易观察，这里稍微设置的大些，使用周期为6小时。