【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

 文章来源地址https://www.toymoban.com/news/detail-849449.html

观测时刻应该由接收机给出，因缺乏该数据，故此处直接使用了星历给出的参考时刻

 

代码如下：

import csv
import math as m
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

with 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 rx


def Rz(fai):
    rz = np.mat([[m.cos(fai), m.sin(fai), 0], [-m.sin(fai), m.cos(fai), 0], [0, 0, 1]])
    return rz

for i in range(len(nfile_lines)):                                     #获取开始行
    if nfile_lines[i].find('C01') != -1:
       start_num = i + 1
       break

for i in range(len(nfile_lines)):                                     #获取终止行
    if nfile_lines[i].find('J01') != -1:
       end_num = i
       break

x_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转为float
        data_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 计算卫星运行的平均角速度n
            GM = 398600441800000
            n0 = m.sqrt(GM) / m.pow(sqrtA, 3)
            n = n0 + delta_n

            #1.2 计算归化时间tk
            if month <= 2:
                year -= 1
                month += 12
            JD = 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秒为BDT
            tp = 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 - toe
            if tk > 302400:
                tk -= 604800
            if tk < -302400:
                tk += 604800  
                
            #1.3 观测时刻卫星平近点角Mk的计算
            Mk = M0 + n * tk

            #1.4 计算偏近点角
            count = 0
            E0 = M0
            Ek = Mk + e * m.sin(E0)
            while abs(Ek - E0) > 1e-10:
                count += 1
                E0 = Ek
                Ek = 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和轨道倾角ik
            uk = Fai_k + sigema_u
            rk = m.pow(sqrtA, 2) * (1 - e * m.cos(Ek)) + sigema_r
            ik = 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-5
            if 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 * tk
                five = -5 * m.pi / 180
                temp = np.mat([[XGk], [YGk], [ZGk]])
                temp = Rz(fi) * Rx(five) * temp
                Xk = 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卫星

到了这里，关于【python】利用广播星历计算BDS卫星的位置的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！