使用最新的 2020 SGP4 例程跟踪给定 TLE 数据的地球卫星。
项目描述
这个 Python 包计算地球轨道卫星的位置和速度,假设卫星的 TLE 轨道元素来自CelesTrak等来源。它实现了最新版本的 SGP4,并定期针对 SGP4 测试套件运行,以确保其卫星位置预测与算法标准分布的预测一致在 0.1 毫米以内。这个误差远小于卫星本身偏离 TLE 文件中描述的理想轨道的 1-3 公里/天。
如果你的平台支持它,这个包编译并使用来自官方 C++ 版本的 SGP4 的逐字源代码。
否则,将使用 SGP4 的较慢但可靠的 Python 实现。
如果不是一次询问单个卫星的位置,而是为这个库提供一个卫星数组和一个时间数组,那么可以使用机器代码处理这些数组,而不是要求您运行慢速 Python循环它们。
请注意,SGP4 传播器在“真赤道平分点”(TEME) 参考系中返回原始x,y,z笛卡尔坐标,该参考系以地球为中心但不随地球旋转——“以地球为中心的惯性”(ECI) 参考系. SGP4 传播器本身并没有实现将这些位置转换为更官方的 ECI 帧(如 J2000 或 ICRF)的数学运算;也不将位置转换为任何像 ITRS 一样的地心地球固定 (ECEF) 框架;也不通过像 WGS84 这样的地球椭球将它们转换为纬度和经度。
要转换为其他坐标系,请查找构建在此之上的综合天文学库,例如Skyfield库:
https://rhodesmill.org/skyfield/earth-satellites.html
用法
该库使用与官方 C++ 代码相同的函数名称,以帮助可能已经熟悉其他语言的 SGP4 的用户。以下是如何计算 2000 年 6 月 29 日 12:50:19 国际空间站的 x、y、z 位置和速度:
>>> from sgp4.api import Satrec >>> >>> s = '1 25544U 98067A 19343.69339541 .00001764 00000-0 38792-4 0 9991' >>> t = '2 25544 51.6439 211.2001 0007417 17.6667 85.6398 15.50103472202482' >>> satellite = Satrec.twoline2rv(s, t) >>> >>> jd, fr = 2458827, 0.362605 >>> e, r, v = satellite.sgp4(jd, fr) >>> e 0 >>> print(r) # True Equator Mean Equinox position (km) (-6102.44..., -986.33..., -2820.31...) >>> print(v) # True Equator Mean Equinox velocity (km/s) (-1.45..., -5.52..., 5.10...)
作为输入,您可以提供:
jd的简单浮点 Julian Date和fr的值 0.0 ,如果您对 64 位浮点数的精度感到满意。请注意,现代儒略日期大于 2,450,000,这意味着 64 位浮点数的近一半精度将被指定日期的整个部分消耗。其余数字将为大约 20.1 µs 的分数提供精度。这对于您的结果的准确性应该没有问题——无论如何,卫星位置通常会偏离几公里,远远小于卫星在 20.1 µs 内移动的时间——但是如果你运行一个求解器,它在搜索上升或设置时间,求解器可能会受到卫星位置每次跳跃之间 20.1 µs 的高原的困扰。
或者,您可以提供一个粗略的日期jd加上一个非常精确的分数 fr来提供其余的值。计算卫星位置的儒略日是这两个值的总和。一种常见的做法是将整数提供为jd,将小数提供为fr;另一个是让jd携带分数 0.5,因为 UTC 午夜发生在每个儒略日期的中途。无论哪种方式,拆分值都可以让求解器一直运行到纳秒级,并且仍然可以看到 SGP4 对微小的日期调整做出平稳的响应,结果卫星位置的微小变化。
以下是如何解释结果:
如果无法计算给定日期的卫星位置,则e将是一个非零错误代码。您可以从 sgp4.api 导入 SGP4_ERRORS以访问将错误代码映射到解释每个代码含义的错误消息的字典。
r在 SGP4 使用的特殊 True Equator Mean Equinox 坐标系中测量距离地球中心公里的卫星位置。
v速度是位置变化的速率,以千米每秒表示。
如果您的应用程序本身不处理儒略日期,您可以使用jday()从日历日期计算jd和fr。
>>> from sgp4.api import jday >>> jd, fr = jday(2019, 12, 9, 12, 0, 0) >>> jd 2458826.5 >>> fr 0.5
OMM
业界正在进行调整,因为固定宽度的 TLE 格式很快就会用完卫星编号。
一些 TLE 文件现在使用新的“Alpha-5”约定,通过使用首字母来扩展卫星编号的范围;例如,“E8493”表示卫星 148493。这个库现在支持 Alpha-5 约定,并且应该在 Python 中返回正确的整数。
一些当局现在正在以“OMM”轨道平均元素消息格式分发卫星元素,以取代 TLE 格式。您可以在 CelesTrak 网站上的TS Kelso 博士的“获取 GP 数据的新方法”中了解 OMM 。
您已经可以尝试对 OMM 的实验性支持:
>>> from sgp4 import omm
读取 OMM 数据需要两个步骤,因为 OMM 支持几种不同的文本格式。首先,解析输入文本以恢复它存储的字段名称和值;其次,根据这些字段值构建一个 Python 卫星对象。例如,从 XML 加载 OMM:
>>> with open('sample_omm.xml') as f: ... fields = next(omm.parse_xml(f)) >>> sat = Satrec() >>> omm.initialize(sat, fields)
或者,从 CSV 加载 OMM:
>>> with open('sample_omm.csv') as f: ... fields = next(omm.parse_csv(f)) >>> sat = Satrec() >>> omm.initialize(sat, fields)
无论哪种方式,卫星对象都应该正确初始化并准备好开始生成位置。
如果您对使用新的 OMM 格式保存卫星参数感兴趣,请阅读下面的“导出”部分。
时代
在给定卫星的生命周期内,随着其轨道的发展,将产生数十或数百个不同的 TLE 记录。每个 TLE 记录都指定了最准确的“纪元日期”。通常,TLE 仅在其纪元日期前后的几周内有用,超过此时间,其预测将变得不可靠。
卫星对象本身提供它们的纪元作为两位数的年份,然后是一年中的小数天数:
>>> satellite.epochyr 19 >>> satellite.epochdays 343.69339541
由于 Sputnik 于 1957 年发射,卫星元素集永远不会指较早的年份,因此 57 到 99 年表示 1957-1999 年,而 0 到 56 年表示 2000-2056 年。TLE 格式大概会在 2057 年过时,必须升级到 4 位数年份。
要将天数及其分数转换为日历日期和时间,请使用days2mdhms()函数。
>>> from sgp4.api import days2mdhms >>> month, day, hour, minute, second = days2mdhms(19, 343.69339541) >>> month 12 >>> day 9 >>> hour 16 >>> minute 38 >>> second 29.363424
SGP4 库还将这两个数字转换为儒略日期和小数儒略日期,因为儒略日期在天文学中更常用。
>>> satellite.jdsatepoch 2458826.5 >>> satellite.jdsatepochF 0.69339541
最后,如果您需要 epoch 日期和时间作为 Python datetime,库中提供了一个便利函数。
>>> from sgp4.conveniences import sat_epoch_datetime >>> sat_epoch_datetime(satellite) datetime.datetime(2019, 12, 9, 16, 38, 29, 363423, tzinfo=UTC)
阵列加速
当你有很多日期时,为了避免 Python 循环的开销,你可以将它们作为数组传递给另一个理解 NumPy 的方法:
>>> import numpy as np >>> np.set_printoptions(precision=2)
>>> jd = np.array((2458826, 2458826, 2458826, 2458826)) >>> fr = np.array((0.0001, 0.0002, 0.0003, 0.0004))
>>> e, r, v = satellite.sgp4_array(jd, fr)
>>> print(e) [0 0 0 0] >>> print(r) [[-3431.31 2620.15 -5252.97] [-3478.86 2575.14 -5243.87] [-3526.09 2529.89 -5234.28] [-3572.98 2484.41 -5224.19]] >>> print(v) [[-5.52 -5.19 1.02] [-5.49 -5.22 1.08] [-5.45 -5.25 1.14] [-5.41 -5.28 1.2 ]]
当您有许多卫星和日期时,为了避免 Python 循环的开销,请从多个单独的卫星构建一个SatrecArray 。它的 sgp4()方法期望jd和fr都是 NumPy 数组,因此如果您只有一个日期,请务必提供长度为 1 的 NumPy 数组。这是 2 颗卫星和 4 个日期的示例计算:
>>> u = '1 20580U 90037B 19342.88042116 .00000361 00000-0 11007-4 0 9996' >>> w = '2 20580 28.4682 146.6676 0002639 185.9222 322.7238 15.09309432427086' >>> satellite2 = Satrec.twoline2rv(u, w)
>>> from sgp4.api import SatrecArray >>> a = SatrecArray([satellite, satellite2]) >>> e, r, v = a.sgp4(jd, fr)
>>> np.set_printoptions(precision=2) >>> print(e) [[0 0 0 0] [0 0 0 0]] >>> print(r) [[[-3431.31 2620.15 -5252.97] [-3478.86 2575.14 -5243.87] [-3526.09 2529.89 -5234.28] [-3572.98 2484.41 -5224.19]] <BLANKLINE> [[ 5781.85 2564. -2798.22] [ 5749.36 2618.59 -2814.63] [ 5716.35 2672.94 -2830.78] [ 5682.83 2727.05 -2846.68]]] >>> print(v) [[[-5.52 -5.19 1.02] [-5.49 -5.22 1.08] [-5.45 -5.25 1.14] [-5.41 -5.28 1.2 ]] <BLANKLINE> [[-3.73 6.33 -1.91] [-3.79 6.3 -1.88] [-3.85 6.28 -1.85] [-3.91 6.25 -1.83]]]
出口
如果你有一个Satrec想要与朋友分享或保存到一个文件中,有一个导出例程可以将它转回 TLE:
>>> from sgp4 import exporter >>> line1, line2 = exporter.export_tle(satellite) >>> line1 '1 25544U 98067A 19343.69339541 .00001764 00000-0 38792-4 0 9991' >>> line2 '2 25544 51.6439 211.2001 0007417 17.6667 85.6398 15.50103472202482'
令人高兴的是,这正是我们用来创建此卫星对象的两条 TLE 行:
>>> (s == line1) and (t == line2) True
另一个导出例程可用,它生成由新 OMM 格式定义的字段(请参阅上面的“OMM”部分):
>>> from pprint import pprint >>> fields = exporter.export_omm(satellite, 'ISS (ZARYA)') >>> pprint(fields) {'ARG_OF_PERICENTER': 17.6667, 'BSTAR': 3.8792e-05, 'CENTER_NAME': 'EARTH', 'CLASSIFICATION_TYPE': 'U', 'ECCENTRICITY': 0.0007417, 'ELEMENT_SET_NO': 999, 'EPHEMERIS_TYPE': 0, 'EPOCH': '2019-12-09T16:38:29.363423', 'INCLINATION': 51.6439, 'MEAN_ANOMALY': 85.6398, 'MEAN_ELEMENT_THEORY': 'SGP4', 'MEAN_MOTION': 15.501034720000002, 'MEAN_MOTION_DDOT': 0.0, 'MEAN_MOTION_DOT': 1.764e-05, 'NORAD_CAT_ID': 25544, 'OBJECT_ID': '1998-067A', 'OBJECT_NAME': 'ISS (ZARYA)', 'RA_OF_ASC_NODE': 211.2001, 'REF_FRAME': 'TEME', 'REV_AT_EPOCH': 20248, 'TIME_SYSTEM': 'UTC'}
重力
SGP4 算法在一组指定地球引力强度的常数之上运行。关于 SGP4 的最新官方论文(见下文)指定“我们使用 WGS-72 作为默认值”,因此这个 Python 模块使用相同的默认值。但是,如果您想使用旧版本的 WGS-72 常量,或者非标准但更现代的 WGS-84 常量,则twoline2rv()构造函数采用可选参数:
>>> from sgp4.api import WGS72OLD, WGS72, WGS84 >>> satellite3 = Satrec.twoline2rv(s, t, WGS84)
如果您选择 WGS-84,您通常会得到不太准确的结果。尽管它反映了更近期和更准确的地球测量结果,但整个行业的卫星 TLE 最有可能以 WGS-72 为基础生成。如果您使用与生成 TLE 相同的基础重力常数,您生成的位置将更好地与每颗卫星的实际位置一致。
提供自己的元素
如果您想直接指定轨道元素而不是解析 TLE,则可以将它们作为浮点数传递给卫星对象的sgp4init()方法。例如,以下是我们在上面的第一个代码示例中从 TLE 加载的相同国际空间站轨道的构建方法:
>>> satellite2 = Satrec() >>> satellite2.sgp4init( ... WGS72, # gravity model ... 'i', # 'a' = old AFSPC mode, 'i' = improved mode ... 25544, # satnum: Satellite number ... 25545.69339541, # epoch: days since 1949 December 31 00:00 UT ... 3.8792e-05, # bstar: drag coefficient (1/earth radii) ... 0.0, # ndot: ballistic coefficient (revs/day) ... 0.0, # nddot: mean motion 2nd derivative (revs/day^3) ... 0.0007417, # ecco: eccentricity ... 0.3083420829620822, # argpo: argument of perigee (radians) ... 0.9013560935706996, # inclo: inclination (radians) ... 1.4946964807494398, # mo: mean anomaly (radians) ... 0.06763602333248933, # no_kozai: mean motion (radians/minute) ... 3.686137125541276, # nodeo: R.A. of ascending node (radians) ... )
这些数字看起来与 TLE 中的数字不同,因为底层sgp4init()例程使用不同的单位:弧度而不是度数。但这是相同的轨道,将产生相同的位置。
请注意, SGP4 传播器会忽略ndot和nddot,因此您可以将它们保留为0.0,而不会对生成的卫星位置产生任何影响。但它们至少会保存到卫星对象中,并且如果您将参数写入 TLE 或 OMM 文件(请参阅上面的“导出”部分),它们会被写出。
要计算“纪元”参数,请取纪元的儒略日期并减去 2433281.5 天。
虽然底层sgp4init()例程未设置属性 epochyr、epochdays、jdsatepoch和jdsatepochF,但该库继续使用您提供的 epoch 为您设置它们。
有关卫星记录初始化后可用属性的完整列表,请参见下一节。
属性
有几十个Satrec属性公开来自底层 C++ SGP4 记录的数据。它们分为以下几类。
鉴别
这些是直接从 TLE 记录复制的,但不被传播数学使用。
轨道元素
这些是轨道参数,从 TLE 记录的文本中逐字复制。它们描述了 TLE 纪元时刻的轨道,因此即使卫星记录被一遍又一遍地用于传播不同时间的位置,它们仍然保持不变。
您还可以将纪元作为儒略日期访问:
计算轨道属性
这些是在首次加载卫星时计算的,以方便可能对它们感兴趣的呼叫者。SGP4 传播器本身不使用它们。
传播者模式
最近传播的结果
可能的错误代码是:
没有错误。
平均偏心率在 0 ≤ e < 1 范围之外。
平均运动已降至零以下。
扰动偏心率在 0 ≤ e ≤ 1 范围之外。
轨道的半圆直肠的长度已降至零以下。
(不再使用。)
轨道已衰减:计算的位置在地下。(仍然返回位置,以防矢量对可能正在搜索重新进入的时刻的软件有帮助。)
最近传播的平均元素
在每次传播的中途,SGP4 例程保存一组“单平均平均元素”,这些元素描述了计算位置时的轨道形状。它们是相对于平均异常进行平均的,包括长期重力、大气阻力以及——在深空模式下——来自太阳和月球的那些扰动的影响,SGP4 对每个天体的整个旋转进行平均。他们忽略了 SGP4 在计算每个位置之前应用的来自太阳和月亮的短期和长期周期性扰动。
重力模型参数
初始化卫星记录时,您选择的重力模型会导致将八个常量的列表复制到:
打印卫星属性
如果你想打印一个卫星,这个库提供了一个方便的“属性转储”例程,它接收一个卫星并生成列出其属性的行:
from sys import stdout from sgp4.conveniences import dump_satrec stdout.writelines(dump_satrec(satellite))
如果要比较两个卫星,则只需传递第二个参数;第二颗卫星的属性将打印在第一颗卫星旁边的第二列中。
stdout.writelines(dump_satrec(satellite, satellite2))
根据官方算法进行验证
此实现通过了 Vallado 等人在 2010 年 8 月发布的 SGP4 参考实现中的所有自动化测试,他们最初于 2006 年发布了他们的 SGP4 修订版:
Vallado、David A.、Paul Crawford、Richard Hujsak 和 TS Kelso,“重温 Spacetrack 报告 #3”,发表在 AIAA/AAS 天体动力学专家会议上,科罗拉多州 Keystone,2006 年 8 月 21-24 日。
如果您想查看该论文,可以在线获取。您可以随时在AIAA-2006-6753.zip下载他们的代码的最新版本,以便与这个 Python 模块(或其他实现)进行比较。
对于开发人员
开发人员可以从 GitHub 上查看这个完整的项目:
https://github.com/brandon-rhodes/python-sgp4
要运行其单元测试,请安装 Python 2、Python 3 和tox 测试工具。在 Python 2 中运行的测试将运行例程的后备纯 Python 版本,而 Python 3 运行快速的新 C++ 加速代码:
cd python-sgp4 tox
旧版 API
在这个库转向包装 Valalado 的官方 C++ 代码并仅在纯 Python 中运行之前,它有一个稍微古怪的 API,仍然支持与旧客户端兼容。您可以通过阅读 1.4 或更早版本的文档来了解它:
变更日志
2022-04-06 — 2.21
将dump_satrec ()添加到sgp4.conveniences模块。
修复了Satrec属性.error,该属性以前从内存中的错误数据构建一个无意义的整数。
从 Python Satrec中删除了.whichconst,以帮助用户避免编写在 C++ 扩展可用时会中断的代码。
2021-07-01 — 2.20
教导sgp4init ()将epochdays和jdsatepochF 舍入 到用于 TLE 中日期分数的相同 8 位小数,如果用户提供的epoch本身在小数点后面有 8 位或更少位。这应该可以更轻松地构建以完美精度往返 TLE 格式的卫星。
修复了export_tle()在 BSTAR 字段的值(如果以科学计数法编写)具有正指数时如何格式化该字段。
修复了由sgp4init ()分配的纪元,因此 2000 年之前的年份有两位而不是三位数(例如,1980 年产生的纪元为 80 而不是 980)。
2021-04-22 — 2.19
扩展了 Python 包索引和模块文档字符串中的文档,因此它列出了该库公开的每个Satrec属性;即使是更晦涩的那些也可能对致力于分析卫星轨道的人们有用。
2021-03-08 — 2.18
如果 TLE 卫星编号缺少所需的 5 位数字, twoline2rv()现在为底层 C++ 库提供一点帮助,因此它仍然可以正确解析分类和国际指示符。
Satrec属性jdsatepoch、jdsatepochF、 epochyr和epochdays现在是可写的,因此用户可以手动调整它们的值——这应该可以弥补 sgp4init()方法无法将它们设置为全浮点精度的事实。