利用北斗BDS广播星历，进行卫星位置计算

2024-03-25 11:30
文章标签 进行 广播 计算 位置 卫星 北斗 bds 星历

本文主要是介绍利用北斗BDS广播星历，进行卫星位置计算，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

目录

一、北斗星历数据获取

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

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

二、toc时刻计算

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

3.1 修正Tsv的卫星钟偏

3.2 计算规划时间

3.3 计算平近点角

3.4 迭代求解偏近点角

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

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

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

4.1 MEO/IGSO卫星位置计算

4.2 GEO卫星位置计算

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

六、计算卫星钟偏

七、MATLAB代码&运行示例

一、北斗星历数据获取

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

下面以RINEX格式文件为例：

这里给出的北斗卫星分别是：C05  GEO卫星、C16  MEO/IGSO卫星。这是两种不同轨道的北斗卫星：GEO地球静止轨道卫星、MEO中圆地球轨道卫星、IGSO倾斜地球同步轨道卫星。其GEO卫星位置的计算方法与MEO/IGSO卫星不同，下面过程注意区分。

     3.05           N: GNSS NAV DATA    C: BDS              RINEX VERSION / TYPE
BNC 2.13.0          BKG Frankfurt       20231102 101016 UTC PGM / RUN BY / DATE
Concatenated RINEX files (8)                                COMMENTEND OF HEADER
C05 23 11 02 08 00 00.0 1.385881332680e-04-6.030731469760e-11 0.000000000000e+001.000000000000e+00-2.278125000000e+02-5.663093033410e-09-6.708219322780e-01-7.327646017070e-06 7.320673903450e-04-1.202011480930e-05 6.493444162370e+033.744000000000e+05-5.634501576420e-08-6.255581541050e-01 3.259629011150e-089.409420036110e-02 3.626718750000e+02-1.792511645600e+00 6.869929017260e-092.378670509750e-10                    9.300000000000e+02                   2.000000000000e+00 0.000000000000e+00 2.999999970670e-10-9.100000000000e-093.779100000000e+05 0.000000000000e+00                                      
C16 23 11 02 08 00 00.0-1.551879104230e-04-1.854782993860e-11 0.000000000000e+001.000000000000e+00-3.535937500000e+01 1.561493613910e-09-2.192800555490e+00-1.096632331610e-06 6.275521009230e-03 8.363742381330e-06 6.493531814580e+033.744000000000e+05-2.700835466380e-08 2.343958791130e+00 5.215406417850e-089.618269834520e-01-8.421875000000e+00-2.286522345190e+00-1.925437345050e-09-7.286017777600e-11                    9.300000000000e+02                   2.000000000000e+00 0.000000000000e+00-2.499999984810e-09 4.800000000000e-093.744000000000e+05 0.000000000000e+00

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

1、星历参考时刻 toc：第一行第一个数据   23 11 02 08 00 00.0

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

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

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

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

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

7、平近点角 M0：第二行第四个数据   -6.708219322780e-01

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

9、偏心率 e：第三行第二个数据   7.320673903450e-04

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

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

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

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

14、升交点赤经 omega0：第四行第三个数据   -6.255581541050e-01

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

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

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

18、近地点角距 omega：第五行第三个数据   -1.792511645600e+00

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

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

21、卫星设备时延 TGD1：第七行第三个数据   2.999999970670e-10

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

1、圆周率 pi：3.1415926535898

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

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

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

二、toc时刻计算

从RINEX文件中获取到的toc时刻是年月日时分秒的形式，toc时刻要将其化成周内秒的形式。比如：23 11 02 08 00 00.0代表的是2023年11月2日即周四的 8:00，其周内秒 toc = 4*24*3600+8*3600，toc = 374400s。我们发现 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 轨道参数计算，并考虑摄动改正项

升交点角距改正项：

轨道半径改正项：

轨道倾角改正项：

升交点角距：

轨道半径：

轨道倾角：

升交点赤经：

如果是GEO卫星升交点赤经计算公式如下：

如果是MEO/IGSO卫星升交点赤经计算公式如下：

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

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

4.1 MEO/IGSO卫星位置计算

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

4.2 GEO卫星位置计算

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

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

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

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

六、计算卫星钟偏

七、MATLAB代码&运行示例

主程序 main.m：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算BD卫星的位置和速度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;% 卫星发射时刻
tsv = 3.744000000000e+05 + 20;% MEO卫星的星历参数
MEO_eph.toc = 3.744000000000e+05;
MEO_eph.af2 = 0.000000000000e+00;
MEO_eph.af1 = -1.854782993860e-11;
MEO_eph.af0 = -1.551879104230e-04;
MEO_eph.TGD1 = -2.499999984810e-09;
MEO_eph.sqrtA = 6.493531814580e+03;
MEO_eph.toe = 3.744000000000e+05;
MEO_eph.deltan = 1.561493613910e-09;
MEO_eph.M0 = -2.192800555490e+00;
MEO_eph.e = 6.275521009230e-03;
MEO_eph.i0 = 9.618269834520e-01;
MEO_eph.IDOT = -7.286017777600e-11;
MEO_eph.omega = -2.286522345190e+00;
MEO_eph.Cuc = -1.096632331610e-06;
MEO_eph.Cus = 8.363742381330e-06;
MEO_eph.Crc = -8.421875000000e+00;
MEO_eph.Crs = -3.535937500000e+01;
MEO_eph.Cic = -2.700835466380e-08;
MEO_eph.Cis = 5.215406417850e-08;
MEO_eph.omega0 = 2.343958791130e+00;
MEO_eph.omegaDot = -1.925437345050e-09;% GEO卫星的星历参数
GEO_eph.toc = 3.744000000000e+05;
GEO_eph.af2 = 0.000000000000e+00;
GEO_eph.af1 = -6.030731469760e-11;
GEO_eph.af0 = 1.385881332680e-04;
GEO_eph.TGD1 = 2.999999970670e-10;
GEO_eph.sqrtA = 6.493444162370e+03;
GEO_eph.toe = 3.744000000000e+05;
GEO_eph.deltan = -5.663093033410e-09;
GEO_eph.M0 = -6.708219322780e-01;
GEO_eph.e = 7.320673903450e-04;
GEO_eph.i0 = 9.409420036110e-02;
GEO_eph.IDOT = 2.378670509750e-10;
GEO_eph.omega = -1.792511645600e+00;
GEO_eph.Cuc = -7.327646017070e-06;
GEO_eph.Cus = -1.202011480930e-05;
GEO_eph.Crc = 3.626718750000e+02;
GEO_eph.Crs = -2.278125000000e+02;
GEO_eph.Cic = -5.634501576420e-08;
GEO_eph.Cis = 3.259629011150e-08;
GEO_eph.omega0 = -6.255581541050e-01;
GEO_eph.omegaDot = 6.869929017260e-09;% 分别计算GEO、MEO卫星位置
[MEO_positionX, MEO_positionY, MEO_positionZ, MEO_Vx, MEO_Vy, MEO_Vz, MEO_sat_Clk_error] = MEO_state_tsv(tsv,MEO_eph);
[GEO_positionX, GEO_positionY, GEO_positionZ, GEO_sat_Clk_error] = GEO_state_tsv(tsv,GEO_eph);% 打印输出
fprintf('GEO卫星位置：[%f, %f, %f]\n',GEO_positionX,GEO_positionY,GEO_positionZ);
fprintf('GEO卫星钟偏：%f\n',GEO_sat_Clk_error);
fprintf('\n');
fprintf('MEO卫星位置：[%f, %f, %f]\n',MEO_positionX,MEO_positionY,MEO_positionZ);
fprintf('MEO卫星速度：[%f, %f, %f]\n',MEO_Vx,MEO_Vy,MEO_Vz);
fprintf('MEO卫星钟偏：%f\n',MEO_sat_Clk_error);

函数 MEO_state_tsv.m：

function [positionX,positionY,positionZ,Vx,Vy,Vz,sat_Clk_error] = MEO_state_tsv(tsv,EPH)
%计算MEO卫星的位置%%%%%%%%%%%%%%%%
%数据准备
%%%%%%%%%%%%%%%%
radv = 7.2921151467e-5;
GM = 3.986005e14;
F = -4.442807633e-10;
toc = EPH.toc;
Tsv = tsv;
af0 = EPH.af0;
af1 = EPH.af1;
af2 = EPH.af2;
TGD1 = EPH.TGD1;
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 - TGD1*1e-9;
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 1Ei = M + e*sin(E0);delt_E = Ei - E0;delt_E = rem(delt_E,2*pi);if abs(delt_E) < 1e-12E = Ei;break;elseE0 = Ei;end
end%求解真近点角和升交点角距
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);%更新卫星钟偏
dtr = F*e*sqrtA*sin(E);
sat_Clk_error = af2*delt^2 + af1*delt + af0 - TGD1*1e-9 + dtr;end

函数 GEO_state_tsv.m：

function [positionX,positionY,positionZ,sat_Clk_error] = GEO_state_tsv(tsv,EPH)
%计算GEO卫星的位置%%%%%%%%%%%%%%%%
%数据准备
%%%%%%%%%%%%%%%%
radv = 7.2921151467e-5;
GM = 3.986005e14;
F = -4.442807633e-10;
toc = EPH.toc;
Tsv = tsv;
af0 = EPH.af0;
af1 = EPH.af1;
af2 = EPH.af2;
TGD1 = EPH.TGD1;
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 - TGD1*1e-9;
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 1Ei = M + e*sin(E0);delt_E = Ei - E0;delt_E = rem(delt_E,2*pi);if abs(delt_E) < 1e-12E = Ei;break;elseE0 = Ei;end
end%求解真近点角和升交点角距
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*tk-radv*toe;
OMEGA = mod(OMEGA,2*pi);%卫星坐标WGS-84
X = cos(u)*r*cos(OMEGA) - sin(u)*r*cos(i)*sin(OMEGA);
Y = cos(u)*r*sin(OMEGA) + sin(u)*r*cos(i)*cos(OMEGA);
Z = sin(u)*r*sin(i);%GEO卫星的位置
f = -5/180*pi;
p = radv*tk;
Rx = [1 0 0; 0 cos(f) sin(f); 0 -sin(f) cos(f)];
Rz = [cos(p) sin(p) 0; -sin(p) cos(p) 0; 0 0 1];
position = Rz*Rx*[X;Y;Z];
positionX = position(1);
positionY = position(2);
positionZ = position(3);%更新卫星钟偏
dtr = F*e*sqrtA*sin(E);
sat_Clk_error = af2*delt^2 + af1*delt + af0 - TGD1*1e-9 + dtr;end

运行结果：

这篇关于利用北斗BDS广播星历，进行卫星位置计算的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/844869。

相关文章

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯，先赞后看哦~🎈🎈 🏆 作者简介：景天科技苑 🏆《头衔》：大厂架构师，华为云开发者社区专家博主，阿里云开发者社区专家博主，CSDN全栈领域优质创作者，掘金优秀博主，51CTO博客专家等。 🏆《博客》：Python全栈，前后端开发，小程序开发，人工智能，js逆向，App逆向，网络系统安全，数据分析，Django，fastapi
阅读更多...
业务中14个需要进行A/B测试的时刻[信息图]

业务中14个需要进行A/B测试的时刻[信息图]

在本指南中，我们将全面了解有关 A/B测试 的所有内容。 我们将介绍不同类型的A/B测试，如何有效地规划和启动测试，如何评估测试是否成功，您应该关注哪些指标，多年来我们发现的常见错误等等。 什么是A/B测试？ A/B测试（有时称为“分割测试”）是一种实验类型，其中您创建两种或多种内容变体——如登录页面、电子邮件或广告——并将它们显示给不同的受众群体，以查看哪一种效果最好。 本质上，A/B测
阅读更多...
poj 1113 凸包+简单几何计算

poj 1113 凸包+简单几何计算

题意： 给N个平面上的点，现在要在离点外L米处建城墙，使得城墙把所有点都包含进去且城墙的长度最短。 解析： 韬哥出的某次训练赛上A出的第一道计算几何，算是大水题吧。 用convexhull算法把凸包求出来，然后加加减减就A了。 计算见下图： 好久没玩画图了啊好开心。 代码： #include <iostream>#include <cstdio>#inclu
阅读更多...
uva 1342 欧拉定理（计算几何模板）

uva 1342 欧拉定理（计算几何模板）

题意： 给几个点，把这几个点用直线连起来，求这些直线把平面分成了几个。 解析： 欧拉定理： 顶点数 + 面数 - 边数= 2。 代码： #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#inc
阅读更多...
uva 11178 计算集合模板题

uva 11178 计算集合模板题

题意： 求三角形行三个角三等分点射线交出的内三角形坐标。 代码： #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vector>#include <
阅读更多...
XTU 1237 计算几何

XTU 1237 计算几何

题面： Magic Triangle Problem Description: Huangriq is a respectful acmer in ACM team of XTU because he brought the best place in regional contest in history of XTU. Huangriq works in a big compa
阅读更多...
POJ1269 判断2条直线的位置关系

POJ1269 判断2条直线的位置关系

题目大意：给两个点能够确定一条直线，题目给出两条直线（由4个点确定），要求判断出这两条直线的关系：平行，同线，相交。如果相交还要求出交点坐标。 解题思路： 先判断两条直线p1p2, q1q2是否共线， 如果不是，再判断 直线 是否平行， 如果还不是， 则两直线相交。  判断共线：  p1p2q1 共线 且 p1p2q2 共线 ，共线用叉乘为 0  来判断，  判断 平行：  p1p
阅读更多...
音视频入门基础：WAV专题（10）——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

音视频入门基础：WAV专题（10）——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

一、引言 从文章《音视频入门基础：WAV专题（6）——通过FFprobe显示WAV音频文件每个数据包的信息》中我们可以知道，通过FFprobe命令可以打印WAV音频文件每个packet（也称为数据包或多媒体包）的信息，这些信息包含该packet的pts、dts： 打印出来的“pts”实际是AVPacket结构体中的成员变量pts，是以AVStream->time_base为单位的显
阅读更多...
遮罩，在指定元素上进行遮罩

遮罩，在指定元素上进行遮罩

废话不多说，直接上代码： ps:依赖 jquer.js 1.首先，定义一个 Overlay.js  代码如下： /*遮罩 Overlay js 对象*/function Overlay(options){//{targetId:'',viewHtml:'',viewWidth:'',viewHeight:''}try{this.state=false;//遮罩状态 true 激活，f
阅读更多...
利用matlab bar函数绘制较为复杂的柱状图，并在图中进行适当标注

利用matlab bar函数绘制较为复杂的柱状图，并在图中进行适当标注

示例代码和结果如下：小疑问：如何自动选择合适的坐标位置对柱状图的数值大小进行标注？😂 clear; close all;x = 1:3;aa=[28.6321521955954 26.2453660695847 21.69102348512086.93747104431360 6.25442246899816 3.342835958564245.51365061796319 4.87
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2