卫星星历外推实战：TLE、SGP4、导航电文与切比雪夫拟合

简介：《卫星星历推算工具》面向卫星通信、导航定位、天文观测及空间科学教学中的轨道预测需求，提供可导入公开星历、结合动力学模型外推卫星位置的工具环境，可计算特定时刻的高度、经纬度及相对地面站的方位角与俯仰角，服务于天线指向、信号接收与卫星过境安排，也可辅助完成过境预报和地面站跟踪准备。压缩包共317个文件，以txt星历与说明、lng多语言、htm帮助文档、bmp界面素材、dat/owm配置数据及exe主程序为主，整体约1.55MB，结构紧凑，便于在普通桌面环境部署试用。目前已有1368人学习。配套的Orbitron支持多格式星历解析、实时位置显示、未来轨道预测和二维/三维轨道图绘制，读者可据此理解轨道动力学参数，完成通信链路预判、观测计划编排与课堂演示，兼顾入门操作与进阶研究的实际需求，适合作为轨道动力学与卫星应用的辅助材料。

1. 卫星星历推算：从两行根数外推 72 小时，位置为什么能飘出几十公里

把一条两行根数丢进脚本，算出 3 天后某颗低轨卫星的位置——听起来十几行代码就够，真跑出来常发现算出的点跟实测轨迹差了十几甚至几十公里。卫星轨道的星历外推从来不是套公式，它是一串前置约定：手里的根数是平均根数还是瞬时根数、时间系统走的是 UTC 还是系统时、坐标系停在惯性系还是转到了地固系、外推时长有没有越过模型本身的适用边界。这份卫星星历推算工具要解决的正是这条链路——把 TLE、导航电文、SP3 采样点这三类输入统一到“给定时刻 → 位置速度”的输出上，并把外推精度衰减这件事显式暴露出来。做地面站过境预报、遥感成像几何、载荷仿真、轨道可视化的人都能直接用；只想要个漂亮三维轨迹的人，也能用，但得先接受它的误差条。

2. 三条输入路线选型：TLE、导航电文、SP3 采样点差在哪

2.1 平均根数、瞬时根数、采样点：三种“星历”根本不是一回事

很多人把 TLE 当成卫星的实时轨道参数，这是第一层误解。TLE 里那六个元素是平均根数——把短周期摄动（地球非球形引力的周期项、日月引力的短周期项）扣掉之后剩下的一套平滑元素，它只对配套的解析模型有意义。SGP4 干的活就是把这些平均根数再“加回去”，用一套固定的大气阻力模型和引力场模型重建位置。所以你不能把 TLE 的倾角直接拿去做几何计算，也不能把它丢进一个纯开普勒递推器里，那样出来的结果连方向都可能是错的。

导航电文里的参数是另一个物种：它是瞬时开普勒根数加摄动改正项，参考历元明确到秒，配的是“先算无摄动开普勒轨道，再叠加谐波改正”这套解法。它的元素在参考历元附近是真实的瞬时值，不是平均值，所以外推几小时还算靠谱，超过半天误差就会开始爬。

SP3 更直接——它就是等间隔采样出来的三维坐标序列，没有元素、没有模型，只有点。你对它做的是插值或拟合，不是外推。这一点很关键：很多人拿 SP3 去做长时间外推，本质上是在用多项式往外“猜”，精度衰减比模型外推还快。

我的选型习惯是先问一句：你要的时间点，落在参考历元之后多久？这个数字基本决定了你该走哪条路。

2.2 精度天花板对照表：先知道能到哪，再决定投多少工

下面这张表是我自己在做方案评估时用的，数字是量级，具体值跟轨道高度、太阳活动水平、外推时长强相关，别当承诺值用。

输入形式	元素性质	配套模型	参考历元附近	外推可用跨度	原生坐标系
TLE 两行根数	平均根数	SGP4 / SDP4	公里级	1～3 天	TEME
导航电文	瞬时根数 + 摄动改正	开普勒递推	米级到十米级	数小时	地固系
SP3 采样点	离散坐标	切比雪夫/拉格朗日插值	分米级	受采样间隔约束	地固系或惯性系

太阳活动高年，大气密度模型误差会被放大，同一颗低轨卫星的外推误差可能是平静期的两三倍。这一点在表格里体现不出来，但做任务规划时必须留余量。我一般会在 24 小时外推的误差预算里直接按表格量级乘 1.5 到 2。

2.3 外推时长是唯一硬指标，其它都是次要项

把选型逻辑压成一句话：参考历元之后 2 小时以内，三条路线都能用，选你手上现成的那条；2 小时到 2 天，用 TLE + SGP4；超过 2 天，别指望单模型，得靠多源数据融合或者干脆重新取根数。

这里有个反直觉的点：外推 1 天和 3 天，误差不是线性增长，而是接近三次方到指数式增长。原因是解析模型里的大气阻力项对时间高度敏感，密度模型本身有 10%～30% 的相对误差，这个误差在外推过程中被积分两次，直接进了半长轴。所以“外推 3 天只比 1 天差 3 倍”这种估计方法一定会翻车。

另一个常被忽略的因素是轨道类型。低轨（高度几百公里）受大气阻力主导，外推衰减最快；中高轨主要受日月引力和地球非球形引力高阶项影响，衰减慢得多，同样外推 3 天，误差可能只有低轨的十分之一。所以选型时先看轨道高度，再看外推时长，两个一起决定。

2.4 先把时间和坐标系的约定立死，后面才不会互相甩锅

我在项目里见过太多次“结果对不上”，最后发现是时间系统混了。UTC、TAI、TT、UT1 之间差的是几十秒量级，对地固系位置来说，几十秒的地球自转就是上百米到几公里的位移。工具内部我一般统一成两个量：连续儒略日（用于动力学递推）和系统时秒计数（用于导航电文）。

坐标系同理，必须写死。SGP4 出来的是 TEME（真赤道平春分点），导航电文递推出来的是地固系，SP3 通常也是地固系。这三者之间要转动，而转动矩阵依赖时间。工具的规矩是：内部计算一律在各自的“原生系”里算完，出口处统一转到一个指定的输出系，绝不在计算中途来回转。中途转一次就多一次时间对齐的机会，也就多一个错误来源。

提示：这两条约定要写进工具的接口文档里，而不是留在代码注释里。接口签名上直接标明输入时间系统、输出坐标系，后面接手的人能少踩一半坑。

3. 把 TLE 喂进 SGP4：解析、调用、批量落盘的一段可复现代码

3.1 两行半的字段拆解与校验

TLE 的格式是定长字段，靠列位置切，不靠空格切。这一点很多人第一次写解析器时会栽——卫星名字里带空格、负号占位、指数格式（如30777-3表示 0.30777e-3）都会让split()失效。

关键字段位置（第 1 行）：轨道编号在第 3～7 列，国际编号 9～17 列，历元在 19～32 列（YYDDD.DDDDDDDD），平均运动一阶导数 34～43 列，二阶导数 45～52 列，BSTAR 阻力项 54～61 列。

第 2 行：轨道倾角 9～16 列，升交点赤经 18～25 列，偏心率 27～33 列（隐含小数点，前面补0.），近地点幅角 35～42 列，平近点角 44～51 列，平均运动 53～63 列。

最后一位是模 10 校验和：把该行前 68 个字符里的数字累加，负号算 1，其它字符算 0，对 10 取模。校验和不对的 TLE 我建议直接拒绝，因为错一位数字带来的可能是整条轨道跑偏。

def parse_tle(line2: str) -> dict: """按列位置解析 TLE 第 2 行，返回开普勒根数（角度制）""" return { "incl": float(line2[8:16]), # 轨道倾角 deg "raan": float(line2[17:25]), # 升交点赤经 deg "ecc": float("0." + line2[26:33].strip()), # 偏心率，隐含小数点 "argp": float(line2[34:42]), # 近地点幅角 deg "ma": float(line2[43:51]), # 平近点角 deg "n": float(line2[52:63]), # 平均运动 rev/day }

[8:16]这类切片用的是 Python 的零基索引，跟 TLE 规范里的“第 9 列”差 1，写的时候要心里有数，别两边混着用。偏心率那行必须用字符串拼接的方式补小数点，直接float(line2[26:33])会拿到一个错误的整数。

3.2 调用 SGP4 拿到 TEME 下的位置速度

解析完不用自己实现 SGP4——这套模型有大量分支（近地点高度低于 220 km 时要切到另一种阻力处理，深空目标要切到 SDP4 的分支），手写基本等于自找麻烦。常用做法是调成熟库，把whichconst显式指定成 WGS72，跟 TLE 的生成基准保持一致。

from sgp4.api import Satrec, WGS72, jday import numpy as np # 下面两行是示意数据，实际使用时替换成你自己的 TLE LINE1 = "1 99999U 24001A 24001.50000000 .00000000 00000-0 10000-3 0 0001" LINE2 = "2 99999 97.5000 120.0000 0010000 90.0000 270.0000 15.20000000 00001" sat = Satrec.twoline2rv(LINE1, LINE2, WGS72) # 目标时刻：2024-01-03 00:00:00 UTC jd, fr = jday(2024, 1, 3, 0, 0, 0.0) err, r_teme, v_teme = sat.sgp4(jd, fr) if err != 0: raise RuntimeError(f"SGP4 返回错误码 {err}") print("TEME 位置 km:", r_teme) # [x, y, z] print("TEME 速度 km/s:", v_teme)

jday返回的是整日儒略日和日内小数两部分，分开传是为了避免大数相减丢精度——儒略日已经到 246 万量级，双精度浮点在这个量级上的有效位只剩 10 微秒左右，合并成一个数会把精度吃掉一部分。这一点在长时间序列批量计算时影响明显。

err非零要当成硬错误处理，常见的是 1（均值偏心率越界）、2（均值运动越界）、6（卫星已再入）。库不会抛异常，只给码，不判就会静默输出垃圾坐标。

3.3 从 TEME 转到地固系：一个够用的 GMST 旋转

TEME 是准惯性系，做地面站可见性、成像条带计算都要转到地固系。严格做法要走岁差、章动、极移一整串，但对公里级精度的任务，用格林尼治平恒星时（GMST）绕 Z 轴转一次就够了，误差在几十米量级，比 SGP4 本身的误差小一个数量级以上。

def teme_to_ecef(r_teme, jd_ut1): """TEME -> 地固系，仅做 GMST 绕 Z 轴旋转，适合公里级精度场景""" d = jd_ut1 - 2451545.0 T = d / 36525.0 gmst_deg = (280.46061837 + 360.98564736629 * d + 0.000387933 * T**2 - T**3 / 38710000.0) % 360.0 th = np.deg2rad(gmst_deg) Rz = np.array([[ np.cos(th), np.sin(th), 0.0], [-np.sin(th), np.cos(th), 0.0], [0.0, 0.0, 1.0]]) return Rz @ np.asarray(r_teme)

jd_ut1这里应该用 UT1，手上只有 UTC 时直接拿来用，差值在 0.9 秒以内，对应地表约 400 米的位移。对过境预报够用，对精密定轨不够，这一点要在接口里标清楚，别让下游误以为是精密结果。

3.4 批量外推与结果写入：一次算完，别在循环里反复建对象

批量场景（比如算一颗星 3 天内的每分钟位置，或者一个星座 24 小时的全部过境）有两个性能点值得注意：Satrec对象只建一次，循环里只调sgp4；时间序列用numpy一次性生成，别在 Python 循环里做浮点累加。

import csv from datetime import datetime, timedelta def batch_propagate(sat, t0: datetime, minutes: int, step_s: int, out_csv: str): rows = [] for k in range(0, minutes * 60 + 1, step_s): t = t0 + timedelta(seconds=k) jd, fr = jday(t.year, t.month, t.day, t.hour, t.minute, t.second + t.microsecond * 1e-6) err, r, v = sat.sgp4(jd, fr) if err != 0: continue ecef = teme_to_ecef(r, jd + fr) rows.append([t.isoformat(), *r, *v, *ecef]) with open(out_csv, "w", newline="") as f: w = csv.writer(f) w.writerow(["utc", "x_teme", "y_teme", "z_teme", "vx", "vy", "vz", "x_ecef", "y_ecef", "z_ecef"]) w.writerows(rows) return len(rows)

step_s的选择看用途：画轨迹 60 秒足够，做过境仰角曲线 10～30 秒，做遮挡分析要压到 1 秒以下。步长越小不是越准，SGP4 是解析模型，任意时刻独立求值，步长只影响曲线平滑度，不影响单点精度——这一点跟数值积分器完全不同，别拿数值积分的直觉来调这里。

注意：批量输出一定带上 UTC 时间戳和坐标系标记。我见过导出文件只留三列坐标、事后没人知道是哪个系的，整批数据作废重算。

4. 导航电文外推：开普勒递推加摄动改正的实操拆解

4.1 导航文件里到底要取哪些参数

RINEX 导航文件（版本 3 之后按星座分文件）每个历元一块，每块 8 行左右，取的是参考历元 toe、开普勒六根数（sqrtA、e、i0、Ω0、ω、M0）、以及六个谐波改正系数（Cuc、Cus、Crc、Crs、Cic、Cis）、两个长期项（Δn、IDOT），再加上一个星期数用于时间归算。

这些参数里最容易取错的是sqrtA——文件里给的是半长轴的平方根，单位是 $\sqrt{m}$，必须平方回去才是半长轴。直接当半长轴用会让轨道尺度差出三个数量级，位置直接飞出屏幕。第二个坑是toe是周内秒，跟绝对时间差一个周跳，跨周计算时要做 ±604800 秒的修正。

4.2 偏近点角迭代与地固系位置解算

开普勒方程 $M = E - e\sin E$ 对偏心率很小的导航卫星收敛很快，直接用不动点迭代，二三十次内必然到 1e-13 量级。收敛之后再算真近点角、升交角距，最后叠加上六个谐波改正。

import numpy as np MU = 3.986005e14 # 地球引力常数 m^3/s^2 OMEGA_E = 7.2921151467e-5 # 地球自转角速度 rad/s def broadcast_ecef(p: dict, t_sow: float, week_roll: float = 0.0) -> np.ndarray: """p 为导航电文字典，t_sow 为目标时刻周内秒，返回地固系位置(m)""" A = p["sqrtA"] ** 2 n0 = np.sqrt(MU / A ** 3) n = n0 + p["dn"] tk = t_sow - p["toe"] + week_roll if tk > 302400: tk -= 604800 if tk < -302400: tk += 604800 Mk = p["M0"] + n * tk Ek = Mk for _ in range(30): Ek_new = Mk + p["e"] * np.sin(Ek) if abs(Ek_new - Ek) < 1e-13: Ek = Ek_new break Ek = Ek_new sinE, cosE = np.sin(Ek), np.cos(Ek) vk = np.arctan2(np.sqrt(1 - p["e"]**2) * sinE, cosE - p["e"]) phi = vk + p["omega"] du = p["Cus"] * np.sin(2*phi) + p["Cuc"] * np.cos(2*phi) dr = p["Crs"] * np.sin(2*phi) + p["Crc"] * np.cos(2*phi) di = p["Cis"] * np.sin(2*phi) + p["Cic"] * np.cos(2*phi) u = phi + du r = A * (1 - p["e"] * cosE) + dr i = p["i0"] + di + p["IDOT"] * tk x_orb, y_orb = r * np.cos(u), r * np.sin(u) Om = p["OMEGA0"] + (p["OMEGADOT"] - OMEGA_E) * tk - OMEGA_E * p["toe"] x = x_orb * np.cos(Om) - y_orb * np.cos(i) * np.sin(Om) y = x_orb * np.sin(Om) + y_orb * np.cos(i) * np.cos(Om) z = y_orb * np.sin(i) return np.array([x, y, z])

week_roll那个参数是给跨周场景留的口子。tk的 ±302400 秒判别是导航电文算法里明确规定的，不是可选优化——不做这一步，跨半周的目标时刻会算出完全错误的位置。

Om那一行里的- OMEGA_E * toe很多人会漏。它来自“地固系随地球转动”的坐标变换，漏掉之后整条轨道会绕 Z 轴整体偏转一个跟 toe 成正比的角度，而且这个角度每周还不一样，看起来像随机误差。

4.3 外推半天以后，误差是怎么爬上去的

导航电文的参数是以 2 小时为周期播发的，参考历元通常取整 2 小时。外推到历元后 2 小时，误差还在米级到十米级；到 4 小时就开始明显爬；超过 6 小时，基本只能当粗定位用。

误差来源有三块：一是谐波改正项本身只拟合到二阶，高阶项被丢掉了；二是OMEGADOT、IDOT这类长期项是线性近似，时间一长就偏离；三是参数本身是广播值，不是精密值。三块叠在一起，就是那条误差曲线。

实测过的做法是给每个外推结果配一个“年龄”字段（目标时刻减参考历元），下游按年龄分级使用：2 小时内正常用，2～4 小时降权，4 小时以上只用于粗筛。这个字段比误差估计值更实用，因为它不需要你预先知道误差模型。

另外一个细节是信号传播的时间改正。如果算的是“接收机在某时刻看到卫星在哪”，要解算信号发射时刻，做一个 60～80 毫秒的迭代；如果要算地固系位置，还要在卫星位置里叠加上地球自转引起的 Sagnac 改正。做几何仿真时这两步漏一步，误差就在几十米量级，跟轨道误差同量级，很容易被误判成“模型不准”。

提示：导航电文里toe是周内秒，t_sow也是。任何一方用绝对时间，另一方用周内秒，都会得到差出一个星期的结果。接口里把两者的单位写死在参数名上，比如toe_sow、t_sow。

5. 避坑与排查：五条把结果算飞的血泪记录

5.1 现象：轨迹整体平移几公里，形状却完全正确

原因：时间系统混用。TLE 的历元字段是 UTC，SGP4 内部按 UT1 处理，导航电文按系统时，SP3 头文件里明确写了是哪种时标。从不同来源取的数据拼在一起，没做统一，几十秒的偏差就变成几公里的地表位移。

解决：在工具入口做一次强制归一，所有时间进来先转成“连续儒略日 + 时标标记”这一个结构体，出口再按需转出去。中间任何计算都不允许直接吃裸的 datetime。这条规矩立下来之后，这类问题再没出现过。

5.2 现象：算出来的星下点经纬度整体偏了一个固定角度

原因：坐标系没转。SGP4 出来的是 TEME，直接拿去做经纬度反算，等于把惯性系坐标当成了地固系。误差量级是地球自转角乘以时间，一天之内的目标能偏出一整圈。

解决：在函数签名上强制区分r_teme和r_ecef，名字不同、类型不同（用 numpy 数组 + 一个显式的 dataclass 包装），传错直接运行时报错。靠命名约定防不住，靠类型能防住。

5.3 现象：插值出来的轨迹在采样点之间抖动，残差比采样点本身还大

原因：插值阶数跟采样间隔不匹配。SP3 采样间隔 15 分钟，用 15 阶拉格朗日插值，两端会出现剧烈振荡（龙格现象），中间段也可能因为采样点太稀而放大噪声。

解决：改用分段切比雪夫拟合，每段 8～12 个采样点、10 阶左右，段与段之间留重叠区做平滑过渡。判定标准不看阶数，看拟合残差的 RMS，超过采样精度的 3 倍就降阶或者缩段长。

5.4 现象：误差在某个时刻突然跳变，前后都不大，就那一段离谱

原因：外推越界，触发了模型内部的保护分支。SGP4 在目标时刻距离历元过远时，内部会返回错误码或者退化结果；导航电文跨周计算没做周跳修正，tk突然变成七天的量级。

解决：批量计算前先做一遍“年龄检查”，超过设定阈值的时刻直接标记为不可用，不输出数值。宁可留空，也别输出一个看起来正常、实际已经失效的点——这种点在后续统计里最难查。

5.5 现象：同一套代码，换个数据源结果就全错，角度差了 57 倍

原因：角度制与弧度制混用。TLE 的倾角、升交点赤经是度，开普勒递推里的三角函数要弧度，两套数据走同一个函数时没人做转换。

解决：所有内部计算一律用弧度，只在解析和输出两个边界做转换，转换函数名里带deg或rad后缀。代码里禁止出现裸的np.sin(x)，必须能看到 x 是怎么来的。这条是血泪经验，改起来最便宜，漏起来最贵。

6. 长弧外推换条路：分段切比雪夫拟合加残差自检

当真需要覆盖几天的连续弧段，又不想每隔几小时重新取一次数据时，我的做法是不递推，改拟合。思路是：先用 SGP4 或导航电文在密集网格上生成一批“训练点”，再用分段切比雪夫多项式把这些点压成系数。运行时只查系数，既快又稳，而且误差可控——因为拟合残差当场就能算出来。

6.1 为什么拟合比递推更稳

递推的误差是累积的，每一步都带着上一步的偏差往前走；拟合的误差是局部的，只跟这一段里的采样点有关，段与段之间不传播。代价是你得先生成采样点，而采样点本身还是由递推来的——所以拟合不会提高精度，它只是把精度锁定在采样时刻的水平，并让中间的查询变得便宜且可复现。

另一个好处是数据量。3 天、1 秒步长的位置序列是 26 万个点，几十兆；压成每 6 小时一段、每段 10 阶的三维切比雪夫系数，只有几百个浮点数，嵌到嵌入式设备里都够。

import numpy as np from numpy.polynomial import chebyshev as Ch def fit_arc(ts: np.ndarray, xs: np.ndarray, deg: int = 10): """ts 为秒，xs 形状 (N,3)。返回 (t0, half, coefs)""" t0 = 0.5 * (ts.min() + ts.max()) half = 0.5 * (ts.max() - ts.min()) u = (ts - t0) / half # 归一化到 [-1, 1] coefs = np.array([Ch.chebfit(u, xs[:, k], deg) for k in range(3)]) return t0, half, coefs def eval_arc(t0, half, coefs, t): u = (np.asarray(t) - t0) / half return np.stack([Ch.chebval(u, coefs[k]) for k in range(3)], axis=-1) def arc_residual(ts, xs, t0, half, coefs): pred = eval_arc(t0, half, coefs, ts) return np.sqrt(np.mean((pred - xs) ** 2, axis=0)) # 每轴 RMS，单位 km

归一化那一步是必须的。切比雪夫多项式定义在 [-1, 1] 上，直接把秒数丢进去，高阶项的数值会炸到 1e20 量级，拟合彻底失效。half一定要跟着系数一起存下来，不然评估时没法还原。

6.2 段长和阶数怎么定：看残差，不看经验值

我一般从“段长 6 小时、10 阶”起步，然后看每段的 RMS 残差。低轨卫星在太阳活动高年，6 小时段可能残差就到几百米，这时候缩到 3 小时；中高轨用 12 小时段都绰绰有余。判定阈值我设成 SGP4 自身量级误差的 1/5——拟合不该成为误差主项。

轨道类型	起步段长	阶数	典型 RMS	调整方向
低轨（<600 km）	3 小时	10	百米级	残差大就缩段长
中轨	6 小时	10	十米级	一般不用调
高轨	12 小时	8	米级	可适当降阶

超过阈值时的处理顺序是：先缩段长，再升阶，最后才考虑换采样源。顺序反了会陷入“阶数越升越振”的循环。

6.3 交叉验证：用一半数据拟合，另一半做检验

最后一步是自检。把同一段弧的采样点按奇偶分成两组，一组拟合、一组验证，看验证组的残差是不是跟训练组同量级。如果验证残差明显更大，说明段长太长、模型在“背答案”，得缩段。

这套流程我现在每次都强制走一遍，起因是早年一次交付，拟合残差漂亮得不像话，结果下游拿去算过境，仰角曲线在段边界处有肉眼可见的折角——段与段之间没有重叠区，边界上系数不连续。从那以后我每次做分段拟合都强制留 20% 重叠、在重叠区做线性加权过渡，并且跑一遍奇偶交叉验证才敢往外发。希望帮到你。

本文还有配套的精品资源，点击获取