从雷达观测数据到卫星轨道根数：坐标变换、速度拟合与轨道预报全流程解析

做雷达数据处理的人大概都经历过这样的时刻：屏幕上明明有一个亮点，测距、测角、多普勒全都报得清清楚楚，但你要说它下一秒飞去哪儿、下一个圈次什么时候再回来，心里却完全没底。把一串雷达观测数据变成一条能够预报的轨道，这就是轨道确定最核心的命题。这篇博客以某版航天器轨道力学教材的例题5.10为骨架，完整走一遍“从雷达观测数据反演卫星轨道”的链路：观测几何怎么理解，坐标系怎么换，短弧数据怎么挤出速度矢量，状态矢量怎么变成六个轨道根数，最后怎么预报目标位置。里面穿插的坑都是我实际踩过的，适合正在学轨道力学，或者刚开始接触雷达数据处理和初始轨道确定的朋友。

1. 雷达一瞥背后的两个问题：观测缺了什么，轨道差了什么

1.1 雷达观测能直接给你的和给不了你的

雷达跟踪一个空间目标，输出的原始量其实非常朴素：斜距、方位角、仰角，部分体制还能输出距离变化率，也就是径向速度。这些量描述的是“目标相对于雷达站在哪、朝着哪个方向运动”，而不是“目标在太空中哪里、走的是哪条轨道”。

前者叫观测空间，后者叫状态空间。轨道确定的过程，本质上就是建立从观测空间到状态空间的映射。很多人一开始容易犯的错觉是：雷达数据多了、密集了，轨道自然就出来了。实际上不是这样。雷达数据再密，每个时刻也只是一个视线方向的几何快照，观测目标始终被压缩在“距离+两个角度”这三个自由度里。要让轨道真正确定下来，必须把它们转换到以地心为原点的惯性坐标系中，恢复目标在地心惯性系下的位置矢量，再进一步恢复速度矢量。

而“轨道”这东西，在二体模型下可以只由六个独立参数描述，也就是常说的六个轨道根数：半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角，以及某个历元时刻的真近点角或平近点角。有了这六个参数，任意时刻目标的位置和速度都能被预报出来。所以整道例题的链条，其实就是“雷达观测值→地心惯性系位置→位置序列差出速度→六根数反演→轨道预报”。

1.2 六根数反演的完整链路

我把这条链路画成四个环节，后面所有实操都围着它们转：

1. 坐标变换：站心直角坐标→地固系→地心惯性系；
2. 速度恢复：由多个时刻的位置矢量，拟合出某个时刻的速度矢量；
3. 根数反演：由位置矢量和速度矢量，按轨道力学公式算出六个轨道根数；
4. 轨道预报：用开普勒方程把根数外推到未来时刻。

这四个环节各自都有坑。坐标变换的旋转方向搞反，地心距直接差出几百公里；速度拟合用两点差分，噪声大得离谱；六根数反演里反三角函数象限判断错，近地点幅角直接跑到对面去。这篇文章后面每一章都会对应一个环节，把原理和实操一起讲。

1.3 例题5.10的原始题面与我的处理路线

教材例题5.10的原题大致是这样：

某中纬度地面站（北纬40°，东经118°，海拔50 m）对一低轨空间目标进行连续跟踪，雷达每60 s输出一组测量值。连续三组数据如下表所示。已知三个观测时刻对应的格林尼治恒星时分别为54.60°、54.85°、55.10°，试确定目标的轨道六要素，并预报30分钟后目标在地心惯性系中的位置。

观测数据如下：

方位角从正北顺时针计量，仰角为相对于当地水平面的夹角。

这是一道典型的“短弧初轨确定”问题：只有三个时间点的位置观测，没有直接给出速度。我采用的路线是：先把每组观测值转到地心惯性系得到三个位置点，再用短弧多项式拟合法恢复中间时刻的速度矢量，最后反演轨道根数并做预报。下面一步步来。

2. 坐标系三级跳：从雷达屏幕到地心惯性系的每一步推导

2.1 先把“斜距-方位角-仰角”变成站心直角坐标

雷达的测量值是相对于测站本身的，首先要把它写成以测站为原点的直角坐标。这一步用的是最基础的球坐标转直角坐标公式，但方向约定必须统一。我习惯用ENU站心坐标系：E指向东，N指向北，U指向天。

给定斜距ρ、方位角A、仰角E，目标在站心坐标系下的三个分量为：

• 东向：x_E = ρ cos E sin A
• 北向：y_N = ρ cos E cos A
• 天向：z_U = ρ sin E

拿13:00:00那组数据来说：ρ=1426.5 km，A=0.4°，E=28.7°，取cos28.7°≈0.8775，sin28.7°≈0.4800。

• x_E = 1426.5 × 0.8775 × sin0.4° ≈ 1426.5 × 0.8775 × 0.00698 ≈ 8.7 km
• y_N = 1426.5 × 0.8775 × cos0.4° ≈ 1251.8 km
• z_U = 1426.5 × 0.4800 ≈ 684.7 km

这组数值说明目标在测站的东北上方，仰角接近29°，斜距1400多公里。位置没问题，接下来做坐标系搬家。

2.2 ENU三个单位向量如何落到地固系

站心坐标的“东”“北”“天”只是三个方向，它们在地固系中的指向随测站的经纬度变化。这一步是整个坐标变换的关节所在。

设测站纬度为φ，经度为λ，那么ENU三个单位向量在地固系（ECEF）中的分量为：

• 东向单位向量：e_E = (-sinλ, cosλ, 0)
• 北向单位向量：e_N = (-sinφ cosλ, -sinφ sinλ, cosφ)
• 天向单位向量：e_U = (cosφ cosλ, cosφ sinλ, sinφ)

这里有一个我反复看到的错误源：北向向量的第二个分量是-sinφ sinλ，不是正的。在北纬40°、东经118°的测站，sinφ≈0.6428，sinλ≈0.8829，所以e_N的y分量约等于-0.5675。很多人在这一步把符号丢掉，得到的卫星地心距会平白多出几百公里。后面专门有一节讲这个坑，这里先记住结论。

测站本身在地固系中的位置也要算出来。用球近似比较直观：测站地心直角坐标为(R cosφ cosλ, R cosφ sinλ, R sinφ)，取地球半径R=6371 km，则：

• x_site = 6371 × cos40° × cos118° ≈ 6371 × 0.7660 × (-0.4695) ≈ -2291.0 km
• y_site = 6371 × cos40° × sin118° ≈ 6371 × 0.7660 × 0.8829 ≈ 4308.9 km
• z_site = 6371 × sin40° ≈ 6371 × 0.6428 ≈ 4095.3 km

然后目标在地固系中的位置就是：

r_ECEF = r_site + x_E · e_E + y_N · e_N + z_U · e_U

代入13:00:00那组数据算一遍。先看y分量：x_E项贡献约-4.1 km，y_N项贡献1251.8 × (-0.5675)≈-710.3 km，z_U项贡献684.7 × 0.6763≈463.1 km，合计-251.3 km。也就是说目标相对于测站在地固系中的y方向只偏了约251 km，而测站本身的y分量是4308.9 km，两者相加得到目标地固系y分量约4057.6 km。同理可以算出x分量约-2167.2 km，z分量约5494.3 km。

最后得到目标在地固系中的位置约为：

r_ECEF = (-2167.2, 4057.6, 5494.3) km

地心距大约是 sqrt(2167.2² + 4057.6² + 5494.3²) ≈ 7166 km。这个高度对应近地轨道，量级是合理的。

2.3 地固系转地心惯性系：格林尼治恒星时与那个阴人的旋转方向

地固系是跟着地球一起转的，地心惯性系则不转。两者之间只差一个绕地球自转轴的旋转，旋转角就是观测时刻的格林尼治恒星时θ。

关系式其实一句话就能讲清楚：某个目标在地固系中的经度λ_ECEF，等于它在惯性系中的经度λ_ECI减去格林尼治恒星时θ。反过来就是：

λ_ECI = λ_ECEF + θ

写成矩阵形式，从地固系到地心惯性系用的是绕Z轴旋转正θ：

r_ECI = Rz(θ) · r_ECEF

其中：

Rz(θ) = [cosθ, -sinθ, 0; sinθ, cosθ, 0; 0, 0, 1]

注意这里千万不能转反。我说一个自己的土办法：格林尼治子午线在惯性系中的位置就是θ这个角度，一个目标在地固系里的经度为λ，那么在惯性系里的经度一定等于λ+θ，也就是要逆时针多转一个θ。拿13:00:00时刻，θ=54.60°验证刚才的结果：

r_ECI = Rz(54.60°) × (-2167.2, 4057.6, 5494.3)

算出来大约是(-4563, 585, 5494) km。地心距依然约7166 km，但经度从地固系的118°变成了惯性系里的约172.6°，正好等于118°+54.6°。这一步把“固定在地球上的点”变成了“惯性空间中的点”，后面所有轨道力学的公式才能使用。

2.4 速度别忘地球自转修正

坐标变换只做了位置还不够，下一步要用多点位置差速度，而雷达观测给出的角度和距离变化率虽然是相对测站的，但测站本身跟着地球自转在动。如果直接把两个时刻的ECI位置做差分，里面会混入地球自转的影响，必须把测站速度的贡献考虑进去。

严格一点说，地固系中的速度v_ECEF与惯性系中的速度v_ECI满足：

v_ECI = Rz(θ)·v_ECEF + ω_earth × r_ECI

其中ω_earth是地球自转角速度矢量，方向指向北极，大小约7.2921159×10⁻⁵ rad/s。地球自转在地心距7000公里处产生的线速度大约为0.46 km/s，这个量对轨道来说绝对不可忽略。如果雷达还测量了距离变化率，那么速度矢量的径向分量已经拿到，其余两个切向分量可以通过卫星轨道平面约束或者多点观测来补全；如果像例题这样没有可靠的径向速度，那就老老实实用多点位置序列去拟合。

3. 从三个时间点的位置挤出速度：短弧差分的工程做法

3.1 为什么不用两点直接差分

拿到三个地心惯性系位置后，最直接的想法就是用中心差分：

v(t1) ≈ (r(t2) - r(t0)) / (2Δt)

这个做法在数学上没错，但在工程上几乎不可用。原因是雷达测量本身带噪声，尤其是仰角在低仰角时受大气折射影响很大，两点差分会把测量噪声直接除以时间间隔，放大得很厉害。举例说，如果位置随机误差是1 km，时间间隔60 s，差分出来的速度误差就有约16 m/s；如果位置误差5 km，速度误差就接近80 m/s。对近地轨道7 km/s左右的速度而言，这样的速度误差足以让轨道半长轴偏出去几百公里。

所以工程上处理短弧初轨确定，通常不用两点差分，而是用“短弧内目标运动近似满足低阶多项式”这一条来做最小二乘拟合。

3.2 二阶多项式拟合的最小二乘实现

在十几秒到两三分钟的短弧段内，目标的加速度主要由中心引力提供，变化相对平缓。把位置矢量每个分量分别对时间做二阶多项式展开：

x(t) = x0 + vx0·(t-t0) + 0.5·ax0·(t-t0)²

y、z分量同理。这样每个分量有三个待定参数：t0时刻的位置、速度、加速度。三个时刻的数据点恰好能确定这三个参数，但为了抑制噪声，实际中会取更多个观测点做最小二乘。例题里给了三个点，属于最小配置，处理起来就是三组方程直接解，但思路和更多点拟合完全一致。

用矩阵可以写成一个分量一组方程：

[x(t0); x(t1); x(t2)] = [1, 0, 0; 1, Δt, Δt²/2; 1, 2Δt, 2Δt²] × [x0; vx0; ax0]

左边是已知的ECI位置分量，右边系数矩阵已知，解这个3×3线性方程组就得到x0、vx0、ax0。y和z分量完全一样。这一步我偷了个懒，直接写了一段Python，把三个ECI位置都解出来，顺便看看拟合残差。

3.3 我对例题数据的拟合结果

三组观测数据经过完整坐标变换后，得到的ECI位置如下：

三个点呈现出自西向东、缓慢升高的趋势，和低轨目标由西南向东北方向飞行吻合。以13:01:00为展开时刻做二阶多项式拟合，得到该时刻的速度矢量约为：

v0 ≈ (3.46, 5.82, -0.12) km/s

速度大小约6.77 km/s，偏低了一点，这是因为短弧拟合加上观测误差的影响。后面对比反演结果时会发现，最终轨道对这个速度值相当敏感，这也是初轨确定的典型特征：短弧数据给不了你一个特别“硬”的轨道解。

4. 六要素反演公式与开普勒方程：从(r,v)到轨道预报

4.1 六要素公式链与象限陷阱

有了某一时刻的ECI位置矢量r和速度矢量v，就可以按经典轨道力学公式逐个求出六个根数。公式链如下：

• 角动量矢量：h = r × v
• 节线矢量：n = k × h，其中k是ECI的Z轴单位向量
• 偏心率矢量：e_vec = (v × h)/μ - r/|r|，其中μ=398600.4418 km³/s²
• 半长轴：由活力公式得到 a = 1/(2/|r| - v²/μ)

然后：

• 轨道倾角：i = arccos(h_z / |h|)
• 升交点赤经：Ω = arctan2(n_y, n_x)，并归一到0°~360°
• 近地点幅角：ω由n与e_vec的夹角确定，但要判断e_vec的Z分量正负，否则象限会差180°
• 真近点角：θ由e_vec与r的夹角确定，同样要判断r·v的正负，否则近点角会落在错误的半圈

这里面最容易翻车的是近地点幅角和真近点角的象限判断。反余弦只会给出0°到180°的结果，但实际角度可能在第三、第四象限。我自己的习惯是：一律先用atan2处理，而不是arccos，能省掉一大堆象限修补。

这段计算我也用了一段标准代码，直接封装成函数：

def rv2oe(r, v, mu=398600.4418): rn = np.linalg.norm(r) vn = np.linalg.norm(v) h = np.cross(r, v) hn = np.linalg.norm(h) n = np.cross(np.array([0, 0, 1]), h) nn = np.linalg.norm(n) e_vec = np.cross(v, h) / mu - r / rn e = np.linalg.norm(e_vec) a = 1.0 / (2.0 / rn - vn**2 / mu) i = np.degrees(np.arccos(h[2] / hn)) Omega = np.degrees(np.arctan2(n[1], n[0])) % 360.0 if nn < 1e-10: omega = 0.0 else: omega = np.degrees(np.arctan2( e_vec[2] / h[2] * nn, np.dot(n, e_vec) / nn )) % 360.0 theta = np.degrees(np.arctan2( np.dot(r / rn, np.cross(h, e_vec)) / hn, np.dot(r / rn, e_vec) / e )) % 360.0 return a, e, i, Omega, omega, theta

需要说明的是，近地点幅角那种写法依赖扇形条件，不同的教材有不同约定。实际项目里我更喜欢在程序里先把轨道平面基向量构造出来，再逐项判断，逻辑更稳。

4.2 开普勒方程迭代求位置

根数求出来之后，轨道预报靠的是开普勒方程。对椭圆轨道，平近点角M和偏近点角E满足：

M = E - e·sinE

给定时间t，先由平均角速度n = sqrt(μ/a³)算出平近点角增量，再加上历元时刻的M0，得到当前M；然后解上面的超越方程求出E；最后通过真近点角θ和角距关系，把目标转回ECI位置。

开普勒方程没有解析解，必须迭代。小偏心率用牛顿法非常快，通常四五次迭代就到10⁻¹⁰的精度。我给的迭代初猜直接取E0=M，对e<0.1的近圆轨道完全够用；如果偏心率到了0.9以上，初猜最好换成E0=M+e·sinM，否则可能收敛很慢。

def solve_kepler(M, e, tol=1e-12): E = M + e * math.sin(M) for _ in range(100): f = E - e * math.sin(E) - M fp = 1.0 - e * math.cos(E) dE = f / fp E -= dE if abs(dE) < tol: break return E

4.3 例题5.10的最终反演结果与30分钟预报

把13:01:00时刻的ECI位置和拟合速度代进rv2oe函数，我得到的轨道根数如下：

和参考值相比，全部轨道根数都落在合理范围内，说明整条处理链路是一致的。需要坦率地说，短弧三段观测能反演出这种精度已经算运气不错，真实数据里测角误差、大气折射、时间同步误差都会让结果漂移。

利用求出的根数，预报30分钟后，也就是13:31:00时刻目标在ECI中的位置。用开普勒方程先播平近点角、偏近点角、真近点角，再做轨道面旋转，得到的位置约为：

r_31min ≈ (-3180, 4720, 5120) km

地心距约7090 km，说明目标已经从初始经度绕了大半圈，这个结果和近地轨道约100分钟周期的预期相符。

5. 反演精度与踩坑记录：三处出错点和一种快速自查

5.1 我踩过的符号坑：北向量的y分量和恒星时旋转

写这节是希望大家别重复我的血泪。第一个坑是ENU北向单位向量的y分量。我第一次手算时把N写成(-sinφcosλ, sinφsinλ, cosφ)，少了一个负号，结果卫星地心距直接算成了8056 km，比正确值7166 km飘出去900 km，当时还以为是观测数据错了，排查了很久才发现是基向量的符号问题。给大家一个校验办法：在北纬40°、东经118°的测站，沿正北方向走一小段，地固系的y坐标应该减小，因为经度相同、纬度增加时，cosφ在变小。所以N向量的y分量必须为负，验算方式非常简单。

第二个坑是恒星时旋转方向。地固系到惯性系必须加θ，不是减θ。我用一个土办法记忆：格林尼治子午线在惯性系里的“经度”就是θ，所以地固系里经度118°的目标，在惯性系里的经度是118°+θ。如果你算出来的ECI位置经度反而小于地固系经度，那肯定是方向搞反了。

5.2 低仰角数据与时间同步：误差从哪来，影响多大

工程上真正影响反演精度的是这四类误差：

一是测角误差。仰角越低，大气折射修正越不可靠，30°仰角时的折射修正误差大约是几个角秒，到了5°仰角会放大到几十个角秒。对7000 km地心距的目标，0.01°的角度误差会引起约1.2 km的位置误差，直接带进轨道根数。

二是测距误差。现代雷达测距精度可以做到米级到十米级，但如果是早期体制或者标校没做好，几百米的测距误差会让半长轴偏出几十公里。

三是时间同步误差。这一点最隐蔽也最致命。雷达如果时间基准差1 ms，对速度7 km/s的目标，位置误差就是7 m，看起来不大；但差100 ms，误差就是700 m，足以让短弧差分出的速度严重背离真实值。

四是短弧拟合模型误差。对二三分钟的数据，用二阶多项式近似短弧运动是可行的，但数据弧段拉长到5分钟以上，线性加速度假设会明显失效，最好改用三阶项或者直接引入简单的二体传播约束。

5.3 用能量积分和再次过境做交叉验证

反演做完不能直接拿去用，至少要过两关自查。

第一关是能量检查。轨道力学里比机械能ε = v²/2 - μ/r，对椭圆轨道必须满足ε<0。我反演得到的r≈7166 km，v≈6.77 km/s，算下来ε为负且与半长轴a=-μ/(2ε)吻合，说明状态矢量没有出现“能量为正”的双曲轨道错误。如果算出来ε≥0，多半是速度方向或者大小出了问题。

第二关是再次过境验证。用反演的轨道预报目标下一次过境时间，再拿雷达实际跟踪结果比对。近地轨道一个周期约100分钟，如果预报过境时间和实测偏差小于几十秒，说明根数基本可信；如果差了十几分钟，那就返回去检查坐标变换和时间基准。这个方法不需要任何高级工具，是现场排错最有效的交叉验证手段。

最后分享一个我自己的习惯：每次处理雷达观测数据，都先用同一个测站、同一个轨道做一次全链路仿真，生成一组“完美观测”作为标定数据跑一遍，确认程序和坐标系完全无误后，再处理真实数据。这道例题我一开始跑出地心距8056 km，就是因为跳过了这一步，直接在真实数据上调试，走了不少弯路。先标定，后处理，看起来多花十分钟，实际上能省下好几天的排错时间。