python pyephem库 卫星轨道计算及两行轨道参数TLE

python pyephem库 卫星轨道计算及两行轨道参数TLE

作者微信公众号:小卫星

 操作系统: Ubuntu18.04 LTS

IDE: PyCharm Community 2018.1 

 Python 2.7

1、安装

$ pip install pyephem

http://rhodesmill.org/pyephem/

2、测试

先找颗卫星,事实上可以用python自动弄下来,进北美防空司令部的:

http://celestrak.com/NORAD/elements/beidou.txt

可以看到北斗卫星的两行轨道参数(TLE),把前三行复制进来,自己的经纬度可以在google上找:

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import ephem
#http://rhodesmill.org/pyephem/quick.html
me = ephem.Observer()
me.lon, me.lat, me.elevation = '108.5000', '34.5000', 800.0
line1 = "BEIDOU G1"
line2 = "1 36287U 10001A   18194.84994825 -.00000297  00000-0  00000-0 0  9999"
line3 = "2 36287   1.4745   3.9768 0004214 156.8560 216.8921  1.00270561 31144"
sat = ephem.readtle(line1, line2, line3)

me.date = '2018/7/14'#ephem.now()
sat.compute(me)

print(sat.az * 180.0 / 3.1416) #卫星的方位角
print(sat.alt * 180.0 / 3.1416 ) #卫星的仰角
print(sat.range_velocity) #卫星相对于观察者的运动速率,为正,表示卫星正在远离观察者。

输出为:

/home/wy/PycharmProjects/wP2/venv/bin/python /home/wy/.PyCharmCE2018.1/config/scratches/wyEPH1.py
131.486721673
38.9986189838
-4.07362747192

Process finished with exit code 0

表示卫星仰角为39度,方位角为131.5度。这是个GEO卫星,所以改变时间,结果变化不大。但是因为星历有时效性,时间设的和星历时间太远,则程序计算误差极大(部分原因是卫星本身会进行轨道调整)。再设的过大,程序不会计算。

如2018年7月14日download的星历,设时间为‘2008/7/10’,会报错:

me.date = '2008/7/14'#ephem.now()
/home/wy/PycharmProjects/wP2/venv/bin/python /home/wy/.PyCharmCE2018.1/config/scratches/wyEPH1.py
Traceback (most recent call last):
  File "/home/wy/.PyCharmCE2018.1/config/scratches/wyEPH1.py", line 13, in <module>
    sat.compute(me)
ValueError: TLE elements are valid for a few weeks around their epoch, but you are asking about a date 3651 days from the epoch

Process finished with exit code 1

3、卫星轨道参数来源

卫星轨道来源其实有很多,主要有两个:

(1)彗星、小行星等自然天体轨道

International Astronomical Union (IAU)的Minor Planet Center (MPC),以接收和分发自然天体轨道为目的,网站:

https://www.minorplanetcenter.net/iau/Ephemerides/Soft03.html

例如:https://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft03Cmt.txt

# From MPC106342
C/1995 O1 (Hale-Bopp),e,89.0947,283.1330,130.8144,183.1685,0.0003976,0.99490338,0.0000,03/29.3733/1997,2000,g -2.0,4.0

其格式为xephem格式,需要用readdb函数读取:

ephwy = ephem.readdb("C/1995 O1 (Hale-Bopp),e,89.0947," + \
                     "283.1330,130.8144,183.1685,0.0003976,0.99490338,0.0000," + \
                     "03/29.3733/1997,2000,g -2.0,4.0")
ephwy.compute('2018/7/14')
print(ephwy.name)
print("%s %s" % (ephwy.ra, ephwy.dec))
print("%s %s" % (ephem.constellation(ephwy), ephwy.mag))

输出为:

C/1995 O1 (Hale-Bopp)
0:30:39.87 -85:43:48.2
('Oct', 'Octans') 22.21


(2)人造卫星、船、器

两行轨道数据(TLE,Two-Line Orbital Element),由美国celestrak发明创立,是用于描述太空飞行体位置和速度的表达式。

http://celestrak.com/NORAD/elements/

例如:

http://celestrak.com/NORAD/elements/beidou.txt

最后一行(2018.07.14基线)

BEIDOU-3 M8             
1 43246U 18029B   18194.72004803 -.00000064  00000-0  10000-3 0  9991
2 43246  55.0642  39.6410 0004409 332.7876  76.6516  1.86232943  1964
第0行,将第1行视为0行,是 卫星通用名称,最长为24个字符。
第1行和第2行是标准的卫星星历格式(TLE格式),每行69个字符,包括0~9,A~Z(大写)、空格、点和正/负号。

卫星星历编号含义:

第1行:

(1)43246U。43246是北美防空司令部(NORAD)自己给的卫星编号。U代表不保密。C、S都是NORAD自己才能看到,是保密的。 

(2)18029B。国际编号,18表示2018年发射,029表示是这一年的第29次发射,B则表示是一箭多星,依次有A、B、C....  。北斗的MEO卫星都是一箭双星,所以只有A、B。

(3)18194.72004803。表示这组轨道数据的时间点,18是2018年,194表示第194天,也就是7月14日。72004803表示这一天里的时刻。

(4).00069181  13771-5  44016-2。为轨道模型参数。

(5)0。为轨道模型类型,SGP4/SDP4轨道模型。

(6)999。表示是该空间飞行器的第999组TLE参数,我错略把各种卫星看了下,发现都是999。

(7)7。为校验位。

第2行:

(1)43246。NORAD卫星编号。

(2)55.0642。轨道倾角。

(3)39.6410。升交点赤经。

(4)0004409。轨道偏心率,实际值为0.0004409。

(7)332.7876。近地点幅角。

(8)76.6516。平近点角。表示这组TLE对应的时刻时,卫星在轨道的什么位置。

(9)1.86232943。每天环绕地球的圈数,其倒数就是周期。

(10)196。发射以来飞行的圈数,可以推算,196/1.86232943=105.24454天就是多少时间之前发射的。

查下新闻就知道:2018年3月30日1时56分,中国在西昌卫星发射中心用长征三号乙运载火箭(及远征一号上面级),以“一箭双星”方式成功发射第三十、三十一颗北斗导航卫星。这两颗卫星属于中圆地球轨道卫星,是中国北斗三号第七、八颗组网卫星。是吻合的。

(11)4 。校验位。


这个东西强大到发射后几天就有数据,而且精确度极高,但是使用时务必注意其中的时间都是UTC时间,中国用的是UTC+8,因此时间加8小时才能对应上。

要在osgEarth 3.2中根据TLE参数绘制卫星轨道,你可以使用以下代码片段: ```cpp #include <iostream> #include <fstream> #include <string> #include <sstream> #include <iomanip> #include <cmath> #include <osgEarth/MapNode> #include <osgEarthUtil/Controls> #include <osgEarthUtil/EarthManipulator> #include <osgEarthUtil/Sky> #include <osgEarthUtil/ExampleResources> #include <osgEarthAnnotation/FeatureNode> #include <osgEarthSymbology/Style> #include <osgEarthDrivers/tms/TMSOptions> #include <osgEarthDrivers/wms/WMSOptions> #include <osgEarthDrivers/feature_ogr/OGRFeatureOptions> #include <osgEarthFeatures/Feature> #include <osgEarthFeatures/FeatureModelLayer> #include <osgEarthUtil/Controls> #include <osgEarthUtil/ActivityMonitorTool> #include <osgEarthUtil/RTTPicker> #include <osgEarthUtil/MouseCoordsTool> #include <osgEarthUtil/EarthManipulator> #include <osgEarthUtil/ViewFitter> #include <osgEarthUtil/ExampleResources> using namespace osgEarth; using namespace osgEarth::Features; using namespace osgEarth::Symbology; using namespace osgEarth::Util; using namespace osgEarth::Annotation; // 用于将字符串转换为浮点数的函数 double to_double(const std::string& str) { std::istringstream ss(str); double value; ss >> value; return value; } // 用于将时间字符串转换为时间戳的函数 time_t to_time_t(const std::string& time_str) { std::tm timeinfo = {}; std::istringstream ss(time_str); ss >> std::get_time(&timeinfo, "%Y-%m-%d %H:%M:%S"); return std::mktime(&timeinfo); } // 计算卫星在给定时间的位置 osg::Vec3d calculate_satellite_position(double time_since_epoch, double inclination, double ascending_node_longitude, double eccentricity, double argument_of_perigee, double mean_anomaly, double semi_major_axis) { // 计算卫星轨道上的平均角速度 double mean_motion = 2.0 * osg::PI / (86400.0 / mean_anomaly); // 计算卫星在给定时间的平均角度 double mean_angle = mean_anomaly + mean_motion * time_since_epoch; // 计算卫星在给定时间的偏心率角 double eccentric_anomaly = mean_angle; for (int i = 0; i < 10; i++) { double delta = (eccentric_anomaly - eccentricity * sin(eccentric_anomaly) - mean_angle) / (1.0 - eccentricity * cos(eccentric_anomaly)); eccentric_anomaly -= delta; if (fabs(delta) < 1e-12) { break; } } // 计算卫星在给定时间的真近点角 double true_anomaly = 2.0 * atan(sqrt((1.0 + eccentricity) / (1.0 - eccentricity)) * tan(0.5 * eccentric_anomaly)); // 计算卫星在给定时间的升交点经度 double raan = ascending_node_longitude * osg::PI / 180.0; double arg_perigee = argument_of_perigee * osg::PI / 180.0; double inclination_rad = inclination * osg::PI / 180.0; double x = semi_major_axis * (cos(raan) * cos(true_anomaly + arg_perigee) - sin(raan) * sin(true_anomaly + arg_perigee) * cos(inclination_rad)); double y = semi_major_axis * (sin(raan) * cos(true_anomaly + arg_perigee) + cos(raan) * sin(true_anomaly + arg_perigee) * cos(inclination_rad)); double z = semi_major_axis * sin(true_anomaly + arg_perigee) * sin(inclination_rad); return osg::Vec3d(x, y, z); } int main(int argc, char** argv) { // 读取TLE文件 std::ifstream tle_file("tle.txt"); if (!tle_file.is_open()) { std::cerr << "Failed to open TLE file!" << std::endl; return -1; } // 读取第一行TLE参数 std::string line1; std::getline(tle_file, line1); // 读取第二行TLE参数 std::string line2; std::getline(tle_file, line2); // 解析TLE参数 std::string name = line1.substr(0, 24); double inclination = to_double(line2.substr(8, 8)); double ascending_node_longitude = to_double(line2.substr(17, 8)); double eccentricity = to_double("0." + line2.substr(26, 7)); double argument_of_perigee = to_double(line2.substr(34, 8)); double mean_anomaly = to_double(line2.substr(43, 8)); double mean_motion = to_double(line2.substr(52, 11)); double semi_major_axis = pow(398600.4418 / pow(mean_motion, 2), 1.0 / 3.0) * 1000.0; // 创建地球模型 osgViewer::Viewer viewer; osg::ref_ptr<osgEarth::MapNode> mapNode = osgEarth::Util::MapNodeHelper().load(arguments, &viewer); viewer.setSceneData(mapNode); // 计算卫星在当前时间的位置 double time_since_epoch = difftime(time(0), to_time_t("2022-01-01 00:00:00")); osg::Vec3d satellite_position = calculate_satellite_position(time_since_epoch, inclination, ascending_node_longitude, eccentricity, argument_of_perigee, mean_anomaly, semi_major_axis); // 创建轨道线段的几何体 osg::ref_ptr<osg::Vec3Array> vertices = new osg::Vec3Array; vertices->push_back(osg::Vec3d(0, 0, 0)); vertices->push_back(satellite_position); osg::ref_ptr<osg::Geometry> geometry = new osg::Geometry; geometry->setVertexArray(vertices.get()); geometry->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, 0, 2)); // 创建轨道线段的节点 osg::ref_ptr<Feature> feature = new Feature(geometry.get(), mapNode->getMapSRS(), Feature::UNDEFINED); osg::ref_ptr<FeatureNode> featureNode = new FeatureNode(mapNode, feature.get()); // 添加轨道线段节点到地球上 mapNode->addChild(featureNode.get()); viewer.run(); return 0; } ``` 这段代码使用 osgEarth 创建了一个地球模型,并使用TLE参数计算卫星在当前时间的位置,并在地球上绘制了一条卫星轨道线。你可以根据自己的需要修改TLE文件路径,调整轨道的位置和方向。
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值