简介：面向卫星轨道仿真与STK对照验证需求的Python工程代码，围绕SGP4简化摄动模型与HPOP高精度轨道预报算法展开，适合航天专业学生、科研人员及编程爱好者用于轨道计算学习与算法验证。压缩包内共4个文件，包含一个可独立运行的Python仿真脚本、两张48小时轨道对比图与误差分析图、一份代码说明文档，整个RAR包约401KB，下载后即可按说明运行；目前已有402人学习下载。通过运行脚本可生成卫星轨道数据，并与STK可视化结果进行对比，直观观察不同方法在48小时内的轨道差异和误差变化；代码说明文档还详细解释了算法原理、输入参数定义、运行环境要求等内容，能够帮助读者理解卫星轨道建模的关键流程，快速复现仿真结果。同时，两张对比图像分别展示了轨道整体走势与逐点误差分布，便于从视觉和数值两个层面评估仿真精度，为课程设计、毕业设计或航天项目预研中的轨道仿真任务提供可直接参考的样例与优化依据。 搜“卫星轨道模型 python 仿真”的人，多半已经受够了网上那些零散代码片段——要么是拿现成库一行出图，要么是公式一大堆却没有能跑起来的完整实现。真正能把卫星轨道模型仿真跑通，同时每一步都清楚自己在算什么，才是这件事最核心的价值。这篇文章我会用纯 Python 从零写一套卫星轨道仿真代码，覆盖开普勒方程求解、轨道根数转位置速度、数值积分、J2 摄动和三维可视化。适合刚接触轨道力学、想让理论“落地成代码”的学生，也适合需要快速做轨道估计或任务预研的工程师。
1. 先搞清楚你在算什么：二体模型与轨道六要素
1.1 二体问题：为什么教科书总拿它开刀
卫星轨道仿真听起来高大上，但绝大多数入门场景其实都在解决同一个问题：给出一组轨道参数，算出某个时刻卫星在哪里、速度是多少。最基础的模型是“二体问题”——把地球和卫星都抽象成质点，只考虑两者之间的万有引力，忽略大气阻力、太阳光压、地球非球形引力等所有干扰。
二体问题的优雅之处在于它有解析解，卫星相对地球的运动轨迹被严格限定在一个平面内，形状是圆锥曲线。对环绕地球运行的航天器来说，轨道就是椭圆。这意味着我们不需要每一步都做复杂的数值积分，也能得到精确的位置，这在后面调试代码时会非常有用：先用解析解验证逻辑，再去搞数值积分。
生活里类比一下，二体问题就像“真空中的球形鸡”——虽然是极大的简化，但它把轨道运动最基本的骨架定义清楚了。没有这个骨架，后面加再多摄动项都是空中楼阁。
1.2 轨道六要素到底存的是什么信息
描述一条环绕地球的椭圆轨道，需要 6 个独立参数，也就是轨道六要素（Orbital Elements）：
- 半长轴 a：决定轨道大小，也通过开普勒第三定律决定轨道周期。
- 偏心率 e：决定轨道扁的程度，0 是正圆，0~1 之间是椭圆。
- 轨道倾角 i：轨道平面相对赤道平面的夹角，决定卫星能飞到多高的纬度。
- 升交点赤经 Ω：轨道平面在惯性空间里的朝向，即升交点相对春分点的经度。
- 近地点幅角 ω：轨道面内近地点相对升交点的角度。
- 平近点角 M：某个时刻卫星在轨道上的位置，通常取 0 表示过近地点。
这 6 个数合在一起，就能唯一确定一条轨道。前五个描述轨道“长什么样、朝哪个方向”，第六个描述“卫星现在走到哪了”。我最早写代码时总喜欢把 M 和“真近点角”混用，结果画出来的轨道位置永远对不上，后来才意识到这俩之间隔着一道开普勒方程。
1.3 解析法还是数值积分，先做选择
写代码前要想清楚：你要的是“某几个时刻的精确位置”，还是要“一段连续时间的轨道演化”。前者用解析法——直接由轨道根数算出任意时刻的位置速度，速度快、精度高，适合做轨道预报的初值；后者用数值积分——把加速度积分成速度、再积分成位置，适合研究摄动、轨道机动、编队飞行等复杂场景。
我建议初学者的路径是：先用解析法实现轨道根数到位置速度的转换，验证几个关键角度、坐标转换是否正确，然后再上 RK4 数值积分。顺序反了的话，你会被一堆叠加的误差搞得完全不知道是自己代码写错了，还是积分步长太大了。
2. 开普勒方程：轨道面上卫星位置怎么解
2.1 平近点角、偏近点角、真近点角
要算卫星在轨道面上的位置，核心是解开普勒方程。卫星沿椭圆轨道运动，速度不是均匀的——近地点快、远地点慢。为了统一描述这种变速运动，天体力学家引入三个角度：
- 平近点角 M：假想的均匀角速度，“平均”意义上的位置，随时间线性增加。
- 偏近点角 E：通过几何投影构造出来的辅助角度。
- 真近点角 ν：卫星实际相对近地点的真实角度，我们最终要的就是它。
三者之间通过开普勒方程连接：
M = E - e * sin(E)已知 M 和偏心率 e，求 E 是一个超越方程，没有初等解析解，必须用数值迭代。
2.2 迭代求解开普勒方程
工程上最常用的方法是牛顿迭代。定义 f(E)=E-e·sin(E)-M，迭代公式：
E_{n+1} = E_n - f(E_n) / f'(E_n) 其中 f'(E) = 1 - e·cos(E)写成 Python 函数特别简洁：
import numpy as np def kepler_solve(M, e, tol=1e-12): """牛顿迭代求解开普勒方程：M = E - e*sin(E)""" E = M if e < 0.8 else np.pi # 高偏心轨道给个更好的初值 for _ in range(100): f = E - e * np.sin(E) - M df = 1.0 - e * np.cos(E) dE = f / df E -= dE if abs(dE) < tol: break return E注意高偏心轨道的初值选择。若 e 接近 0.9 以上，直接从 M 开始迭代容易收敛慢甚至震荡，我习惯给 π 作为初值，实测几轮就能收敛。
得到偏近点角 E 后，真近点角 ν 可以通过半角公式计算：
E = kepler_solve(M, e) nu = 2.0 * np.arctan2( np.sqrt(1 + e) * np.sin(E / 2), np.sqrt(1 - e) * np.cos(E / 2) )这里用arctan2而不是arctan，否则象限会搞错。我吃过这个亏：用atan算出来角度总落在 -90°~90°，轨道前半段和后半段完全错乱。
2.3 最容易翻车的单位与象限
开普勒方程里的角和偏心率都必须是无量纲的。角度请一律用弧度，别在公式里混入角度制。Pyhton 里math.sin、np.sin默认都吃弧度，但很多人从文件里读出的是角度制轨道根数，忘了转弧度就传进函数，出来的位置错到离谱。
另一个坑是arctan2的参数顺序：是arctan2(y, x)，不是arctan2(x, y)。写半角公式时我一开始写成atan2(np.cos(...), np.sin(...))，结果真近点角永远差 90 度。这个错误非常隐蔽，因为轨道形状看起来是正常的，只是整体旋转了一个角度。
3. 从轨道面到地心惯性系：三次旋转别搞反
3.1 三个旋转角分别干什么
上一步算出的 r 还停留在“轨道平面坐标系”（PQW 系），坐标原点在地心，X 轴指向近地点，Z 轴指向轨道面法向。要让卫星位置变成地心惯性系（ECI）下的三维坐标，需要做三次旋转：
- 绕 Z 轴旋转
-ω：把近地点方向从 X 轴转出去，对准升交点方向。 - 绕 X 轴旋转
-i：把轨道平面“掰”到倾角 i 的位置。 - 绕 Z 轴旋转
-Ω：把升交点从春分点方向转到正确经度。
注意是-ω、-i、-Ω这个顺序，不能交换。旋转矩阵乘法的顺序就是坐标变换的顺序，搞反了轨道会在天上乱飞。
3.2 用旋转矩阵转出位置和速度
先定义两个基础旋转矩阵：
def rot_z(theta): c, s = np.cos(theta), np.sin(theta) return np.array([ [c, -s, 0], [s, c, 0], [0, 0, 1] ]) def rot_x(theta): c, s = np.cos(theta), np.sin(theta) return np.array([ [1, 0, 0], [0, c, -s], [0, s, c] ])然后通过轨道根数计算 PQW 系下的位置和速度。取 p=a(1-e²)，h=sqrt(μp)：
def coe2rv(a, e, i, Omega, omega, M): """轨道六要素 -> ECI 位置速度向量""" E = kepler_solve(M, e) nu = 2.0 * np.arctan2(np.sqrt(1 + e) * np.sin(E / 2), np.sqrt(1 - e) * np.cos(E / 2)) r_pqw = np.array([ a * (np.cos(E) - e), a * np.sqrt(1 - e**2) * np.sin(E), 0.0 ]) p = a * (1 - e**2) h = np.sqrt(mu * p) v_pqw = np.array([ -np.sqrt(mu / p) * np.sin(nu), np.sqrt(mu / p) * (e + np.cos(nu)), 0.0 ]) R = rot_z(-Omega) @ rot_x(-i) @ rot_z(-omega) return R @ r_pqw, R @ v_pqw这段代码值得反复看的重点是速度公式。速度不是把位置导一导就出来的，它由轨道力学推导而来：轨道面内速度的 X 分量和 Y 分量都跟真近点角 ν 直接相关。如果你想加深理解，可以试着从位置公式对时间求导，会发现结果正是这样。
3.3 对拍验证：拿真实轨道参数算一算
写完转换函数必须验证。拿一个简单的极轨道卫星，轨道根数设成 i=90°、Ω=0°、ω=0°，近地点取在 X 轴正方向。当 M=0 时，卫星应该在近地点，也就是 ECI 坐标系 X 轴正方向附近。跑一下代码：
mu = 3.986004418e14 a = 7000e3 # 半长轴 7000 km e = 0.01 i = np.deg2rad(90) Omega = np.deg2rad(0) omega = np.deg2rad(0) M = np.deg2rad(0) r, v = coe2rv(a, e, i, Omega, omega, M) print("r =", r) print("距离地心 =", np.linalg.norm(r))你会发现位置向量的 Z 分量接近 0（因为 M=0 时还在近地点附近），而轨道倾角 90° 决定了轨道面经过极地。距离约等于 a(1-e)，也就是近地点距离。这一步如果数值对不上，后面所有仿真都不要继续，先回头查旋转矩阵或角度单位。
3.4 一个小坑：坐标旋转的符号习惯
不同教材旋转矩阵的符号习惯不一样。有的定义绕 Z 轴转正角度是逆时针，有的默认顺指针，还有的用Rz(Ω) @ Rx(i) @ Rz(ω)而不是负数。关键不是死记公式，而是用“已知轨道面位置反推”的方式验证：比如 Ω=0、i=0、ω=0 时，轨道面就在赤道面，近地点方向就是 X 轴，旋转矩阵应该变成单位矩阵。这组退化条件能快速暴露符号问题。
4. 数值积分：让卫星真正“长”在引力场里跑
4.1 从根数到初始状态
解析法适合“问某个时刻卫星在哪”，但如果想观察轨道在 J2 摄动、大气阻力等影响下怎么慢变，就必须做数值积分。第一步还是用coe2rv生成初始时刻的位置和速度，然后把问题转化为初始值问题：
dr/dt = v dv/dt = a(r)其中加速度来自地球中心引力场：
a(r) = -mu * r / |r|^3这个形式极简，但包含了一个重要性质：加速度始终指向地心，大小和距离平方成反比。
4.2 RK4 积分器实现
数值积分方法很多，效率最高的入门方案是四阶龙格-库塔法（RK4）。它比欧拉法精度高得多，又不像高阶自适应方法那样复杂。代码实现如下：
def accel(r): """二体引力加速度""" r_norm = np.linalg.norm(r) return -mu * r / r_norm**3 def rk4_step(state, dt): """state = [x, y, z, vx, vy, vz]""" r = state[:3] v = state[3:] def f(s): r_, v_ = s[:3], s[3:] return np.concatenate([v_, accel(r_)]) k1 = f(state) k2 = f(state + 0.5 * dt * k1) k3 = f(state + 0.5 * dt * k2) k4 = f(state + dt * k3) return state + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)主循环里不断调用rk4_step即可。注意 state 的拼接顺序要和 f 函数里一致，否则速度被当成位置积分，轨道会在几秒内爆炸。
4.3 步长与能量守恒：仿真可信度的试金石
步长 dt 怎么选？工程经验是：低轨卫星（周期约 90 分钟）用 1~10 秒没问题；高轨（周期约 24 小时）可以放大到 60 秒甚至更大。但别盲目贪大，RK4 虽然精度高，步长太大会导致轨道严重漂移。
判断仿真是否可信，最好的办法是盯能量。二体问题中比机械能守恒：
epsilon = |v|^2/2 - mu/|r|不管轨道怎么走，这个值应该保持不变。我每次仿真都会顺手算一下：
epsilon = 0.5 * np.dot(v, v) - mu / np.linalg.norm(r)如果几步以后能量变化超过万分之几，说明步长太大或者代码有 bug。这个方法能帮你快速区分“算法不对”和“模型太简化”，调试体验完全不一样。
5. J2 摄动与三维可视化：仿真不只是画条椭圆
5.1 J2 摄动到底改了什么
真实地球不是完美球体，赤道部分隆起，导致对卫星的引力中多出一个高阶项，其中最主要的就是 J2 项。J2 摄动最直观的影响是让轨道面发生长期漂移：升交点赤经 Ω 和近地点幅角 ω 会随时间缓慢变化，这就是“轨道面进动”。对低轨卫星来说，这种效应非常显著，ISS 的轨道面每天都在东移。
物理上理解 J2，可以想象卫星在地球“腰带”的额外引力下，受到一个轻微的赤道面拉力。这个力让轨道法向方向慢慢转动。如果你做长期仿真完全忽略 J2，几十圈后轨道预报误差会是几百公里量级。
5.2 代码层面加 J2 有多简单
只需在加速度函数里加一个 J2 修正项。标准公式如下，RE 是地球赤道半径，r 是卫星到地心距离：
J2 = 1.08262668e-3 RE = 6378137.0 def accel_j2(r): r_norm = np.linalg.norm(r) x, y, z = r factor = 1.5 * J2 * mu * RE**2 / r_norm**5 zr2 = 5.0 * z**2 / r_norm**2 ax = -mu * x / r_norm**3 - factor * x * (1 - zr2) ay = -mu * y / r_norm**3 - factor * y * (1 - zr2) az = -mu * z / r_norm**3 - factor * z * (3 - zr2) return np.array([ax, ay, az])对比二体加速度，J2 项的量级很小，在低轨（约 500~800 km）大约是二体引力的千分之一，但它像“温水煮青蛙”一样持续作用，长期积分下轨道面进动效果非常明显。把accel换成accel_j2，其他代码一行都不用改，再跑几十圈轨道就能看到升交点位置在缓慢移动。
5.3 三维可视化与动态轨迹
最后一步是让仿真结果“看得见”。matplotlib 的三维投影足够用：
import matplotlib.pyplot as plt def plot_orbit(r_history): r_history = np.array(r_history) fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111, projection='3d') # 画地球示意 u = np.linspace(0, 2 * np.pi, 50) v = np.linspace(0, np.pi, 50) xs = RE * np.outer(np.cos(u), np.sin(v)) ys = RE * np.outer(np.sin(u), np.sin(v)) zs = RE * np.outer(np.ones_like(u), np.cos(v)) ax.plot_surface(xs, ys, zs, color='b', alpha=0.3) # 画轨道轨迹 ax.plot(r_history[:, 0], r_history[:, 1], r_history[:, 2], 'r-') ax.set_box_aspect([1, 1, 1]) plt.show()实际绘图时最重要的一行是set_box_aspect([1,1,1])，不设置的话三维坐标轴比例会自动拉伸，圆形轨道看起来像个大饼，很影响判断。
动态轨迹还要考虑内存问题：如果积分几千步，每一步都存 r 没问题，但如果做几万步的高频输出，建议按固定间隔采样，否则可视化会卡死。我一般先在积分循环里只保留位置、不保留速度，需要速度时再单独存，能省一半内存。
这轮写下来，我个人最大的体会是：轨道仿真的代码量其实不多，难的是每一步都得知道自己在算什么。开普勒方程解的是“轨道面上的几何位置”，旋转矩阵负责“把几何位置搬到惯性空间”，数值积分则是“让位置随时间演化”，三层逻辑各司其职。如果你照着上面的代码跑通一次，再把 J2 关掉对比轨道差异，基本上就算摸到轨道仿真的门槛了。最后分享一个小技巧：调试阶段把所有角度的中间值都打印出来，比如E、nu、omega+nu，用肉眼确认它们递增或变化趋势是否合理，这比对着坐标数值猜问题快得多。
本文还有配套的精品资源，点击获取
