发表更新astronomy / 天文学15 分钟读完 (大约2238个字)0次访问

基于Skyfield&Orekit中SGP4轨道位置计算及卫星网络仿真中应用

本文介绍如何基于Skyfield&Orekit中SGP4轨道位置计算并简单讨论其在卫星网络仿真中应用。

SGP4

TLE(Two-line element set)轨道根数对应的计算模型是简化普适模型(simplified perturbations models,包括SGP, SGP4, SDP4, SGP8 and SDP8),SGP(Simplified General Perturbations)简化摄动模型适用轨道周期小于225分钟的近地物体,SDP(Simplified Deep Space Perturbations)简化深空摄动模型适用于轨道周期大于 225 分钟的物体。

SGP4 模型星历误差大约1km并且每天增长大约1–3km误差。最初的 SGP 模型由 Kozai 于 1959 年开发，由 Hilton & Kuhlman 于 1966 年完善，最初由National Space Surveillance Control Center（以及后来的United States Space Surveillance Network)）用于跟踪轨道上的物体，SDP4模型在历元上的误差为10公里。

SGP4模型具体算法请自行研究，本人不做天文研究仅利用该模型做位置预测，卫星网络仿真可能能应用上。

Skyfield&Orekit

Skyfield是一个优雅的天文学Python软件包，Skyfield 可用于计算恒星、行星、 以及绕地球轨道运行的卫星。 Skyfield中提供了SGP4模型相关API，基于相关api可用来预测卫星位置。Skyfield根据MIT License开源分发。

Orekit)是一个用 Java 编写的低级空间动力学库，自 2008 年以开源许可证发布以来，获得了广泛的认可。Orekit 旨在为飞行动力学应用的开发提供准确、高效的低级组件。Orekit 根据Apache License version 2.0开源分发。另外社区维护了一个 Python 包装器orekit-python-wrapper，可基于Python使用Orekit API。

坐标系

对卫星进行轨道位置计算需要使用到如下天文相关坐标系统：
ECI(Earth-Centered, Inertial/地心惯性坐标系)坐标系统，ECEF(Earth-centered, Earth-fixed/地心地固坐标系)坐标系统及WGS84经纬度。

ECEF与ECI坐标系有着相同的坐标原点和z轴定义,但ECEF坐标系是与地球保持同步旋转。

ITRS(International Terrestrial Reference System/国际地球参考系统)是一种ECEF，不受地球自转、地球赤道扭曲等影响，因此可以更精确地反映地球表面的物理位置。

ITRF(International Terrestrial Reference Frame/国际地球参考框架是ITRS的一个实现。

地表物理位置唯一的WGS84经纬度坐标可转化为唯一的ITRS坐标，通常可用于定位卫星地面站。

PS：具体坐标系相关详细标准请自行查阅。

Skyfield SGP4轨道位置计算

如下以STARLINK-1007卫星为例，计算当前卫星ECI坐标，ECEF坐标及WGS84经纬度。

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#!/usr/bin/python
import time 
from sgp4.api import Satrec, WGS84 
from skyfield.api import EarthSatellite
from skyfield.api import load as skyfield_load 
from skyfield.api import wgs84
from skyfield import framelib

ts = skyfield_load.timescale(builtin=True)

class Satellite():
    id = 0
    name = ""
    tle_line1 = ""
    tle_line2 = ""
    sat = None
    sw_name = ""
    sw_dpid = ""

    def __init__(self, id, line0, line1, line2):
        self.id = id
        self.name = line0
        self.tle_line1 = line1
        self.tle_line2 = line2

        satrec = Satrec.twoline2rv(line1, line2, WGS84)
        self.sat = EarthSatellite.from_satrec(satrec, ts)
        print(self.sat)
      
    def frame_now(self):
        sat_obj = self.sat.at(ts.now())
        # ECI coords are True Equator and Equinox of date (Geocentric Celestial Reference System)
        position = sat_obj.position.km
        velocity = sat_obj.velocity.km_per_s

        # Add ECEF values to instrument.
        # ECEF coords are International Terrestrial Reference System (ITRS)
        ecef = sat_obj.frame_xyz_and_velocity(framelib.itrs)

        # Convert to geocentric latitude, longitude, altitude.
        lat, lon = wgs84.latlon_of(sat_obj)
        alt = wgs84.height_of(sat_obj).km

        return {'position_eci_x': position[0],
                'position_eci_y': position[1],
                'position_eci_z': position[2],
                'velocity_eci_x': velocity[0],
                'velocity_eci_y': velocity[1],
                'velocity_eci_z': velocity[2],
                'position_ecef_x': ecef[0].km[0],
                'position_ecef_y': ecef[0].km[1],
                'position_ecef_z': ecef[0].km[2],
                'velocity_ecef_x': ecef[1].km_per_s[0],
                'velocity_ecef_y': ecef[1].km_per_s[1],
                'velocity_ecef_z': ecef[1].km_per_s[2],
                'geod_latitude': lat.degrees,
                'geod_longitude': lon.degrees,
                'geod_altitude': alt}

if __name__ == "__main__":
    line0 = "STARLINK-1007"
    line1 = "1 44713U 19074A   24120.93093052  .00013004  00000+0  89012-3 0  9990"
    line2 = "2 44713  53.0547 353.6824 0001286  89.2500 270.8636 15.06406908246617"

    sat = Satellite(0, line0, line1, line2)
    for i in range(0, 10):
        frame_data = sat.frame_now()

        x = frame_data['position_ecef_x']
        y = frame_data['position_ecef_y']
        z = frame_data['position_ecef_z']

        vx = frame_data['velocity_ecef_x']
        vy = frame_data['velocity_ecef_y']
        vz = frame_data['velocity_ecef_z']

        lon = frame_data['geod_longitude']
        lat = frame_data['geod_latitude']
        alt = frame_data['geod_altitude']

        timestamp = time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(time.time()))
        print("tm", timestamp, "\n\tx:", x, "y:", y, "z:", z, "\n\tvx:", vx, "vy:", vy, "vz:", vz, "\n\tlon:", lon, "lat:", lat, "alt:", alt)
        time.sleep(10)

计算结果与n2yo.com对比：
skyfield-n2yo

Orekit轨道位置计算

同样以STARLINK-1007卫星为例，计算当前卫星ECEF坐标(ITRS)及WGS84经纬度，基于orekit 11.2和hipparchus 3.0版本。

Java Demo

maven pom依赖如下：

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
31
32
33
34
35
36
37
38
<dependencies>  
    <dependency>  
        <groupId>org.orekit</groupId>  
        <artifactId>orekit</artifactId>  
        <version>11.2</version>  
    </dependency>  
  
    <dependency>  
        <groupId>org.hipparchus</groupId>  
        <artifactId>hipparchus-core</artifactId>  
        <version>3.0</version>  
    </dependency>  
  
    <dependency>  
        <groupId>org.hipparchus</groupId>  
        <artifactId>hipparchus-fitting</artifactId>  
        <version>3.0</version>  
    </dependency>  
  
    <dependency>  
        <groupId>org.hipparchus</groupId>  
        <artifactId>hipparchus-optim</artifactId>  
        <version>3.0</version>  
    </dependency>  
  
    <dependency>  
        <groupId>org.hipparchus</groupId>  
        <artifactId>hipparchus-ode</artifactId>  
        <version>3.0</version>  
    </dependency>  
  
    <dependency>  
        <groupId>org.hipparchus</groupId>  
        <artifactId>hipparchus-geometry</artifactId>  
        <version>3.0</version>  
    </dependency>  
  
</dependencies>

java SGP4计算demo:

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

import org.hipparchus.geometry.euclidean.threed.Vector3D;  
import org.hipparchus.util.FastMath;  
  
import java.io.File;  
import java.util.Calendar;  
import java.util.Date;  
import java.util.TimeZone;  
import java.util.concurrent.TimeUnit;  
  
import org.orekit.bodies.GeodeticPoint;  
import org.orekit.bodies.OneAxisEllipsoid;  
import org.orekit.data.DataContext;  
import org.orekit.data.DataProvidersManager;  
import org.orekit.data.DirectoryCrawler;  
import org.orekit.errors.OrekitException;  
import org.orekit.frames.Frame;  
import org.orekit.frames.FramesFactory;  
  
import org.orekit.propagation.SpacecraftState;  
import org.orekit.propagation.analytical.tle.TLE;  
import org.orekit.propagation.analytical.tle.TLEPropagator;  
import org.orekit.time.AbsoluteDate;  
import org.orekit.time.TimeScalesFactory;  
import org.orekit.utils.Constants;  
import org.orekit.utils.IERSConventions;  
import org.orekit.utils.PVCoordinates;  
  
  
public class OrekitSgp4Demo {  
    // orekit-data文件名,这里将orekit-data解压放到resources文件目录下  
	public static final String OREKIT_PATH_NAME = "orekit-data";  
  
    private static void initOreKit() throws OrekitException {  
        String resPath = OrekitSgp4Demo.class.getClassLoader().getResource("").getPath();  
        String orekitPhysicalDataLocation = resPath + "/" + OREKIT_PATH_NAME;  
  
        System.out.println("orekit path: " + orekitPhysicalDataLocation);  
  
        File orekitData = new File(orekitPhysicalDataLocation);  
        if (!orekitData.exists()) {  
            throw new RuntimeException(  
                    "Unrecoverable error: Could not locate orekit physical data at default location of '"  
                            + orekitPhysicalDataLocation + "'!");  
        }  
        DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();  
        manager.addProvider(new DirectoryCrawler(orekitData));  
    }  
  
    public static TLE buildTLEFromLines(String tleLine1, String tleLine2) throws OrekitException {  
        // Always test if TLE is good first  
        if (!TLE.isFormatOK(tleLine1, tleLine2)) {  
            throw new RuntimeException("Bad TLE format!\n" + tleLine1 + '\n' + tleLine2 + '\n');  
        }  
        return new TLE(tleLine1, tleLine2);  
    }  
  
    public static String tle1Line0 = "STARLINK-1007";  
    public static String tle1Line1 = "1 44713U 19074A   24120.93093052  .00013004  00000+0  89012-3 0  9990";  
    public static String tle1Line2 = "2 44713  53.0547 353.6824 0001286  89.2500 270.8636 15.06406908246617";  
  
    public static void main(String[] args) throws InterruptedException {  
        // init orekit  
        initOreKit();  
  
        // build tle  
        TLE tle = buildTLEFromLines(tle1Line1, tle1Line2);  
        TLEPropagator tleProp = TLEPropagator.selectExtrapolator(tle);  
  
        // get ITRF frame  
        Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);  
        OneAxisEllipsoid oae = new OneAxisEllipsoid(  
                Constants.WGS84_EARTH_EQUATORIAL_RADIUS,  
                Constants.WGS84_EARTH_FLATTENING,  
                itrf);  
  
        for (int i = 0; i < 10; i++) {  
  
            Date now = Calendar.getInstance(TimeZone.getTimeZone("GMT+:08:00")).getTime();  
            AbsoluteDate absoluteDate = new AbsoluteDate(now, TimeScalesFactory.getUTC());  
            SpacecraftState state = tleProp.propagate(absoluteDate);  
            PVCoordinates pvInITRF = state.getPVCoordinates(itrf);  
            Vector3D position = pvInITRF.getPosition();  
            Vector3D velocity = pvInITRF.getVelocity();  
  
            GeodeticPoint gp = oae.transform(position, itrf, absoluteDate);  
            double longitude = FastMath.toDegrees(gp.getLongitude());  
            double latitude = FastMath.toDegrees(gp.getLatitude());  
            double altitude = gp.getAltitude();  
  
            System.out.println(  
                    "time:" + now +  
                            " \n\tx:"  
                            + position.getX()  
                            + "m, y:"  
                            + position.getY()  
                            + "m, z:"  
                            + position.getZ()  
                            + "m, \n\tvx:"  
                            + velocity.getX()  
                            + "m/s, vy:"  
                            + velocity.getY()  
                            + "m/s, vz:"  
                            + velocity.getZ()  
                            + "m/s, \n\tlon:"  
                            + longitude  
                            + "°, lat:"  
                            + latitude  
                            + "°, alt:"  
                            + altitude  
                            + "m"  
                              );  
  
            TimeUnit.SECONDS.sleep(10);  
        }  
    }  
}

计算结果与n2yo.com对比：
orekit-java-n2yo

Python Demo

基于Python调用Orekit最好使用conda安装环境，如下：

1
2
3
4
conda create -y -n orekit_env -c conda-forge orekit=11.2
conda activate orekit_env

pip install numpy
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#!/usr/bin/python
import orekit
from orekit.pyhelpers import setup_orekit_curdir

from org.orekit.frames import FramesFactory, TopocentricFrame
from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint

from org.orekit.utils import IERSConventions, Constants

import time
import numpy as np

from datetime import datetime, timezone

from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
from org.orekit.time import AbsoluteDate, TimeScalesFactory

vm = orekit.initVM()
setup_orekit_curdir("./orekit-data.zip")

itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                         Constants.WGS84_EARTH_FLATTENING, 
                         itrf)


if __name__ == "__main__":
    line0 = "STARLINK-1007"
    line1 = "1 44713U 19074A   24120.93093052  .00013004  00000+0  89012-3 0  9990"
    line2 = "2 44713  53.0547 353.6824 0001286  89.2500 270.8636 15.06406908246617"

    tle = TLE(line1,line2)

    print (tle)
    print ('Epoch :', tle.getDate())


    propagator = TLEPropagator.selectExtrapolator(tle)

    for i in range(0, 10):
        date = datetime.utcnow().replace(tzinfo=timezone.utc)
        now = date.astimezone()
        seconds = float(date.second) + float(date.microsecond) / 1e6

        absoluteDate = AbsoluteDate(date.year,
                                    date.month,
                                    date.day,
                                    date.hour,
                                    date.minute,
                                    seconds,
                                    TimeScalesFactory.getUTC())

        state = propagator.propagate(absoluteDate)


        pv = state.getPVCoordinates(itrf)

        p = pv.getPosition()
        v = pv.getVelocity()

        gp = earth.transform(p, itrf, absoluteDate)

        print("tm", now, "\n\tx:", p.getX(), "y:", p.getY(), "z:", p.getZ(), "\n\tvx:", v.getX(), "vy:", v.getY(), "vz:", v.getZ(),
            "\n\tlon:", np.rad2deg(gp.getLongitude()), "lat:", np.rad2deg(gp.getLatitude()), "alt:", gp.getAltitude())

        time.sleep(10)

计算结果与n2yo.com对比：
orekit-python-n2yo

卫星网络仿真应用

卫星网络如Starlink星座，卫星移动模型可以基于sgp4模型(Skyfield/orekit/…)进行位置预测，坐标系上自然需要使用考虑不受地球自传影响的ITRS坐标系，地面站通常需要WGS84经纬度转换为ITRS坐标系下固定坐标(相当于地面站移动模型是固定不动的)。

WGS84转ITRS坐标系可参考wgs84-geocentric进行转换，如下为北京地区转换：
wgs84-geocentric

卫星网络中卫星及地面站通过ITRS坐标系建立3D网络拓扑图，根据连接关系和位置距离可初步对星地链路建模(如波速分配覆盖/连接可用性/带宽/延迟等)和星间链路建模(如带宽/延迟/丢包率等)，根据全局动态拓扑时间切片后拓扑也可该时间切片下最短路径等计算，并根据路径配置路由进行网络仿真。

sat-net-sim

基于Skyfield&Orekit中SGP4轨道位置计算及卫星网络仿真中应用

https://xiaohuiluo.github.io/2025/01/23/skyfield-and-orekit/

作者

HUI

发布于

2025-01-23

更新于

2025-01-23

许可协议

#

评论

关注我

分类

最新文章

归档

标签

H31
linux1
lpm1
mptcp1
multipath1
netcfg1
onos3
orekit1
sdn3
sgp41
skyfield1
subsystem1
tle - 两行轨道数据1
多路径1
最长匹配1

订阅更新

×