利用GPS广播星历，进行卫星位置计算

KD学长

已于 2023-12-14 16:16:42 修改

阅读量8.2k

点赞数 45

分类专栏： GNSS 文章标签：算法笔记经验分享

于 2023-12-14 16:11:18 首次发布

本文链接：https://blog.csdn.net/qq_45568498/article/details/134950531

版权

GNSS 专栏收录该内容

10 篇文章

订阅专栏

三、根据卫星发射时刻Tsv，计算卫星轨道参数

一、GPS星历数据获取

GPS广播星历获取途径：RINEX格式文件、导航电文解码。

下面以RINEX格式文件为例：

     3.03           NAVIGATION DATA     M (Mixed)           RINEX VERSION / TYPE
BCEmerge            congo               20190315 004602 GMT PGM / RUN BY / DATE 
Merged GPS/GLO/GAL/BDS/QZS/SBAS/IRNSS navigation file       COMMENT 
based on CONGO and MGEX tracking data                       COMMENT 
DLR: O. Montenbruck;  TUM: P. Steigenberger                 COMMENT 
GAUT  1.8626451492e-09-8.881784197e-16 259200 2044          TIME SYSTEM CORR    
GPGA  9.3132257462e-10 1.820765760e-14 345600 2044          TIME SYSTEM CORR    
GPUT  9.3132257462e-10 0.000000000e+00 589824 2044          TIME SYSTEM CORR    
QZUT  9.3132257462e-10 0.000000000e+00 528384 2044          TIME SYSTEM CORR    
    18                                                      LEAP SECONDS        
                                                            END OF HEADER       
G01 2019 03 14 02 00 00-1.803217455745e-04-7.503331289627e-12 0.000000000000e+00
     5.600000000000e+01-4.325000000000e+01 4.276249551530e-09 1.676641337630e+00
    -2.229586243629e-06 8.502699085511e-03 9.126961231232e-06 5.153659761429e+03
     3.528000000000e+05-2.011656761169e-07-1.602363456915e+00 8.009374141693e-08
     9.744228427026e-01 2.130625000000e+02 6.832534868102e-01-7.954974213749e-09
     1.032185851827e-10 1.000000000000e+00 2.044000000000e+03 0.000000000000e+00
     2.000000000000e+00 0.000000000000e+00 5.587935447693e-09 5.600000000000e+01
     3.456000000000e+05 4.000000000000e+00

1.1 计算卫星位置需要用到的星历数据：

1、星历参考时刻 toc：第一行第一个数据 2019 03 14 02 00 00

2、卫星钟差 af0：第一行第二个数据 -1.803217455745e-04

3、卫星钟偏 af1：第一行第三个数据 -7.503331289627e-12

4、卫星钟偏移 af2：第一行第四个数据 0.000000000000e+00

5、轨道半径改正项 Crs：第二行第二个数据 -4.325000000000e+01

6、角速度改正项 deltan：第二行第三个数据 4.276249551530e-09

7、平近点角 M0：第二行第四个数据 1.676641337630e+00

8、升交点角距改正 Cuc：第三行第一个数据 -2.229586243629e-06

9、偏心率 e：第三行第二个数据 8.502699085511e-03

10、升交点角距改正 Cus：第三行第三个数据 9.126961231232e-06

11、轨道长半轴 sqrtA：第三行第四个数据 5.153659761429e+03

12、周内秒 toe：第四行第一个数据 3.528000000000e+05

13、轨道倾角改正项 Cic：第四行第二个数据 -2.011656761169e-07

14、升交点赤经 omega0：第四行第三个数据 -1.602363456915e+00

15、轨道倾角改正项 Cis：第四行第四个数据 8.009374141693e-08

16、轨道倾角 i0：第五行第一个数据 9.744228427026e-01

17、轨道半径改正项 Crc：第五行第二个数据 2.130625000000e+02

18、近地点角距 omega：第五行第三个数据 6.832534868102e-01

19、赤经变化 omegaDot：第五行第四个数据 -7.954974213749e-09

20、轨道倾角变率 IDOT：第六行第一个数据 1.032185851827e-10

21、卫星设备时延 TGD：第七行第三个数据 5.587935447693e-09

1.2 计算卫星位置需要用到的常数：

1、圆周率 pi：3.1415926535898

2、地球自转角速度 radv：7.2921151467e-5

3、地球标准引力常数 GM：3.986005e14

4、地球引力常数 F：-4.442807633e-10

二、toc时刻计算

从RINEX文件中获取到的toc时刻是年月日时分秒的形式，toc时刻要将其化成周内秒的形式。比如：2019 03 14 02 00 00 代表的是2019年3月14日即周四的 2:00，其周内秒 toc = 4*24*3600+2*3600，toc = 352800s。我们发现 toc = toe，toc时刻与toe时刻是同步的。

三、根据卫星发射时刻Tsv，计算卫星轨道参数

3.1 修正Tsv的卫星钟偏

卫星的真实发射时刻Tsv由初始Tsv修正卫星钟偏和卫星设备时延后得到，并且这里的Tsv也是周内秒的形式：

其中sat_Clk_error表示卫星钟偏。

3.2 计算规划时间

星历中的toe时刻，表示星历中的所有卫星数据是toe时刻瞬间的。发射时间Tsv的卫星数据，需要通过 Tsv-toe 的差值进行计算，这个差就叫规划时间tk。

3.3 计算平近点角

轨道长半轴：

平均角速度：

平近点角：

3.4 迭代求解偏近点角

偏近点角的计算公式：

为了偏近点角的准确性，采用迭代的过程求解。

3.5 求解真近点角和升交点角距

真近点角：

升交点角距：

3.6 轨道参数计算，并考虑摄动改正项

升交点角距改正项：

轨道半径改正项：

轨道倾角改正项：

升交点角距：

轨道半径：

轨道倾角：

升交点赤经：

四、计算Tsv时刻的卫星位置

根据上面的计算，已经得到了Tsv时刻的卫星轨道参数：升交点角距、轨道半径、轨道倾角、升交点赤经。

卫星位置的ECEF坐标，可由下面公式计算得到：

五、计算Tsv时刻的卫星速度

根据上面的计算，已经得到了Tsv时刻的卫星轨道速度参数：平均角速度 n、偏近点角 E、偏心率 e、升交点角距 phi 等。

卫星 (X,Y,Z) 三个方向上的速度，可根据下面公式得到：

六、计算卫星钟偏

七、MATLAB代码&运行示例

主程序 main.m：

clear;
clc;

% 卫星发射时刻
tsv = 3.528000000000e+05 + 20;

% GPS星历参数
eph.toc = 3.528000000000e+05;
eph.af2 = 0.000000000000e+00;
eph.af1 = -7.503331289627e-12;
eph.af0 = -1.803217455745e-04;
eph.TGD = 5.587935447693e-09;
eph.sqrtA = 5.153659761429e+03;
eph.toe = 3.528000000000e+05;
eph.deltan = 4.276249551530e-09;
eph.M0 = 1.676641337630e+00;
eph.e = 8.502699085511e-03;
eph.i0 = 9.744228427026e-01;
eph.IDOT = 1.032185851827e-10;
eph.omega = 6.832534868102e-01;
eph.Cuc = -2.229586243629e-06;
eph.Cus = 9.126961231232e-06;
eph.Crc = 2.130625000000e+02;
eph.Crs = -4.325000000000e+01;
eph.Cic = -2.011656761169e-07;
eph.Cis = 8.009374141693e-08;
eph.omega0 = -1.602363456915e+00;
eph.omegaDot = -7.954974213749e-09;

% 计算GPS卫星位置、速度、钟差
[positionX, positionY, positionZ, Vx, Vy, Vz, sat_Clk_error] = sv_state_tsv(tsv,eph);


% 命令行打印结果
fprintf('卫星位置：[%f, %f, %f]\n',positionX,positionY,positionZ);
fprintf('卫星速度：[%f, %f, %f]\n',Vx,Vy,Vz);
fprintf('卫星钟偏：%f\n',sat_Clk_error);

函数 sv_state_tsv.m：

function [positionX, positionY, positionZ, Vx, Vy, Vz, sat_Clk_error] = sv_state_tsv(tsv,eph)
% %输入：卫星星历
% %输出：卫星的位置WGS-84、速度、卫星钟偏

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%常量和星历数据准备
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
radv = 7.2921151467e-5;
GM = 3.986005e14;
F = -4.442807633e-10;
toc = eph.toc;
Tsv = tsv;
af2 = eph.af2;
af1 = eph.af1;
af0 = eph.af0;
TGD = eph.TGD;
sqrtA = eph.sqrtA;
toe = eph.toe;
deltan = eph.deltan;
M0 = eph.M0;
e = eph.e;
i0 = eph.i0;
IDOT = eph.IDOT;
omega = eph.omega;
Cuc = eph.Cuc;
Cus = eph.Cus;
Crc = eph.Crc;
Crs = eph.Crs;
Cic = eph.Cic;
Cis = eph.Cis;
omega0 = eph.omega0;
omegaDot = eph.omegaDot;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算卫星位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算规划时间
delt = Tsv - toc;
sat_Clk_error = af2*delt^2 + af1*delt + af0 - TGD;
Tsv = Tsv - sat_Clk_error;                              %信号发射时刻
tk = Tsv - toe;                                         %规划时间

%计算平近点角
a = sqrtA^2;
n = sqrt(GM/(a^3)) + deltan;
M = M0 + n*tk;
M = mod(M,2*pi);

%计算偏近点角
E0 = M;
while 1
    E = M+e*sin(E0);
    delt_E = E - E0;
    delt_E = rem(delt_E,2*pi);
    if abs(delt_E) < 1e-12
        break;
    else
        E0 = E;
    end
end
E = mod(E,2*pi);

%计算真近点角和升交点角距
nu = atan2(sqrt(1-e^2)*sin(E),cos(E)-e);
phi = nu + omega;
phi = rem(phi,2*pi);

%计算摄动修正后的升交点角距、轨道半径、轨道倾角
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 = phi + delt_u;
r = a*(1-e*cos(E))+delt_r;
i = i0 + IDOT*tk + delt_i;

%计算升交点赤经
OMEGA = omega0 + (omegaDot-radv)*tk - radv*toe;
OMEGA = mod(OMEGA,2*pi);

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算卫星速度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E_dot = n/(1-e*cos(E));
Vk_dot = (sqrt(1-e^2)*E_dot)/(1-e*cos(E));

Uk_dot = 2*Vk_dot*(Cus*cos(2*phi)-Cuc*sin(2*phi));
Rk_dot = 2*Vk_dot*(Crs*cos(2*phi)-Crc*sin(2*phi));
Ik_dot = 2*Vk_dot*(Cis*cos(2*phi)-Cic*sin(2*phi));

U_dot = Vk_dot + Uk_dot;
R_dot = a*e*E_dot*sin(E) + Rk_dot;
I_dot = IDOT + Ik_dot;

OMEGA_dot = omegaDot - radv;

x = R_dot*cos(u) - r*U_dot*sin(u);
y = R_dot*sin(u) + r*U_dot*cos(u);

Vx = -positionY*OMEGA_dot - (y*cos(i)-positionZ*I_dot)*sin(OMEGA) + x*cos(OMEGA);
Vy = positionX*OMEGA_dot + (y*cos(i)-positionZ*I_dot)*cos(OMEGA) + x*sin(OMEGA);
Vz = y*sin(i) + y*I_dot*cos(i); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算卫星钟偏
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delt = Tsv - toc;
dtr = F*e*sqrtA*sin(E);
sat_Clk_error = af2*delt^2 + af1*delt + af0 - TGD + dtr;


end

运行结果：