简介这份资源聚焦地心惯性坐标系ECI与地心固定坐标系ECEF之间的转换面向从事卫星轨道计算、定位导航及航天器姿态分析的工程人员与相关专业学生。内容围绕地球自转对坐标的影响展开涉及儒略日期到格林尼治平恒星时的换算、恒星时与世界协调时的处理以及结合地球自转角速度求解ECEF坐标的完整思路并提及WGS84椭球模型在精度修正中的作用。压缩包共3个文件均为.m脚本体积约3KB分别承担时间参数换算与坐标转换等核心功能结构精简、便于直接嵌入MATLAB工程调用。目前已有925人学习下载适合需要快速实现ECI到ECEF转换、理解时间系统与坐标框架关系的读者参考也可作为轨道仿真与导航算法开发中的基础工具脚本使用。1. ECI 与 ECEF 坐标系转换从卫星轨道到地面站指向的落地链路做卫星轨道计算、地面站天线指向或者星上载荷几何校正的人迟早会撞上 ECI 和 ECEF 这两个坐标系。标题里那串ECITOECEF_ECI_ECI坐标系转换ECEF坐标系_eci坐标看着像随手敲的但它指向的问题非常具体给定一个在 ECI 惯性系下描述的位置或速度矢量怎么把它换算到随地球一起转的 ECEF 地固系里反过来又怎么回去。这件事在纸面上就是一个旋转矩阵但真正落地时会牵扯到时间系统、岁差章动、极移、地球自转角以及你用的到底是哪种 ECIJ2000 还是 TEME。新手容易以为套一个公式就完事熟手知道误差往往出在时间标签和参考系定义上。这篇笔记按“概念先立住、再动手复现、最后看坑”的顺序把这条转换链路拆开讲清楚适合做轨道仿真、遥感几何处理、地面站跟踪的工程师照着落地。2. ECI 与 ECEF 到底差在哪旋转矩阵背后的物理量2.1 两个坐标系的定义与选型理由ECIEarth-Centered Inertial是以地球质心为原点、坐标轴在惯性空间中基本不转的坐标系。常见的有 J2000 平赤道平春分点系和 TEMETrue Equator Mean Equinox系。J2000 的 X 轴指向 J2000.0 时刻的平春分点Z 轴指向该时刻的平北极适合做长期轨道积分和星历存储。TEME 则常用于 TLE 两行根数的传播结果它的春分点是“真”的和 J2000 之间差着岁差和章动。ECEFEarth-Centered Earth-Fixed以地球质心为原点X 轴指向本初子午线与赤道交点Z 轴指向协议地球极CTP随地球自转一起转。它适合描述地面站位置、地表目标点、卫星星下点轨迹。选型上如果你的输入是星历文件里的惯性系位置输出要给天线指向或地图叠加那 ECI 到 ECEF 就是必经之路。反过来如果你拿到的是 GPS 广播星历算出的 ECEF 位置想和惯性系下的摄动力模型对齐就要做 ECEF 到 ECI。两者不是简单转置关系因为中间夹着时间相关的旋转。2.2 转换的核心三个旋转矩阵的乘积从 ECI 到 ECEF常见做法是r_ECEF W(t) · R(t) · P(t) · N(t) · r_ECI其中 N 是章动矩阵P 是岁差矩阵R 是地球自转角矩阵W 是极移矩阵。如果用的是 J2000 到 ECEF这四项都要考虑如果精度要求不高或者时间跨度很短可以只保留 R把岁差章动极移忽略。很多工程代码里直接用一个 GMST 角构造 R 矩阵误差在角秒级对地面站指向够用对高精度定轨就不够。R 矩阵的形式是绕 Z 轴旋转地球自转角 θR [[cosθ, sinθ, 0], [-sinθ, cosθ, 0], [0, 0, 1]]注意方向从 ECI 到 ECEF 是顺时针转 θ所以矩阵里 sin 的符号和从 ECEF 到 ECI 相反。这个符号搞反是新手最常见的翻车点后面避坑章会细说。2.3 时间系统GMST 怎么算、用哪个时间θ 通常取格林尼治平恒星时 GMST。计算 GMST 需要 UT1 时间而日常拿到的是 UTC。UT1 和 UTC 差一个 DUT1量级在 0.9 秒以内。对大多数地面站指向直接用 UTC 代入 GMST 公式误差对应地表约 400 米天线波束宽的话可以接受对精密定轨必须用 UT1。一个常用的 GMST 近似公式IAU 1982GMST 280.46061837 360.98564736629 * (JD_UT1 - 2451545.0) (度)JD_UT1 是 UT1 的儒略日。算完取模 360。这个公式在 2000 年前后几十年内精度够用再往前或往后要考虑岁差模型更新。提示如果你只有 UTC又不想引入 DUT1 文件可以在代码里把 UTC 当 UT1 用但要在文档里写明这个近似别让下游以为你做了精密处理。3. 用 Python 跑通 ECI 到 ECEF 的最小实现3.1 环境准备与依赖选择我一般用 Python 做原型依赖三个库numpy做矩阵运算astropy做时间转换和坐标系定义sgp4或skyfield做 TLE 传播。如果只是验证旋转逻辑numpy加手写 GMST 就够。astropy的好处是它内置了 J2000、ITRS即 ECEF之间的完整转换链可以拿来当参考真值对比自己手写的结果。安装pip install numpy astropy如果你要处理 TLE再加pip install sgp43.2 手写旋转矩阵的完整代码下面这段代码实现从 J2000 ECI 到 ECEF 的简化转换只考虑地球自转忽略岁差章动极移。适合做教学验证和低精度场景。import numpy as np from datetime import datetime, timezone def julian_date(dt_utc): 把 UTC datetime 转成儒略日这里近似当 UT1 用 # 简化算法精度到秒级足够 y dt_utc.year m dt_utc.month d dt_utc.day frac (dt_utc.hour dt_utc.minute/60 dt_utc.second/3600) / 24.0 if m 2: y - 1 m 12 A y // 100 B 2 - A A // 4 jd int(365.25*(y4716)) int(30.6001*(m1)) d B - 1524.5 frac return jd def gmst_deg(jd_ut1): IAU 1982 GMST 公式返回度 T (jd_ut1 - 2451545.0) / 36525.0 gmst 280.46061837 360.98564736629 * (jd_ut1 - 2451545.0) return gmst % 360.0 def eci_to_ecef(r_eci, dt_utc): 输入 ECI 位置km输出 ECEF 位置km jd julian_date(dt_utc) theta np.radians(gmst_deg(jd)) # 绕 Z 轴顺时针旋转 theta R np.array([ [np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1] ]) return R np.asarray(r_eci, dtypefloat) # 测试假设某时刻 ECI 下卫星在 X 轴 7000 km dt datetime(2024, 1, 1, 0, 0, 0, tzinfotimezone.utc) r_eci [7000.0, 0.0, 0.0] r_ecef eci_to_ecef(r_eci, dt) print(ECEF:, r_ecef)逻辑说明julian_date把 UTC 转儒略日这里没有做 UT1 修正。gmst_deg用 IAU 1982 公式算格林尼治平恒星时。eci_to_ecef构造绕 Z 轴的旋转矩阵注意矩阵里 sin 的符号——从 ECI 到 ECEF 是顺时针所以第一行第二列是正 sin第二行第一列是负 sin。参数上r_eci单位是 km输出也是 kmdt_utc必须带时区否则 datetime 默认 naive算出来的儒略日会偏。3.3 用 astropy 做交叉验证手写容易错我习惯用astropy对一遍。下面代码把同一个 ECI 矢量转到 ITRS即 ECEF和上面结果对比。from astropy.time import Time from astropy.coordinates import CartesianRepresentation, GCRS, ITRS import astropy.units as u t Time(2024-01-01T00:00:00, scaleutc) # GCRS 接近 J2000 ECI gcrs GCRS(CartesianRepresentation(7000*u.km, 0*u.km, 0*u.km), obstimet) itrs gcrs.transform_to(ITRS(obstimet)) print(ITRS:, itrs.cartesian.xyz)对比两者如果差异在几十米到几百米量级说明手写版忽略了岁差章动极移符合预期如果差异是几千公里那多半是旋转方向或 GMST 公式错了。参数上GCRS的 obstime 要和转换时刻一致ITRS也要带 obstime否则 astropy 会用默认时间。3.4 反向转换ECEF 回 ECI 的注意点反向就是乘 R 的转置def ecef_to_eci(r_ecef, dt_utc): jd julian_date(dt_utc) theta np.radians(gmst_deg(jd)) R np.array([ [np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1] ]) return R.T np.asarray(r_ecef, dtypefloat)注意 R 是正交矩阵转置等于逆。但前提是你从 ECI 到 ECEF 用的就是同一个 R。如果中间加了岁差章动反向要把那些矩阵按逆序转置乘回去。很多人反向时直接调eci_to_ecef再取负那是错的。4. 避坑与排查ECI/ECEF 转换里最容易翻车的五件事4.1 旋转方向搞反卫星跑到地球另一边现象ECI 下卫星在 X 轴正方向转成 ECEF 后跑到 X 轴负方向附近星下点经度差了 180 度。原因R 矩阵的 sin 符号写反把顺时针当成逆时针。解决记住从 ECI 到 ECEF 是跟着地球转地球自西向东转惯性系看地固系是顺时针。用 astropy 对一遍或者拿一个已知时刻的恒星时角验证。4.2 把 UTC 当 UT1 用长弧段累积误差现象短时间转换没问题几小时后地面站指向偏了几公里。原因UTC 和 UT1 差 DUT1虽然只有不到 1 秒但地球自转 1 秒对应地表约 465 米几小时累积加上其他误差就明显了。解决精密场景下载 IERS 的 DUT1 文件用astropy.time.Time的ut1尺度低精度场景在文档里写明近似。4.3 TLE 传播结果直接当 J2000 用现象用 sgp4 算出的 TEME 坐标直接套 J2000 到 ECEF 的旋转矩阵位置偏了几十公里。原因TEME 和 J2000 差着岁差章动虽然量级不大但对高精度指向不可忽略。解决要么用 skyfield 的TEME到ITRS完整链要么自己补岁差章动矩阵。别把 TEME 当 J2000。4.4 时间标签时区丢失现象本地时间当 UTC 代入结果整体偏了 8 小时星下点经度差 120 度。原因datetime不带时区时julian_date会按字面值算实际是本地时间。解决所有时间统一用timezone.utc或者用astropy.time.Time显式指定 scale。代码里加断言检查dt.tzinfo is not None。4.5 单位混用km 和 m 串了现象位置数值对了但量级差 1000 倍天线指向完全飞掉。原因ECI 输入用 kmECEF 输出用 m中间没换算。解决在函数签名和注释里写死单位输入输出都标 km需要 m 时在调用层乘 1000。用 astropy 的u.km和u.m做量纲检查能提前发现。5. 进阶技巧用视速度验证转换正确性5.1 位置对了不代表速度对很多人只验证位置忽略速度。ECI 到 ECEF 的速度转换不是简单乘 R因为 ECEF 是旋转系还要减去地球自转带来的牵连速度v_ECEF R · v_ECI - ω × r_ECEF其中 ω 是地球自转角速度矢量在 ECEF 下是[0, 0, 7.292115e-5]rad/s。如果你只转位置不转速度或者速度直接乘 R算出来的多普勒频移和地面站跟踪角速度都会错。5.2 用数值微分做自检一个不依赖公式的自检方法取两个相邻时刻 t 和 tdt分别算 ECEF 位置做差除以 dt得到数值速度和公式算的 v_ECEF 对比。如果差在合理范围说明位置和速度转换一致。def check_velocity(r_eci, v_eci, dt_utc, dt1.0): from datetime import timedelta r1 eci_to_ecef(r_eci, dt_utc) r2 eci_to_ecef(np.array(r_eci) np.array(v_eci)*dt, dt_utc timedelta(secondsdt)) v_num (r2 - r1) / dt # 公式法 theta np.radians(gmst_deg(julian_date(dt_utc))) R np.array([[np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1]]) omega np.array([0, 0, 7.292115e-5]) v_ecef R np.array(v_eci) - np.cross(omega, r1) return v_num, v_ecef参数上dt 取 1 秒足够太小会被浮点误差淹没太大数值微分误差上升。对比两者如果差在 1e-3 km/s 量级说明实现一致。5.3 我自己的习惯我每次写完转换代码第一件事是拿一个已知答案的算例跑通比如 ISS 在某个时刻的 TLE 传播结果用 skyfield 算 ITRS 位置再和自己代码对比。第二件事是检查速度因为速度错往往比位置错更隐蔽。第三件事是在代码里留一个assert检查转换前后矢量模长是否一致——纯旋转不改变模长如果模长变了说明矩阵构造有问题。这个习惯帮我省过很多次后悔药。希望帮到你。本文还有配套的精品资源点击获取 ECI与ECEF坐标系转换：从卫星轨道到地面站指向的工程实践
简介这份资源聚焦地心惯性坐标系ECI与地心固定坐标系ECEF之间的转换面向从事卫星轨道计算、定位导航及航天器姿态分析的工程人员与相关专业学生。内容围绕地球自转对坐标的影响展开涉及儒略日期到格林尼治平恒星时的换算、恒星时与世界协调时的处理以及结合地球自转角速度求解ECEF坐标的完整思路并提及WGS84椭球模型在精度修正中的作用。压缩包共3个文件均为.m脚本体积约3KB分别承担时间参数换算与坐标转换等核心功能结构精简、便于直接嵌入MATLAB工程调用。目前已有925人学习下载适合需要快速实现ECI到ECEF转换、理解时间系统与坐标框架关系的读者参考也可作为轨道仿真与导航算法开发中的基础工具脚本使用。1. ECI 与 ECEF 坐标系转换从卫星轨道到地面站指向的落地链路做卫星轨道计算、地面站天线指向或者星上载荷几何校正的人迟早会撞上 ECI 和 ECEF 这两个坐标系。标题里那串ECITOECEF_ECI_ECI坐标系转换ECEF坐标系_eci坐标看着像随手敲的但它指向的问题非常具体给定一个在 ECI 惯性系下描述的位置或速度矢量怎么把它换算到随地球一起转的 ECEF 地固系里反过来又怎么回去。这件事在纸面上就是一个旋转矩阵但真正落地时会牵扯到时间系统、岁差章动、极移、地球自转角以及你用的到底是哪种 ECIJ2000 还是 TEME。新手容易以为套一个公式就完事熟手知道误差往往出在时间标签和参考系定义上。这篇笔记按“概念先立住、再动手复现、最后看坑”的顺序把这条转换链路拆开讲清楚适合做轨道仿真、遥感几何处理、地面站跟踪的工程师照着落地。2. ECI 与 ECEF 到底差在哪旋转矩阵背后的物理量2.1 两个坐标系的定义与选型理由ECIEarth-Centered Inertial是以地球质心为原点、坐标轴在惯性空间中基本不转的坐标系。常见的有 J2000 平赤道平春分点系和 TEMETrue Equator Mean Equinox系。J2000 的 X 轴指向 J2000.0 时刻的平春分点Z 轴指向该时刻的平北极适合做长期轨道积分和星历存储。TEME 则常用于 TLE 两行根数的传播结果它的春分点是“真”的和 J2000 之间差着岁差和章动。ECEFEarth-Centered Earth-Fixed以地球质心为原点X 轴指向本初子午线与赤道交点Z 轴指向协议地球极CTP随地球自转一起转。它适合描述地面站位置、地表目标点、卫星星下点轨迹。选型上如果你的输入是星历文件里的惯性系位置输出要给天线指向或地图叠加那 ECI 到 ECEF 就是必经之路。反过来如果你拿到的是 GPS 广播星历算出的 ECEF 位置想和惯性系下的摄动力模型对齐就要做 ECEF 到 ECI。两者不是简单转置关系因为中间夹着时间相关的旋转。2.2 转换的核心三个旋转矩阵的乘积从 ECI 到 ECEF常见做法是r_ECEF W(t) · R(t) · P(t) · N(t) · r_ECI其中 N 是章动矩阵P 是岁差矩阵R 是地球自转角矩阵W 是极移矩阵。如果用的是 J2000 到 ECEF这四项都要考虑如果精度要求不高或者时间跨度很短可以只保留 R把岁差章动极移忽略。很多工程代码里直接用一个 GMST 角构造 R 矩阵误差在角秒级对地面站指向够用对高精度定轨就不够。R 矩阵的形式是绕 Z 轴旋转地球自转角 θR [[cosθ, sinθ, 0], [-sinθ, cosθ, 0], [0, 0, 1]]注意方向从 ECI 到 ECEF 是顺时针转 θ所以矩阵里 sin 的符号和从 ECEF 到 ECI 相反。这个符号搞反是新手最常见的翻车点后面避坑章会细说。2.3 时间系统GMST 怎么算、用哪个时间θ 通常取格林尼治平恒星时 GMST。计算 GMST 需要 UT1 时间而日常拿到的是 UTC。UT1 和 UTC 差一个 DUT1量级在 0.9 秒以内。对大多数地面站指向直接用 UTC 代入 GMST 公式误差对应地表约 400 米天线波束宽的话可以接受对精密定轨必须用 UT1。一个常用的 GMST 近似公式IAU 1982GMST 280.46061837 360.98564736629 * (JD_UT1 - 2451545.0) (度)JD_UT1 是 UT1 的儒略日。算完取模 360。这个公式在 2000 年前后几十年内精度够用再往前或往后要考虑岁差模型更新。提示如果你只有 UTC又不想引入 DUT1 文件可以在代码里把 UTC 当 UT1 用但要在文档里写明这个近似别让下游以为你做了精密处理。3. 用 Python 跑通 ECI 到 ECEF 的最小实现3.1 环境准备与依赖选择我一般用 Python 做原型依赖三个库numpy做矩阵运算astropy做时间转换和坐标系定义sgp4或skyfield做 TLE 传播。如果只是验证旋转逻辑numpy加手写 GMST 就够。astropy的好处是它内置了 J2000、ITRS即 ECEF之间的完整转换链可以拿来当参考真值对比自己手写的结果。安装pip install numpy astropy如果你要处理 TLE再加pip install sgp43.2 手写旋转矩阵的完整代码下面这段代码实现从 J2000 ECI 到 ECEF 的简化转换只考虑地球自转忽略岁差章动极移。适合做教学验证和低精度场景。import numpy as np from datetime import datetime, timezone def julian_date(dt_utc): 把 UTC datetime 转成儒略日这里近似当 UT1 用 # 简化算法精度到秒级足够 y dt_utc.year m dt_utc.month d dt_utc.day frac (dt_utc.hour dt_utc.minute/60 dt_utc.second/3600) / 24.0 if m 2: y - 1 m 12 A y // 100 B 2 - A A // 4 jd int(365.25*(y4716)) int(30.6001*(m1)) d B - 1524.5 frac return jd def gmst_deg(jd_ut1): IAU 1982 GMST 公式返回度 T (jd_ut1 - 2451545.0) / 36525.0 gmst 280.46061837 360.98564736629 * (jd_ut1 - 2451545.0) return gmst % 360.0 def eci_to_ecef(r_eci, dt_utc): 输入 ECI 位置km输出 ECEF 位置km jd julian_date(dt_utc) theta np.radians(gmst_deg(jd)) # 绕 Z 轴顺时针旋转 theta R np.array([ [np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1] ]) return R np.asarray(r_eci, dtypefloat) # 测试假设某时刻 ECI 下卫星在 X 轴 7000 km dt datetime(2024, 1, 1, 0, 0, 0, tzinfotimezone.utc) r_eci [7000.0, 0.0, 0.0] r_ecef eci_to_ecef(r_eci, dt) print(ECEF:, r_ecef)逻辑说明julian_date把 UTC 转儒略日这里没有做 UT1 修正。gmst_deg用 IAU 1982 公式算格林尼治平恒星时。eci_to_ecef构造绕 Z 轴的旋转矩阵注意矩阵里 sin 的符号——从 ECI 到 ECEF 是顺时针所以第一行第二列是正 sin第二行第一列是负 sin。参数上r_eci单位是 km输出也是 kmdt_utc必须带时区否则 datetime 默认 naive算出来的儒略日会偏。3.3 用 astropy 做交叉验证手写容易错我习惯用astropy对一遍。下面代码把同一个 ECI 矢量转到 ITRS即 ECEF和上面结果对比。from astropy.time import Time from astropy.coordinates import CartesianRepresentation, GCRS, ITRS import astropy.units as u t Time(2024-01-01T00:00:00, scaleutc) # GCRS 接近 J2000 ECI gcrs GCRS(CartesianRepresentation(7000*u.km, 0*u.km, 0*u.km), obstimet) itrs gcrs.transform_to(ITRS(obstimet)) print(ITRS:, itrs.cartesian.xyz)对比两者如果差异在几十米到几百米量级说明手写版忽略了岁差章动极移符合预期如果差异是几千公里那多半是旋转方向或 GMST 公式错了。参数上GCRS的 obstime 要和转换时刻一致ITRS也要带 obstime否则 astropy 会用默认时间。3.4 反向转换ECEF 回 ECI 的注意点反向就是乘 R 的转置def ecef_to_eci(r_ecef, dt_utc): jd julian_date(dt_utc) theta np.radians(gmst_deg(jd)) R np.array([ [np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1] ]) return R.T np.asarray(r_ecef, dtypefloat)注意 R 是正交矩阵转置等于逆。但前提是你从 ECI 到 ECEF 用的就是同一个 R。如果中间加了岁差章动反向要把那些矩阵按逆序转置乘回去。很多人反向时直接调eci_to_ecef再取负那是错的。4. 避坑与排查ECI/ECEF 转换里最容易翻车的五件事4.1 旋转方向搞反卫星跑到地球另一边现象ECI 下卫星在 X 轴正方向转成 ECEF 后跑到 X 轴负方向附近星下点经度差了 180 度。原因R 矩阵的 sin 符号写反把顺时针当成逆时针。解决记住从 ECI 到 ECEF 是跟着地球转地球自西向东转惯性系看地固系是顺时针。用 astropy 对一遍或者拿一个已知时刻的恒星时角验证。4.2 把 UTC 当 UT1 用长弧段累积误差现象短时间转换没问题几小时后地面站指向偏了几公里。原因UTC 和 UT1 差 DUT1虽然只有不到 1 秒但地球自转 1 秒对应地表约 465 米几小时累积加上其他误差就明显了。解决精密场景下载 IERS 的 DUT1 文件用astropy.time.Time的ut1尺度低精度场景在文档里写明近似。4.3 TLE 传播结果直接当 J2000 用现象用 sgp4 算出的 TEME 坐标直接套 J2000 到 ECEF 的旋转矩阵位置偏了几十公里。原因TEME 和 J2000 差着岁差章动虽然量级不大但对高精度指向不可忽略。解决要么用 skyfield 的TEME到ITRS完整链要么自己补岁差章动矩阵。别把 TEME 当 J2000。4.4 时间标签时区丢失现象本地时间当 UTC 代入结果整体偏了 8 小时星下点经度差 120 度。原因datetime不带时区时julian_date会按字面值算实际是本地时间。解决所有时间统一用timezone.utc或者用astropy.time.Time显式指定 scale。代码里加断言检查dt.tzinfo is not None。4.5 单位混用km 和 m 串了现象位置数值对了但量级差 1000 倍天线指向完全飞掉。原因ECI 输入用 kmECEF 输出用 m中间没换算。解决在函数签名和注释里写死单位输入输出都标 km需要 m 时在调用层乘 1000。用 astropy 的u.km和u.m做量纲检查能提前发现。5. 进阶技巧用视速度验证转换正确性5.1 位置对了不代表速度对很多人只验证位置忽略速度。ECI 到 ECEF 的速度转换不是简单乘 R因为 ECEF 是旋转系还要减去地球自转带来的牵连速度v_ECEF R · v_ECI - ω × r_ECEF其中 ω 是地球自转角速度矢量在 ECEF 下是[0, 0, 7.292115e-5]rad/s。如果你只转位置不转速度或者速度直接乘 R算出来的多普勒频移和地面站跟踪角速度都会错。5.2 用数值微分做自检一个不依赖公式的自检方法取两个相邻时刻 t 和 tdt分别算 ECEF 位置做差除以 dt得到数值速度和公式算的 v_ECEF 对比。如果差在合理范围说明位置和速度转换一致。def check_velocity(r_eci, v_eci, dt_utc, dt1.0): from datetime import timedelta r1 eci_to_ecef(r_eci, dt_utc) r2 eci_to_ecef(np.array(r_eci) np.array(v_eci)*dt, dt_utc timedelta(secondsdt)) v_num (r2 - r1) / dt # 公式法 theta np.radians(gmst_deg(julian_date(dt_utc))) R np.array([[np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1]]) omega np.array([0, 0, 7.292115e-5]) v_ecef R np.array(v_eci) - np.cross(omega, r1) return v_num, v_ecef参数上dt 取 1 秒足够太小会被浮点误差淹没太大数值微分误差上升。对比两者如果差在 1e-3 km/s 量级说明实现一致。5.3 我自己的习惯我每次写完转换代码第一件事是拿一个已知答案的算例跑通比如 ISS 在某个时刻的 TLE 传播结果用 skyfield 算 ITRS 位置再和自己代码对比。第二件事是检查速度因为速度错往往比位置错更隐蔽。第三件事是在代码里留一个assert检查转换前后矢量模长是否一致——纯旋转不改变模长如果模长变了说明矩阵构造有问题。这个习惯帮我省过很多次后悔药。希望帮到你。本文还有配套的精品资源点击获取