两个 Topos() 位置之间的路径:确定达到给定高度的纬度和经度

Path between two Topos() locations: determine latitude and longitude where a given altitude is reached

我对轨道力学领域还很陌生,目前正在努力解决以下问题,这应该很容易用 Skyfield 解决,但我对所有不同的坐标系和它们之间的转换有点不知所措他们。

我在地球上有一个 Topos 位置和一个 LEO 卫星的 Topos 位置。我正在考虑他们之间的视线。我想确定沿着这条路径的位置的纬度和经度,它与大气的特定层相交。

一个例子是中间层和基于纬度和经度给出的 100 公里左右其属性的现有数据集。交叉点可以让我更好地理解这些属性对与卫星通信的相互作用。

我尝试直接使用 Skyfield 进行操作,但最终得到一个 Apparent 对象,我无法将其转换回地球上的纬度和经度。首先我用三角函数求出地球到100km高度的那个点的距离

然后,我在地球上取了位置,用不变的仰角、方位角来保持路径的方向,最后加上计算出的距离,就到了这个位置。我想我需要一个 Geocentric 对象来使用 subpoint() 以获得所需的纬度,这个位置的经度。

这是我目前拥有的:

from skyfield.api import load, Distance
from skyfield.toposlib import Topos
import numpy as np

ts = load.timescale()

earth_position = Topos('52.230039 N', '4.842402 E', elevation_m=10)
space_position = Topos('51.526200 N', '5.347795 E', elevation_m=625 * 1000)

difference = (space_position - earth_position).at(ts.now()).altaz()

distance_to_height = 100 / np.sin(difference[0].radians)

position = earth_position.at(ts.now()).from_altaz(alt_degrees=difference[0].degrees, az_degrees=difference[1].degrees, distance=Distance(km=distance_to_height))

我已多次浏览文档,偶然发现 frame_latlon(frame) 通用 ICRF 对象,但我不确定如何进一步进行。

用纬度和经度完全三角地尝试它也没有产生预期的结果。

不幸的是,我确实没有任何可用于更轻松地解决此问题的经过验证的结果。再用三角函数想象一下,很明显,卫星位置高度的增加会使交叉点的纬度、经度更靠近地球上的位置。降低高度会使这个交叉点更靠近卫星。

这是一个有趣的问题,Skyfield 的 API 没有提供简单的方法来询问;如果您可以通过了解视线与特定高度的交点来概述将要解决的更大问题,那么可以为解决相同问题的未来用户编写解决该问题的例程。

同时:

  1. 为了得到你的脚本 运行 我必须从 api 导入 Distance
  2. 无法识别名称 dis,因此我将其替换为 distance_to_height,希望它是预期的名称。
  3. 调用 ts.now() 时,每次调用的日期和时间都略有不同。虽然脚本运行得如此之快以至于它可能无关紧要,但为了清楚起见,我已经转向在脚本开头仅调用一次 now(),这也比重复调用它稍微快一些。 (实际上在这种情况下它要快得多,因为旋转矩阵只计算一次而不是必须为每个单独的时间对象重新计算一次,但这是一个不容易看到的隐藏细节。)
  4. 我怀疑你的几何有问题:100 / sin() 机动只有在地球是平的情况下才有效,我想?但也许你总是在处理几乎在头顶上方的卫星,所以错误是可控的? (或者也许我对几何的想象有误;如果数学实际上是正确的,请随时提供图表。)
  5. 为了便于阅读,我将给 altaz() 的组件命名而不是数字。

有了这些调整, 我认为答案是您需要手动构造一个 Geocentric 通过将观察者的位置加在一起来定位 以及您创建的相对向量 在观察者和沿着视线的那种 100 公里点之间。 必须采取这样的手动步骤表明一个可能的区域 Skyfield 可以改进的地方。 这是它在代码中的样子:

from skyfield.api import load, Distance
from skyfield.positionlib import Geocentric
from skyfield.toposlib import Topos
import numpy as np

ts = load.timescale()
t = ts.now()

earth_position = Topos('52.230039 N', '4.842402 E', elevation_m=10)
space_position = Topos('51.526200 N', '5.347795 E', elevation_m=625 * 1000)

alt, az, distance = (space_position - earth_position).at(t).altaz()

distance_to_height = 100 / np.sin(alt.radians)

e = earth_position.at(t)
p = e.from_altaz(alt_degrees=alt.degrees, az_degrees=az.degrees, distance=Distance(km=distance_to_height))

g = Geocentric(e.position.au + p.position.au, t=t)
s = g.subpoint()
print(s)
print(s.elevation.km, '<- warning: 100/sin() did not produce exactly 100')

我看到的结果是:

Topos 52deg 06' 30.0" N 04deg 55' 51.7" E
100.02752954478532 <- warning: 100/sin() did not produce exactly 100

对于未来, 我在 Skyfield TODO.rst 文件中添加了一些想法 一起可能会朝着解锁更惯用的方式迈进 将来进行这种计算—— 尽管我怀疑甚至还需要执行更多步骤:

https://github.com/skyfielders/python-skyfield/commit/ba1172a0ccfef84473436d9d7b8a7d7011344cbd