首先：

将GEO卫星，IGSO卫星和MEO卫星进行分类，下列链接提供了相应北斗卫星的PRN号，方便对北斗卫星进行分类。

中国卫星导航系统管理办公室测试评估研究中心

根据其含有的卫星PRN号进行分类处理，其中GEO卫星要进行5°偏差改正。

其中给出了改正矩阵。

其次，给出计算公式如图:

计算公式已经给出，接下来是准备文件。

计算应准备文件：

导航星历文件，这里我截取了只有北斗卫星的卫星导航星历电文，文件格式如下：

其中的END OF HEADER是模拟读取卫星导航星历文件的格式进行编写的。

该代码只进行了GEO和其他卫星的卫星位置计算，还未进行对IGSO、MEO卫星的分类工作，大家可以对代码进行功能添加，可以将数据导入到excel进行筛选，完成对IGSO、MEO卫星的分类。

代码如下：

import math as m
import numpy as np
import csv
with open('C:\\data\\nav\\2022\\145\\daily\\北斗卫星.txt', 'r') as f:if f == 0:print("不能打开文件！")else:print("导航文件打开成功！")nfile_lines = f.readlines()  # 按行读取N文件print(len(nfile_lines))f.close()def start_num():  # 定义数据记录的起始行for i in range(len(nfile_lines)):if nfile_lines[i].find('END OF HEADER') != -1:start_num = i + 1return start_num
def rx(fai):result=np.mat([[1,0,0],[0,m.cos(fai),m.sin(fai)],[0,-1*m.sin(fai),m.cos(fai)]])return result
def rz(fai):result=np.mat([[m.cos(fai),m.sin(fai),0],[-1*m.sin(fai),m.cos(fai),0],[0,0,1]])return resultn_dic_list = []n_data_lines_nums = int((len(nfile_lines) - start_num()) / 8)
print("一共%d组数据" % (n_data_lines_nums))# 第j组，第i行
for j in range(n_data_lines_nums):n_dic = {}for i in range(8):data_content = nfile_lines[start_num() + 8 * j + i]n_dic['数据组数'] = j + 1if i == 0:n_dic['卫星PRN号'] = str(data_content[0:3])n_dic['历元'] = data_content[4:23]n_dic['卫星钟偏差(s)'] = float((data_content.strip('\n')[23:42]))  # 利用字符串切片功能来进行字符串的修改n_dic['卫星钟漂移(s/s)'] = float((data_content.strip('\n')[42:61]))n_dic['卫星钟漂移速度(s/s*s)'] = float((data_content.strip('\n')[61:80]))if i == 1:n_dic['IODE'] = float((data_content.strip('\n')[4:23]))n_dic['C_rs'] = float((data_content.strip('\n')[23:42]))n_dic['n'] = float((data_content.strip('\n')[42:61]))n_dic['M0'] = float((data_content.strip('\n')[61:80]))if i == 2:n_dic['C_uc'] = float((data_content.strip('\n')[4:23]))n_dic['e'] = float((data_content.strip('\n')[23:42]))n_dic['C_us'] = float((data_content.strip('\n')[42:61]))n_dic['sqrt_A'] = float((data_content.strip('\n')[61:80]))if i == 3:n_dic['TEO'] = float((data_content.strip('\n')[4:23]))n_dic['C_ic'] = float((data_content.strip('\n')[23:42]))n_dic['OMEGA'] = float((data_content.strip('\n')[42:61]))n_dic['C_is'] = float((data_content.strip('\n')[61:80]))if i == 4:n_dic['I_0'] = float((data_content.strip('\n')[4:23]))n_dic['C_rc'] = float((data_content.strip('\n')[23:42]))n_dic['w'] = float((data_content.strip('\n')[42:61]))n_dic['OMEGA_DOT'] = float((data_content.strip('\n')[61:80]))if i == 5:n_dic['IDOT'] = float((data_content.strip('\n')[4:23]))n_dic['L2_code'] = float((data_content.strip('\n')[23:42]))n_dic['PS_week_num'] = float((data_content.strip('\n')[42:61]))if i == 6:n_dic['卫星精度(m)'] = float((data_content.strip('\n')[4:23]))n_dic['卫星健康状态'] = float((data_content.strip('\n')[23:42]))n_dic['TGD'] = float((data_content.strip('\n')[42:61]))n_dic['IODC'] = float((data_content.strip('\n')[61:80]))n_dic_list.append(n_dic)
with open('C:\\data\\北斗卫星.csv', 'w', newline='') as f:header = ['数据组数', '卫星PRN号', '历元', '卫星钟偏差(s)', '卫星钟漂移(s/s)', '卫星钟漂移速度(s/s*s)', 'IODE','C_rs', 'n', 'M0', 'C_uc', 'e', 'C_us', 'sqrt_A', 'TEO', 'C_ic', 'OMEGA', 'C_is', 'I_0', 'C_rc', 'w','OMEGA_DOT', 'IDOT', 'L2_code', 'PS_week_num', 'L2_P_code', '卫星精度(m)', '卫星健康状态', 'TGD', 'IODC','X','Y','Z']writer = csv.DictWriter(f, fieldnames=header)writer.writeheader()writer.writerows(n_dic_list)
f.close()
prn_x_y_z=[]
with open('C:\\data\\北斗卫星.csv', 'rt') as csvfile:reader = csv.DictReader(csvfile)for row in reader:PRN = str(row["卫星PRN号"])TIME = row["历元"]year = int(TIME.strip('\n')[2:4])month = int(TIME.strip('\n')[5:7])day = int(TIME.strip('\n')[8:10])hour = int(TIME.strip('\n')[11:13])minute = int(TIME.strip('\n')[14:16])second = float(TIME.strip('\n')[17:19])a_0 = float(row["卫星钟偏差(s)"])a_1 = float(row["卫星钟漂移(s/s)"])a_2 = float(row["卫星钟漂移速度(s/s*s)"])IODE = float(row["IODE"])C_rs = float(row["C_rs"])δn = float(row["n"])M0 = float(row["M0"])C_uc = float(row["C_uc"])e = float(row["e"])C_us = float(row["C_us"])sqrt_A = float(row["sqrt_A"])TEO = float(row["TEO"])C_ic = float(row["C_ic"])OMEGA = float(row["OMEGA"])C_is = float(row["C_is"])I_0 = float(row["I_0"])C_rc = float(row["C_rc"])w = float(row["w"])OMEGA_DOT = float(row["OMEGA_DOT"])IDOT = float(row["IDOT"])L2_code = float(row["L2_code"])PS_week_num = float(row["PS_week_num"])# 卫星精度(m)=2# 卫星健康状态=0TGD = float(row["TGD"])IODC = float(row["IODC"])t1=None# 1.计算卫星运行平均角速度 GM:WGS84下的引力常数 =3.986005e14，a:长半径GM = 398600500000000n_0 = m.sqrt(GM) / m.pow(sqrt_A, 3)n = n_0 + δn# 2.计算归化时间t_k 计算t时刻的卫星位置  UT：世界时 此处以小时为单位UT = hour + (minute / 60.0) + (second / 3600)# GPS时起始时刻1980年1月6日0点   year是两位数 需要转换到四位数if year >= 80:if year == 80 and month == 1 and day < 6:year = year + 2000else:year = year + 1900else:year = year + 2000if month <= 2:year = year - 1month = month + 12  # 1，2月视为前一年13，14月# 需要将当前需计算的时刻先转换到儒略日再转换到GPS时间JD = (365.25 * year) + int(30.6001 * (month + 1)) + day + UT / 24 + 1720981.5WN = int((JD - 2444244.5) / 7)  # WN:GPS_week number 目标时刻的GPS周t_oc = ((JD - 2444244.5) - (7.0 * WN)) * 24 * 3600.0-14# t_GPS:目标时刻的GPS秒 减去14秒为BDT# 对观测时刻t1进行钟差改正,注意：t1应是由接收机接收到的时间t_k = -14# if t1 is None:#     t_k = 0# else:#     δt = a_0 + a_1(t1 - t_oc) + a_2(t1 - t_oc) ^ 2#     t = t1 - δt#     t_k = t - TEO# if t_k > 302400:#     t_k -= 604800# else :#     t_k += 604800# 3.平近点角计算M_k = M_0+n*t_kM_k = M0 + n * t_k  # 实际应该是乘t_k，但是没有接收机的观测时间，所以为了练手设t_k=0# 4.偏近点角计算 E_k  (迭代计算) E_k = M_k + e*sin(E_k)E = 0;E1 = 1;count = 0;while abs(E1 - E) > 1e-10:count = count + 1E1 = EE = M_k + e * m.sin(E)if count > 1e8:print("计算偏近点角时未收敛！")break# 5.计算卫星的真近点角V_k = m.atan((m.sqrt(1 - e * e) * m.sin(E)) / (m.cos(E) - e));# 6.计算升交距角 u_0(φ_k), ω:卫星电文给出的近地点角距u_0 = V_k + w# 7.摄动改正项 δu、δr、δi :升交距角u、卫星矢径r和轨道倾角i的摄动量δu = C_uc * m.cos(2 * u_0) + C_us * m.sin(2 * u_0)δr = C_rc * m.cos(2 * u_0) + C_rs * m.sin(2 * u_0)δi = C_ic * m.cos(2 * u_0) + C_is * m.sin(2 * u_0)# 8.计算经过摄动改正的升交距角u_k、卫星矢径r_k和轨道倾角 i_ku = u_0 + δur = m.pow(sqrt_A, 2) * (1 - e * m.cos(E)) + δri = I_0 + δi + IDOT * (t_k);  # 实际乘t_k=t-t_oe# 9.计算卫星在轨道平面坐标系的坐标,卫星在轨道平面直角坐标系（X轴指向升交点）中的坐标为：x_k = r * m.cos(u)y_k = r * m.sin(u)# 10.观测时刻升交点经度Ω_k的计算，升交点经度Ω_k等于观测时刻升交点赤经Ω与格林尼治恒星时GAST之差  Ω_k=Ω_0+(ω_DOT-omega_e)*t_k-omega_e*t_oeomega_e = 7.292115e-5  # 地球自转角速度OMEGA_k = OMEGA + (OMEGA_DOT - omega_e) * t_k - omega_e * TEO;  # 星历中给出的Omega即为Omega_o=Omega_t_oe-GAST_w# 11.计算卫星在地固系中的直角坐标lX_k = x_k * m.cos(OMEGA_k) - y_k * m.cos(i) * m.sin(OMEGA_k)Y_k = x_k * m.sin(OMEGA_k) + y_k * m.cos(i) * m.cos(OMEGA_k)Z_k = y_k * m.sin(i)# 12.判断卫星是否为GEO卫星，否则不进行极移改正。if PRN in ['C01','C02','C03','C04','C05','C59','C60','C61']:fi=omega_e*t_kfive=180/m.pi*5a=np.mat([X_k,Y_k,Z_k])a=rx(fi)*rz(-1*five)*a.TX_k=str(a[0,0])Y_k=str(a[1,0])Z_k=str(a[2,0])if month > 12:  # 恢复历元year = year + 1month = month - 12print("历元：", year, "年", month, "月", day, "日", hour, "时", minute, "分", second, "秒", "卫星PRN号:",PRN,"平均角速度:", n, "卫星平近点角:", M_k, "偏近点角:", E, "真近点角:", V_k, "升交距角:", u_0,"摄动改正项:", δu,δr, δi, "经摄动改正后的升交距角、卫星矢径和轨道倾角:", u, r, i, "轨道平面坐标X,Y:", x_k, y_k,"观测时刻升交点经度:", OMEGA_k, "地固直角坐标系(极移改正)X:", X_k, "地固直角坐标系Y（极移改正）:", Y_k,"地固直角坐标系Z（极移改正）:",Z_k)prn_x_y_z.append(PRN+",")prn_x_y_z.append(str(X_k)+",")prn_x_y_z.append(str(Y_k)+",")prn_x_y_z.append(str(Z_k))prn_x_y_z.append("\n")else:prn_x_y_z.append(PRN + ",")prn_x_y_z.append(str(X_k) + ",")prn_x_y_z.append(str(Y_k) + ",")prn_x_y_z.append(str(Z_k))prn_x_y_z.append("\n")print("卫星坐标数据计算完成！")
f=open("C:\\data\\北斗卫星位置.txt",'w')
f.writelines(prn_x_y_z)
print("文件写入成功！")
f.close()

在运行结束后，与卫星精密星历进行比对，发现有些卫星的坐标的正负号出现问题，但都是整体性出现的问题，有的卫星坐标XYZ三者都为正确坐标的相反数，不知是何问题，还请广大读者予以指正，该代码仅供参考，还请广大读者进行批评指正。

结果图如下：

学习不易，如有转载，请与博主联系

【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）相关推荐

基于卫星星历计算卫星在CGCS2000大地坐标系中的坐标
目录一.北斗系统概述 1.空间星座 2.坐标系统 3.时间系统二.实验目的三.实验内容四.实验过程五.实验结果一.北斗系统概述 1.空间星座北斗卫星导航系统简称北斗系统,英文缩写为 BD ...
【python】利用广播星历计算BDS卫星的位置
前言本程序为<卫星导航定位基础>大作业之二,功能为实现对广播星历文件的读取和处理,计算出北斗卫星的位置坐标,并绘制出二维和三维的卫星位置分布图.若需要对其他类型卫星数据处理,可根据本程序 ...
利用GPS北斗卫星信号开发设计NTP网络时间服务器
利用GPS北斗卫星信号开发设计NTP网络时间服务器利用GPS北斗卫星信号开发设计NTP网络时间服务器引言准确的时间是天文观测所必需的.天文望远镜在特定时间内的准确指向.CCD曝光时间的控制以及不 ...
时间的基本概念及GPS北斗卫星时钟授时技术
时间的基本概念及GPS北斗卫星时钟授时技术时间的基本概念及GPS北斗卫星时钟授时技术 GPS时间同步的原理和技术 1.有关时间的一些基本概念: 时间与频率之间互为倒数关系,两者密不可分,时间标准的基 ...
揭秘!中国人一定要知道的北斗卫星系统
文章目录简介北斗系统简介北斗一号北斗二号北斗三号北斗系统的三个组成部分北斗系统的三频服务北斗系统提供的服务授时服务怎么通过卫星定位星历表和卫星位置未知的时钟准确的时钟电离层 ...
精准授时，GPS北斗卫星授时同步时钟系统的天花板
精准授时,GPS北斗卫星授时同步时钟系统的天花板精准授时,GPS北斗卫星授时同步时钟系统的天花板对于一个进入信息社会的现代化大国,导航定位和授时系统是最重要的,而且也是最关键的国家基础设施之一.精 ...
【重器】GPS北斗卫星时钟基准与卫星授时服务技术原理
[重器]GPS北斗卫星时钟基准与卫星授时服务技术原理 [重器]GPS北斗卫星时钟基准与卫星授时服务技术原理 1.前言由计算机网络系统组成的分布式系统,若想协调一致进行:IT行业的"整点开拍 ...
赞！北斗卫星助力NTP时钟服务器开启计时服务
赞!北斗卫星助力NTP时钟服务器开启计时服务赞!北斗卫星助力NTP时钟服务器开启计时服务精确时钟自动校准技术,是一种简便的获取北斗卫星精确时间信息的专利技术,具有灵敏度高.不受时间及地域限制等特点 ...
python读取导航电文并计算卫星位置
python简单计算卫星位置前言一.思路总的可分为两个部分:获取导航参数和计算卫星位置. ①获取导航参数:首先讲导航星历中的数据切片,存入csv文件中,再读取csv文件的数据并赋值给各参数 ②计 ...

【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）

首先：

其次，给出计算公式如图:

计算应准备文件：

【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）相关推荐

最新文章

热门文章