前言
- 本代码为GNSS课程设计代码，仅供参考，使用的计算方法与公式均来源于王坚主编的《卫星定位原理与应用（第二版）》。
- 本代码计算结果可以通过下载精密星历进行比照，误差在1-10m左右。
- 实现功能：读取卫星广播星历，并将其计算为WGS-84坐标系下的坐标，每颗卫星，每15分钟输出一次。
广播星历下载
- 有多重方法进行下载，由于网络原因以及使用的便捷程度，建议使用武汉大学的IGS网站进行下载。（http://www.igs.gnsswhu.cn/index.php）
- 文件名不需要填写，在选择时间时注意，即使只选择一天的数据，设置结束时间也要到第二天，否则会显示错误，检索结果中的数字代表该日期是所选年份的第几天。
Python函数库
本代码使用numpy，pandas，gnss_lib_py,matplotlib四个函数库，请提前安装。
安装代码
#两条命令根据使用环境进行选择 pip install gnss-lib-py pandas matplotlib #Python环境安装代码 conda install gnss-lib-py pandas matplotlib -c conda-forge #conda环境安装代码 gnss_lib_py库
GitHub主页：https://github.com/Stanford-NavLab/gnss_lib_py?tab=readme-ov-file
文档主页：https://gnss-lib-py.readthedocs.io/en/latest/index.html
本文主要使用该库的读取以及转化为DataFrame功能，其中参数的命名规则以及时间转换规则可以在文档中找到。
具体代码
建议使用jupyter进行执行，下方代码为分单元块格式，如果没有jupyter环境可以直接粘贴到一个python文件进行运行。
代码块一：
import gnss_lib_py as glp import datetime import pandas as pd import numpy as np import matplotlib.pyplot as plt 代码块二：
# 导入23n文件 file_path = 'brdc2550.23n' data = glp.RinexNav(file_path) data_df = data.pandas_df() 代码块三：
# 寻找最小差值的参考时刻 def find(inweekmilli, refers, insv): filter_sv = refers[refers['gnss_sv_id'] == insv] defference = np.abs(inweekmilli - filter_sv['t_oe']) return defference.idxmin() times = np.array([None] * 24 * 4) gpsmillis = np.array([None] * 24 * 4) n = 0 for hour in range(0, 24): minut = 0 while minut < 60: times[n] = datetime.datetime(2023, 9, 12, hour, minut, 0, tzinfo=datetime.timezone.utc) gpsmillis[n] = glp.datetime_to_gps_millis(times[n]) minut += 15 n += 1 gpsmillis = np.array(gpsmillis) GM=3.986005E+14 sqrtGM = np.sqrt(GM) sv_list = [f'G{str(i).zfill(2)}' for i in range(1, 33)] outdata = pd.DataFrame(columns=['data', 'gnss_sv_id', 'X', 'Y', 'Z'], index=range(24 * 4 * 32)) orbit = pd.DataFrame(columns=['data', 'gnss_sv_id', 'x', 'y'], index=range(24 * 4 * 32)) m = 0 j = 0 print("正在计算，请稍候。") for gpsmilli in gpsmillis: for sv in sv_list: week, milli_week = glp.gps_millis_to_tow(gpsmilli) milli_week=milli_week-18 index = find(milli_week, data_df, sv) print(f'sv:{sv},time:{times[j]},index:{index}') print(milli_week,data_df.iloc[index]['t_oe']) a = np.power(data_df.iloc[index]['sqrtA'], 2) n0 = sqrtGM / np.power(a, 3 / 2) n = n0 + data_df.iloc[index]['deltaN'] tk = milli_week - data_df.iloc[index]['t_oe'] M = data_df.iloc[index]['M_0'] + n * tk e = data_df.iloc[index]['e'] # 打印中间结果M和e print(f"M: {M}, e: {e}") # 解开普勒方程 E = M for _ in range(50): # 使用迭代方法求解E E = M + e * np.sin(E) # 打印中间结果E print(f"E: {E}") f = np.arctan((np.sqrt(1 - e**2) * np.sin(E)) / (np.cos(E) - e)) if E > np.pi*0.5: f=f+np.pi if E < -np.pi*0.5: f=f-np.pi if np.pi*0.5 > E > 0 > f: f=f+np.pi if -np.pi*0.5 < E < 0 < f: f=f-np.pi print(f"arctan({(np.sqrt(1 - e**2) * np.sin(E)) / (np.cos(E) - e)}),f:{f}") u_pie = data_df.iloc[index]['omega'] + f r_pie = a * (1 - e * np.cos(E)) C_uc = data_df.iloc[index]['C_uc'] C_us = data_df.iloc[index]['C_us'] C_rc = data_df.iloc[index]['C_rc'] C_rs = data_df.iloc[index]['C_rs'] C_ic = data_df.iloc[index]['C_ic'] C_is = data_df.iloc[index]['C_is'] delta_u = C_uc * np.cos(2 * u_pie) + C_us * np.sin(2 * u_pie) delta_r = C_rc * np.cos(2 * u_pie) + C_rs * np.sin(2 * u_pie) delta_i = C_ic * np.cos(2 * u_pie) + C_is * np.sin(2 * u_pie) u = u_pie + delta_u r = r_pie + delta_r i = data_df.iloc[index]['i_0'] + delta_i + data_df.iloc[index]['IDOT'] * tk print(f'u:{u}') x = r * np.cos(u) y = r * np.sin(u) w_e = 7.292115E-5 L = data_df.iloc[index]['Omega_0'] + (data_df.iloc[index]['OmegaDot']- w_e )* milli_week - data_df.iloc[index]['OmegaDot']*data_df.iloc[index]['t_oe'] X = x * np.cos(L) - y * np.cos(i) * np.sin(L) Y = x * np.sin(L) + y * np.cos(i) * np.cos(L) Z = y * np.sin(i) orbit.iloc[m,:] = [times[j],sv,x,y] outdata.iloc[m, :] = [times[j], sv, X, Y, Z] m += 1 j += 1 print("由于结果较长，请到Excel中查看，文件位于代码同级目录下outdata.csv。") outdata.to_csv('outdata.csv') print("导出成功。") 代码块四：
# 三维坐标可视化显示 out_sv = 'G20' fig = plt.figure() ax = plt.axes(projection='3d') X = outdata[outdata['gnss_sv_id'] == out_sv]['X'] Y = outdata[outdata['gnss_sv_id'] == out_sv]['Y'] Z = outdata[outdata['gnss_sv_id'] == out_sv]['Z'] ax.plot(X, Y, Z, label=out_sv) ax.legend() plt.show() 
