跳转至

Components

这里选择 SkyField ToC 中的重点章节进行整理

天文学概念一概忽略，重点聚焦卫星网络与无线网络的内容 🌟 🛰️

Data Files

在 skyfield.api 中，load()是用来拉取的 👀

最基础使用

第一次使用，按照下面的代码走，你会得到对应的输出:

Python
1
2
3
4
5
from skyfield.api import load

# fetch from Internet or LocalCache
planets = load('de421.bsp')
print('Ready')

稍等十几秒，得到输出:

Python
1
2
[#################################] 100% de421.bsp
Ready

检查项目文件夹:

Bash
1
2
❯ ls
 de421.bsp  sun_noon.py  venus.py

发现多了一个de421.bsp

再次运行，立马输出:

Python
1
Ready
Note

本质上，load()指令会从 网络库 / 本地缓存文件 中读取，然后加载得到结果

只有第一次从网络库中拉取，后续都会从本地的缓存文件中读取

当然，如果你后续把缓存文件(de421.bsp)删去，显然它还是会从网络库 fetch

指定缓存文件加入的环境

使用 loader() 指令load的路径

Python
1
2
from skyfield.api import Loader
load = Loader('~/skyfield-data')

此文件后续的所有load()对应的加载与缓存，都会在上述提供的目录下进行

Tip

比如可以将 Satellite 的 TLE 文件放进这个文件夹中

进度条设置

Python
1
2
from skyfield.api import Loader
load = Loader('~/skyfield-data', verbose=False)

这里 verbose=False 就是关闭进度条显示的意思

Earth Satellites

🔥 这一部分非常非常重要 🔥

Skyfield 使用 标准SGP4轨道元素 来预测位置并进行一系列仿真

轨道元素来源: Celestrak

What is TLE

TLE（Two-Line Element Set）是描述卫星轨道的一种数据格式，广泛应用于航天和遥感领域。

TLE数据通常由两个文本行组成，包含了卫星的基本轨道参数，这些参数用于计算和预测卫星在地球周围的运动状态。

关键要素:

  1. 轨道倾角(i)
  2. 升交点赤经(ø)
  3. 偏心率(e)
  4. 近地点幅角(w)
  5. 平近点角(M0)
  6. 平均运动(n0)

Format Example:

Text Only
1
2
3
ISS (ZARYA)
1 25544U 98067A   24127.82853009  .00015698  00000+0  27310-3 0  9995
2 25544  51.6393 160.4574 0003580 140.6673 205.7250 15.50957674452123

目前 Skyfield 支持 TLE / JSON / CSV 格式

既可以用“开箱即用的数据供应”，也可以用“特定查询”

Warning

这篇笔记全程以“TLE”文件格式为例，原因很简单，STK的TLE导入非常方便！

Download satellite elements

Either is ok!

  1. 下载数据在本地run之前，后续就是缓存了
  2. 卫星数据实时更新，因此更新换代特别快，记得按需fetch and refresh
  3. 默认按URL来fetch的数据都使用完全相同的文件名，建议 use the filename= optional argument to create a separate local file for each remote URL

下载TLE文件的范式

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
from skyfield.api import load

max_days = 7.0         # download again once 7 days old
name = 'stations.tle'  # custom filename, not 'gp.php'; use TLE format

base = 'https://celestrak.org/NORAD/elements/gp.php'
url = base + '?GROUP=stations&FORMAT=tle' # object && TLE format

if not load.exists(name) or load.days_old(name) >= max_days:
    load.download(url, filename=name) # customize file name to avoid conflicts

效果显示:

alt text

Load satellite elements

下载元素文件后，将它们加载到Skyfield。对于传统的TLE格式：

加载TLE文件的范式

Python
1
2
3
4
5
6
7
8
9
from skyfield.api import load
from skyfield.iokit import parse_tle_file

ts = load.timescale() # time scale is needed!

with load.open('stations.tle') as f:
    satellites = list(parse_tle_file(f, ts))

print('Loaded', len(satellites), 'satellites')

效果显示:

alt text

Load satellite data from a string

如果已经把TLE/CSV读入内存，并不是从文件读入，那么读取方式是:

Python
1
2
3
# For TLE and JSON:
from io import BytesIO
f = BytesIO(byte_string)

Index by name or num

可以用for循环将内容加载进dictionary，然后在dict中按name/num进行索引

按名字

Python
1
2
3
4
5
# insert sat's name into a dict
by_name = {sat.name: sat for sat in satellites}
# get obj by its name
satellite = by_name['ISS (ZARYA)']
print(satellite)

结果:

Text Only
1
ISS (ZARYA) catalog #25544 epoch 2024-05-09 08:48:20 UTC

按数字

Python
1
2
3
4
5
# insert sat's num into a dict
by_number = {sat.model.satnum: sat for sat in satellites}
# get obj by its satnum
satellite = by_number[25544]
print(satellite)

结果:

Text Only
1
ISS (ZARYA) catalog #25544 epoch 2024-05-09 08:48:20 UTC

Load a single TLE set from strings

实际上， 我们的Satellite是拿TLE配置进行实例化的

在之前的流程中，我们是下载了TLE并用它设计对象，其实我们也可以在Py脚本中直接手动输入TLE配置String，然后实例化 Satellite 对象

Python
1
2
3
4
5
6
7
from skyfield.api import EarthSatellite

ts = load.timescale()
line1 = '1 25544U 98067A   14020.93268519  .00009878  00000-0  18200-3 0  5082'
line2 = '2 25544  51.6498 109.4756 0003572  55.9686 274.8005 15.49815350868473'
satellite = EarthSatellite(line1, line2, 'ISS (ZARYA)', ts)
print(satellite)
Initiate a Satellite Object

实例化一个 Satellite 对象 需要提供的参数:

  • arg1: line1 of TLE
  • arg2: line2 of TLE
  • objName: str name of Satellite
  • ts: timescale

后两个参数是optional的

Find when a satellite rises and sets

我们可以在指定时间区间内，搜索 卫星高度超过地平线高度数量 的每种情况:

重点用法

Python
1
satellite.find_events(...)

示例解析

在一天中的高度30°高以上的有多少次

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
from skyfield.api import wgs84

# 导入 wgs84 模块，这是用于处理地球坐标系统的模块
# (arg0: 纬度, arg1: 经度)
bluffton = wgs84.latlon(+40.8939, -83.8917)

# 创建一个地面观测点
t0 = ts.utc(2014, 1, 23) # start_time
t1 = ts.utc(2014, 1, 24) # end_time
t, events = satellite.find_events(bluffton, t0, t1, altitude_degrees=30.0)

event_names = 'rise above 30°', 'culminate', 'set below 30°'

for ti, event in zip(t, events):
    name = event_names[event]
    print(ti.utc_strftime('%Y %b %d %H:%M:%S'), name)

结果:

Text Only
1
2
3
4
5
6
2014 Jan 23 06:25:37 rise above 30°
2014 Jan 23 06:26:58 culminate
2014 Jan 23 06:28:19 set below 30°
2014 Jan 23 12:54:56 rise above 30°
2014 Jan 23 12:56:27 culminate
2014 Jan 23 12:57:58 set below 30°

Generate a satellite position

地心坐标系表示的位置

生成卫星位置的最简单形式是调用 at() 方法, return an (x,y,z) position relative to the Earth's center in the Geocentric Celestial Reference System

Python
1
2
3
4
5
# You can instead use ts.now() for the current time
t = ts.utc(2014, 1, 23, 11, 18, 7)

geocentric = satellite.at(t)
print(geocentric.position.km)

alt text

Satellite latitude, longitude, and height

纬度、经度、高度 表示的位置

  • 返回经纬度: latlon_of()
  • 返回高度: height_of()
  • 返回点的位置: geographic_position_of()
Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from skyfield.api import EarthSatellite, wgs84
from skyfield.api import load

ts = load.timescale()
line1 = '1 25544U 98067A   14020.93268519  .00009878  00000-0  18200-3 0  5082'
line2 = '2 25544  51.6498 109.4756 0003572  55.9686 274.8005 15.49815350868473'
satellite = EarthSatellite(line1, line2, 'ISS (ZARYA)', ts)

t = ts.utc(2014, 1, 23, 11, 18, 7)

print("================= by at()  ===================")

geocentric = satellite.at(t)
print(geocentric.position.km)

print("================= by wgs84 ===================")

lat, lon = wgs84.latlon_of(geocentric)
print('Latitude:', lat)
print('Longitude:', lon)

print("----------------------------------------------")

height = wgs84.height_of(geocentric)
print('Height:', height.km)

print("----------------------------------------------")

pos = wgs84.geographic_position_of(geocentric)
print('Position of point:', pos)

结果:

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
❯ python3 sat_pos.py
================= by at()  ===================
[-3918.87650458 -1887.64838745  5209.08801512]
================= by wgs84 ===================
Latitude: 50deg 14' 37.4"
Longitude: -86deg 23' 23.3"
----------------------------------------------
Height: 420.87372841094077
----------------------------------------------
Position of point: WGS84 latitude +50.2437 N longitude -86.3898 E elevation 420873.7 m

Satellite altitude, azimuth, and distance

卫星相对于“观察者”的位置表示（球坐标系）

本质: 矢量减法

调用方法:

Python
1
2
3
4
5
6
difference = satellite - bluffton
topocentric = difference.at(t)
alt, az, distance = topocentric.altaz()
# alt: 仰角
# az: 方位角
# distance: 观察者与点的直线距离

Satellite right ascension and declination

right ascension: 赤经, RA

declination: 赤纬, DEC

alt text

Python
1
2
3
ra, dec, distance = topocentric.radec()  # ICRF ("J2000")
# or
ra, dec, distance = topocentric.radec(epoch='date')

Satellite’s range distance/rate

行星在做轨道运动，就拿椭圆为例，它的 离中心的距离 和 abs(速度) 一直在变，我们想看看 指定行星/卫星在特定的时间区间内的距离/速率变化 情况:

调用方法

Python
1
2
3
pos.frame_latlon_and_rates(...)
# 第三个返回值: distance变化范围
# 第六个返回值: rate变化范围

示例

Python
1
2
3
4
5
6
7
t = ts.utc(2014, 1, 23, 11, range(17, 23))
pos = (satellite - bluffton).at(t)
_, _, the_range, _, _, range_rate = pos.frame_latlon_and_rates(bluffton)

from numpy import array2string
print(array2string(the_range.km, precision=1), 'km')
print(array2string(range_rate.km_per_s, precision=2), 'km/s')

结果:

Text Only
1
2
[1434.2 1190.5 1064.3 1097.3 1277.4 1553.6] km
[-4.74 -3.24 -0.84  1.9   3.95  5.14] km/s

Find when a satellite is in sunlight

对于 GS(地面站) 而言，只有指定卫星在它对应的仰角内是“sunlight”，才可以观察到

否则就只能采取雷达遥感等追踪技术

Python
1
2
3
4
5
eph = load('de421.bsp')

two_hours = ts.utc(2014, 1, 20, 0, range(0, 120, 20))
sunlit = satellite.at(two_hours).is_sunlit(eph)
print(sunlit)

Find whether the Earth blocks a satellite’s view

有时，对于观测卫星而言，地球会挡住它们想观察的行星/卫星

调用方法

Python
1
sunlit = p.is_behind_earth()

示例

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# 加载的星历被赋值给变量 eph
eph = load('de421.bsp')
# 加载的星历 eph 中提取地球和金星的位置对象
earth, venus = eph['earth'], eph['venus']
# 观察时间: 从2014-1-20-0:00开始的2h
two_hours = ts.utc(2014, 1, 20, 0, range(0, 120, 20))
# 从地球加上卫星定义的位置看到的金星的视位置
p = (earth + satellite).at(two_hours).observe(venus).apparent()
sunlit = p.is_behind_earth()
print(sunlit)

Detect Propagation Errors

遇见error时，将返回 (nan, nan, nan)

示例:

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from skyfield.api import EarthSatellite, load

text = """
GOCE
1 34602U 09013A   13314.96046236  .14220718  20669-5  50412-4 0   930
2 34602 096.5717 344.5256 0009826 296.2811 064.0942 16.58673376272979
"""
lines = text.strip().splitlines()

ts = load.timescale()
sat = EarthSatellite(lines[1], lines[2], lines[0])
print(sat.epoch.utc_jpl())
'''
sat.epoch.utc_jpl() 返回卫星轨道元素的起始时间（epoch），并以 JPL UTC 格式打印
'''

# 计算卫星在特定时间的位置（2013-11-09）
geocentric = sat.at(ts.utc(2013, 11, 9))
print('Before:')
print(geocentric.position.km)
print(geocentric.message)

# 计算卫星在特定时间的位置（2013-11-13）
geocentric = sat.at(ts.utc(2013, 11, 13))
print('\nAfter:')
print(geocentric.position.km)
print(geocentric.message)

结果:

Text Only
1
2
3
4
5
6
7
8
9
❯ python3 pp_error.py
A.D. 2013-Nov-10 23:03:03.9479 UTC
Before:
[nan nan nan]
mean eccentricity is outside the range 0.0 to 1.0

After:
[-5021.82658191   742.71506112  3831.57403957]
mrt is less than 1.0 which indicates the satellite has decayed

本质: propagation error 可以通过比较不同时间点的卫星位置来检测

  1. 位置变化：通过计算卫星在两个不同时间点（2013年11月9日和2013年11月13日）的地心坐标，可以观察到它们之间的变化。
    • 如果这两个位置之间的变化超出了预期范围，可能意味着存在传播误差
  2. 消息输出：使用 geocentric.message 可以获取与计算相关的任何警告或信息。
    • 如果存在传播误差或其他问题，Skyfield 库可能会在这个消息中提供相关信息。
geocentric.message

实际上笔者认为这个geocentric.message就跟logging一样

Python
1
2
3
4
from pprint import pprint

geocentric = sat.at(ts.utc(2013, 11, [9, 10, 11, 12, 13]))
pprint(geocentric.message)

采用 数组(list)导入 的方式，可以清晰地看到message像logging流一样记录状态，并且支持导出

Text Only
1
2
3
4
5
['mean eccentricity is outside the range 0.0 to 1.0',
None,
None,
None,
'mrt is less than 1.0 which indicates the satellite has decayed']

Build a satellite with a specific gravity model

考虑 重力模型

Python
1
2
3
4
5
6
from sgp4.api import Satrec, WGS84
# Satrec 用于表示卫星的轨道元素，而 WGS84 是用于计算卫星位置的重力模型
satrec = Satrec.twoline2rv(line1, line2, WGS84)
# WGS84 参数指定使用 WGS84 重力模型进行计算
sat = EarthSatellite.from_satrec(satrec, ts)
# 初始化对象，使用重力模型构建

Build a satellite from orbital elements

使用 轨道的参数信息 ，构建卫星模型

调用方法

基于 sgp4init() 方法

示例

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
from sgp4.api import Satrec, WGS72
from skyfield.api import EarthSatellite, load

satrec = Satrec()

satrec.sgp4init(
    WGS72,           # gravity model
    'i',             # 'a' = old AFSPC mode, 'i' = improved mode
    5,               # satnum: Satellite number
    18441.785,       # epoch: days since 1949 December 31 00:00 UT
    2.8098e-05,      # bstar: drag coefficient (/earth radii)
    6.969196665e-13, # ndot: ballistic coefficient (radians/minute^2)
    0.0,             # nddot: second derivative of mean motion (radians/minute^3)
    0.1859667,       # ecco: eccentricity
    5.7904160274885, # argpo: argument of perigee (radians)
    0.5980929187319, # inclo: inclination (radians)
    0.3373093125574, # mo: mean anomaly (radians)
    0.0472294454407, # no_kozai: mean motion (radians/minute)
    6.0863854713832, # nodeo: right ascension of ascending node (radians)
)

ts = load.timescale()
sat = EarthSatellite.from_satrec(satrec, ts)

print('Satellite number:', sat.model.satnum)
print('Epoch:', sat.epoch.utc_jpl())

结果:

Text Only
1
2
3
❯ python3 sgp-4.py   
Satellite number: 5
Epoch: A.D. 2000-Jun-27 18:50:24.0000 UTC

根据上述信息，我们可以找到满足上述所有参数特征的对象: (satellite + timestamp)

The result should be a satellite object that behaves exactly as though it had been loaded from TLE lines.

Note

回顾上述内容，构建一个卫星对象，我们有两种方式:

  1. 根据TLE给出的参数，进行 Satellite 对象实例化
  2. 手动给出轨道参数信息，构建 Satellite 对象
    • interact directly with the sgp4 library

使用sgp-4自定义参数进行初始化，很灵活，参数列表: 传送门

TL; DR

事件查询

Python
1
satellite.find_events(...)

地心坐标系位置

Python
1
2
t = ts.utc(2014, 1, 23, 11, 18, 7)
geocentric = satellite.at(t)

经纬度坐标系位置

  • 返回经纬度: latlon_of()
  • 返回高度: height_of()
  • 返回点的位置: geographic_position_of()

俯仰角相对坐标系位置

Python
1
2
3
difference = satellite - bluffton
topocentric = difference.at(t)
alt, az, distance = topocentric.altaz()

赤经纬坐标系位置

Python
1
2
3
ra, dec, distance = topocentric.radec()  # ICRF ("J2000")
# or
ra, dec, distance = topocentric.radec(epoch='date')

距离/速率变化范围

Python
1
pos.frame_latlon_and_rates(bluffton)

是否Sunlight

Python
1
sunlit = satellite.at(two_hours).is_sunlit(eph)

卫星是否被地球遮挡

Python
1
sunlit = p.is_behind_earth()

错误流输出

Python
1
print(geocentric.message)

考虑重力的卫星模型

Python
1
2
3
4
from sgp4.api import Satrec, WGS84
# Satrec 用于表示卫星的轨道元素，而 WGS84 是用于计算卫星位置的重力模型
satrec = Satrec.twoline2rv(line1, line2, WGS84)
# WGS84 参数指定使用 WGS84 重力模型进行计算

基于轨道元素的卫星模型

Python
1
2
3
4
5
6
from sgp4.api import Satrec, WGS72

satrec = Satrec()
satrec.sgp4init(
    ...
)