通过MATLAB模拟24个GPS卫星的轨道运行效果
目录
一、理论基础
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小时。
