【python】利用广播星历计算BDS卫星的位置

前言

本程序为《卫星导航定位基础》大作业之二，功能为实现对广播星历文件的读取和处理，计算出北斗卫星的位置坐标，并绘制出二维和三维的卫星位置分布图。若需要对其他类型卫星数据处理，可根据本程序修改增进。

本文章部分代码借鉴于@学测绘的小杨【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）

获取广播星历文件

可以通过下列链接进行下载

1.ftp://igs.gnsswhu.cn/pub/gnss/mgex/daily/rinex3/2020/

2.ftp://epncb.oma.be/pub/obs/BRDC/2023/

3.ftp://pub:tarc@ftp2.csno-tarc.cn/almanac/2023

注意：如果浏览器不能访问上面的FTP地址，则需用FTP软件进行下载

这里笔者提供另一个方法下载，以第一个链接为例。

1.先复制链接

2.打开“我的电脑”

3.将链接粘贴在路径上运行

4.任选一个文件夹打开

5.将压缩包复制粘贴到存放数据处解压即可

解读广播星历文件的格式

每一行参数代表意义如下：

卫星PRN号卫星钟时间卫星时钟偏差(s) 卫星时钟漂移(s/s) 卫星时钟漂移率(s/s²)

数据、星历发布时间(数据期龄) 轨道半径的正弦调和改正项的振幅(m) 卫星平均运动速率与计算值之差(rad/s) 参考时间的平近点角(rad)

维度幅角的余弦调和改正项的振幅(rad) 轨道偏心率轨道幅角的正弦调和改正项的振幅(rad) 长半轴平方根

星历的参考时刻轨道倾角的余弦调和改正项的振幅(rad) 参考时刻的升交点赤经维度倾角的正弦调和改正项的振幅(rad）

参考时间的轨道倾角(rad) 轨道平径的余弦调和改正项的振幅(m) 近地点角距升交点赤经变化率(rad)

轨道倾角变化率(rad/s) L2频道C/A码标识 GPS周数（toe）L2 P数据标志

卫星精度卫星健康状态 TGD电频机延迟改正数 IODC时钟数据有效期

电文发送时间拟合区间（h）

从广播星历文件中读取所需的数据并进行计算

计算步骤与算法如下：

来源于《GPS测量原理及应用第四版》（武汉大学出版社）

需要注意的是：

计算BDS卫星时，取地心引力常数μ=GM=3.986004418e14，地球自转角速度w=7.292115e-5

代码如下：

import csv
import math as m
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as pltwith open('C:\\Users\\20615\\Desktop\\Statellite\\task\\two\\brdm1260.20p', 'r') as f:if f == 0:print("不能打开文件！")else:print("星历文件已打开！")nfile_lines = f.readlines()                                       #按行读取N文件print(f'文件一共{len(nfile_lines)}行')def Rx(fai):rx = np.mat([[1, 0, 0], [0, m.cos(fai), m.sin(fai)], [0, -m.sin(fai), m.cos(fai)]])return rxdef Rz(fai):rz = np.mat([[m.cos(fai), m.sin(fai), 0], [-m.sin(fai), m.cos(fai), 0], [0, 0, 1]])return rzfor i in range(len(nfile_lines)):                                     #获取开始行if nfile_lines[i].find('C01') != -1:start_num = i + 1breakfor i in range(len(nfile_lines)):                                     #获取终止行if nfile_lines[i].find('J01') != -1:end_num = ibreakx_list = []     #二维X坐标
y_list = []     #二维Y坐标
X_list = []     #三维X坐标
Y_list = []     #三维Y坐标
Z_list = []     #三维Z坐标
satellite_lines = int(((end_num - start_num) + 1) / 8)                #得到总数据数
print(f'一共{satellite_lines}组数据')#第j组 第i行
for j in range(satellite_lines):for i in range(8):                                                  #使用.strip('\n)以去除空格，使数据能够从str转为floatdata_content = nfile_lines[start_num + 8 * j + i - 1]           #如果直接跳过空格读取，会导致负号丢失！！！if i == 0:                                                      PRN = data_content[0:3]                                     #卫星号                                  year = int(data_content[4:8])                               #年month = int(data_content[9:11])                             #月day = int(data_content[12:14])                              #日hour = int(data_content[15:17])                             #时minute = int(data_content[18:20])                           #分second = int(data_content[21:23])                           #秒a0 = float((data_content.strip('\n')[23:42]))               #卫星时钟偏差（s）a1 = float((data_content.strip('\n')[42:61]))               #卫星时钟漂移（s/s）a2 = float((data_content.strip('\n')[61:80]))               #卫星时钟漂移率（s/s²）if i == 1:IODE = float((data_content.strip('\n')[4:23]))              #数据、星历发布时间(数据期龄)Crs = float((data_content.strip('\n')[23:42]))              #轨道半径的正弦调和改正项的振幅（m）delta_n = float((data_content.strip('\n')[42:61]))          #卫星平均运动速率与计算值之差(rad/s)M0 = float((data_content.strip('\n')[61:80]))               #参考时间的平近点角(rad)if i == 2:  Cuc = float((data_content.strip('\n')[4:23]))               #维度幅角的余弦调和改正项的振幅(rad)e = float((data_content.strip('\n')[23:42]))                #轨道偏心率Cus = float((data_content.strip('\n')[42:61]))              #轨道幅角的正弦调和改正项的振幅(rad)sqrtA = float((data_content.strip('\n')[61:80]))            #长半轴平方根if i == 3:toe = float((data_content.strip('\n')[4:23]))               #星历的参考时刻Cic = float((data_content.strip('\n')[23:42]))              #轨道倾角的余弦调和改正项的振幅(rad)OMEGA = float((data_content.strip('\n')[42:61]))            #参考时刻的升交点赤经Cis = float((data_content.strip('\n')[61:80]))              #维度倾角的正弦调和改正项的振幅(rad）if i == 4:i0 = float((data_content.strip('\n')[4:23]))                #参考时间的轨道倾角(rad)Crc = float((data_content.strip('\n')[23:42]))              #轨道平径的余弦调和改正项的振幅(m)omega = float((data_content.strip('\n')[42:61]))            #近地点角距delta_OMEGA = float((data_content.strip('\n')[61:80]))      #升交点赤经变化率(rad)if i == 5:IDOT = float((data_content.strip('\n')[4:23]))              #轨道倾角变化率(rad/s)L2code = float((data_content.strip('\n')[23:42]))           #L2频道C/A码标识week = float((data_content.strip('\n')[42:61]))             #GPS周数（toe）L2Pflag = float((data_content.strip('\n')[61:80]))          #L2 P数据标志if i == 6:sacc = float((data_content.strip('\n')[4:23]))              #卫星精度sHEA = float((data_content.strip('\n')[23:42]))             #卫星健康状态TGD= float((data_content.strip('\n')[42:61]))               #TGD电频机延迟改正数IODC = float((data_content.strip('\n')[61:80]))             #IODC时钟数据有效期#1.计算北斗卫星在轨道平面直角坐标系下的坐标#1.1 计算卫星运行的平均角速度nGM = 398600441800000n0 = m.sqrt(GM) / m.pow(sqrtA, 3)n = n0 + delta_n#1.2 计算归化时间tkif month <= 2:year -= 1month += 12JD = 365.25 * year + int(30.6001 * (month + 1)) + day + 1720981.5 + hour / 24.0 + minute / 1440 + second / 86400  #计算儒略日WN = int((JD - 2444244.5) / 7)                                   # WN:GPS_week number 目标时刻的GPS周toc = ((JD - 2444244.5) - (7.0 * WN)) * 24 * 3600.0-14           # tGPS:目标时刻的GPS秒 减去14秒为BDTtp = 24 * 3600 * day + 3600 * hour + 60 * minute + second        #观测时刻delta_t = a0 + a1 * (tp - toc) + a2 * m.pow(tp - toc, 2)t = tp - delta_t                                                 #观测时刻作卫星钟差改正tk = t - toeif tk > 302400:tk -= 604800if tk < -302400:tk += 604800  #1.3 观测时刻卫星平近点角Mk的计算Mk = M0 + n * tk#1.4 计算偏近点角count = 0E0 = M0Ek = Mk + e * m.sin(E0)while abs(Ek - E0) > 1e-10:count += 1E0 = EkEk = Mk + e * m.sin(E0)if count > 1e8:print("计算偏近点角时未收敛")break#1.5 真近点角Vk的计算Vk = m.atan2((m.sqrt(1 - e * e) * m.sin(Ek)), (m.cos(Ek)) - e)#1.6 升交距角Φk的计算Fai_k = Vk + omega#1.7 摄动改正项δu、δr、δi的计算sigema_u = Cuc * m.cos(2 * Fai_k) + Cus * m.sin(2 * Fai_k)sigema_r = Crc * m.cos(2 * Fai_k) + Crs * m.sin(2 * Fai_k)sigema_i = Cic * m.cos(2 * Fai_k) + Cis * m.sin(2 * Fai_k)#1.8 计算经过摄动改正的升交距角uk、卫星矢径rk和轨道倾角ikuk = Fai_k + sigema_urk = m.pow(sqrtA, 2) * (1 - e * m.cos(Ek)) + sigema_rik = i0 + sigema_i + IDOT * tk#1.9 计算卫星在轨道平面坐标系的坐标xk = rk * m.cos(uk)yk = rk * m.sin(uk)#将卫星平面坐标写入以提供二维绘图数据x_list.append(xk)y_list.append(yk)#2.计算卫星在CGCS2000地固坐标系中的空间直角坐标omega_e = 7.2921150e-5if PRN in ['C01','C02','C03','C04','C05','C59','C60','C61']:        #判断是否为GEO卫星#2.1.1 计算观测时刻升交点精度Ωk(惯性系)OMEGA_k = OMEGA + delta_OMEGA * tk - omega_e * toe              #2.2.1 计算GEO卫星在自定义坐标系中的空间直角坐标XGk = xk * m.cos(OMEGA_k) - yk * m.cos(ik) * m.sin(OMEGA_k)YGk = xk * m.sin(OMEGA_k) + yk * m.cos(ik) * m.cos(OMEGA_k)ZGk = yk * m.sin(ik)#2.3 计算GEO卫星在CGCS2000地固坐标系中的空间直角坐标fi = omega_e * tkfive = -5 * m.pi / 180temp = np.mat([[XGk], [YGk], [ZGk]])temp = Rz(fi) * Rx(five) * tempXk = temp[0,0]Yk = temp[1,0]Zk = temp[2,0]else:#2.1.2 计算观测时刻升交点精度Ωk(地固坐标系)OMEGA_k = OMEGA + (delta_OMEGA - omega_e) * tk - omega_e * toe  #2.2.2 计算MEO和IGSO卫星在CGCS2000地固坐标系中的空间直角坐标Xk = xk * m.cos(OMEGA_k) - yk * m.cos(ik) * m.sin(OMEGA_k)Yk = xk * m.sin(OMEGA_k) + yk * m.cos(ik) * m.cos(OMEGA_k)Zk = yk * m.sin(ik)#将卫星空间坐标写入以提供三维绘图数据X_list.append(Xk)Y_list.append(Yk)Z_list.append(Zk)#计算完毕print('PRN号：%s  X坐标：%-18fY坐标：%-18fZ坐标：%-18f'%(PRN, Xk, Yk, Zk))
print("卫星坐标数据计算完毕！")#3.根据卫星坐标绘制二维与三维图像#3.1 绘制二维图像
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.scatter(x_list, y_list, marker = '.', color = 'red', s = 8)#3.2 绘制三维图像
fig = plt.figure()
ax2 = fig.add_subplot(111, projection='3d')
ax2.scatter(X_list, Y_list, Z_list, color = 'blue', s = 8, alpha = 0.5)
#展示图像
print("已展示卫星位置二维与三维图像")
plt.show()

结果

根据卫星轨道平面坐标绘制的二维图像如下：

根据卫星地固坐标绘制的三维图像如下：

中间的球壳为MEO卫星，旁边8字形半球壳为IGSO卫星，密集分布的点状为GEO卫星