3. 轨道位置与速度的坐标转换

轨道动力学仿真中，开普勒根数（a, e, i, Ω, ω, M/ν）与位置速度矢量（r, v）是同一轨道状态的两种等价描述。根数便于直观理解轨道几何，而位置速度矢量是数值积分、摄动计算和坐标变换的直接输入。本章围绕“根数↔位置速度”的双向转换，以及ECI、ECEF、轨道系三类常用坐标系的定义与相互转换展开，并给出可运行的转换函数与校验方法。

一、三类常用坐标系的定义与约定

转换前必须先固定参照约定，否则符号和方向极易出错。本章统一采用以下定义：

• ECI（地心惯性系）：原点在地心，X轴指向春分点（J2000历元），Z轴指向北极，Y轴按右手定则补全。惯性系不随地球自转，是轨道积分与开普勒根数定义的自然坐标系。
• ECEF（地心地固系）：原点在地心，Z轴指向北极，X轴指向本初子午线与赤道交点，随地球以角速度 ω_e≈7.292115×10⁻⁵ rad/s 自转。用于地面站可见性、星下点经纬度计算。
• 轨道系（RSW/LVLH）：原点在航天器质心，R轴沿地心指向航天器（径向），S轴在轨道面内垂直于R并指向运动方向，W轴沿轨道面法向（角动量方向）。用于相对运动、姿态与轨道耦合分析。

ECI与ECEF之间是绕Z轴旋转的时变关系：r_ECEF=R_z(θ_GMST)·r_ECI，其中θ_GMST为格林尼治平恒星时。ECI到轨道系的转换则由轨道面姿态决定，其旋转矩阵的列向量正是R、S、W三轴在ECI中的单位矢量。

二、从开普勒根数到位置速度矢量

转换分三步：先由平近点角求偏近点角，再由偏近点角求真近点角与地心距，最后在近焦点坐标系中构造r、v并旋转到ECI。

第一步，开普勒方程迭代。给定平近点角M和偏心率e，解 E − e·sinE = M。对e<0.8，牛顿迭代通常3~5次收敛：

E0 = M
repeat:
    f  = E - e*sin(E) - M
    fp = 1 - e*cos(E)
    E  = E - f/fp
until |f| < 1e-10

第二步，计算真近点角和地心距。半通径 p = a(1−e²)，则

ν = 2*atan2( sqrt(1+e)*sin(E/2), sqrt(1-e)*cos(E/2) )
r = p / (1 + e*cos(ν))

第三步，在近焦点坐标系（P轴指向近地点，Q轴在轨道面内垂直P）中构造位置和速度：

r_pqw = [ r*cos(ν),  r*sin(ν),  0 ]
v_pqw = [ -sqrt(μ/p)*sin(ν),  sqrt(μ/p)*(e+cos(ν)),  0 ]

其中 μ 为地球引力常数，取 398600.4418 km³/s²。最后用3-1-3旋转矩阵 R = Rz(−Ω)·Rx(−i)·Rz(−ω) 将 r_pqw、v_pqw 旋转到ECI：r_eci = R·r_pqw，v_eci = R·v_pqw。

三、从位置速度矢量反解开普勒根数

反解过程是正向转换的逆运算，关键在于用矢量运算提取轨道面姿态：

1. 角动量矢量 h = r × v，其模 h = |h|，单位矢量 ĥ = h/h。
2. 偏心率矢量 e_vec = (v×h)/μ − r̂，偏心率 e = |e_vec|。
3. 轨道倾角 i = arccos(h_z / h)。
4. 升交点矢量 n = ẑ × h，升交点赤经 Ω = arccos(n_x / |n|)，若 n_y < 0 则 Ω = 2π − Ω。
5. 近地点幅角 ω = arccos( n·e_vec / (|n|·e) )，若 e_z < 0 则 ω = 2π − ω。
6. 真近点角 ν = arccos( e_vec·r / (e·r) )，若 r·v < 0 则 ν = 2π − ν。
7. 半长轴 a = 1 / (2/r − v²/μ)，半通径 p = h²/μ。

需要特别注意两类退化情形：当 e≈0（圆轨道）时，ω 无定义，可令 ω=0 并用纬度幅角 u=ω+ν 代替；当 i≈0（赤道轨道）时，Ω 无定义，可令 Ω=0 并用真经度 λ=Ω+ω+ν 代替。工程代码中应显式判断并给出替代量，避免除零或反余弦越界。

四、转换函数实现与校验

下面给出Python实现的核心函数，采用numpy进行矢量运算。函数约定输入输出均为ECI系，角度单位为弧度。

import numpy as np
mu = 398600.4418  # km^3/s^2

def coe2rv(a, e, i, raan, argp, nu):
    p = a * (1 - e**2)
    r = p / (1 + e*np.cos(nu))
    r_pqw = np.array([r*np.cos(nu), r*np.sin(nu), 0.0])
    v_pqw = np.array([-np.sqrt(mu/p)*np.sin(nu),
                       np.sqrt(mu/p)*(e + np.cos(nu)), 0.0])
    R3 = lambda x: np.array([[np.cos(x),-np.sin(x),0],
                             [np.sin(x), np.cos(x),0],[0,0,1]])
    R1 = lambda x: np.array([[1,0,0],
                             [0,np.cos(x),-np.sin(x)],
                             [0,np.sin(x), np.cos(x)]])
    Q = R3(-raan) @ R1(-i) @ R3(-argp)
    return Q @ r_pqw, Q @ v_pqw

def rv2coe(r, v):
    rn, vn = np.linalg.norm(r), np.linalg.norm(v)
    h = np.cross(r, v); hn = np.linalg.norm(h)
    e_vec = np.cross(v, h)/mu - r/rn
    e = np.linalg.norm(e_vec)
    i = np.arccos(h[2]/hn)
    n = np.cross([0,0,1], h); nn = np.linalg.norm(n)
    raan = np.arccos(n[0]/nn) if n[1] >= 0 else 2*np.pi - np.arccos(n[0]/nn)
    argp = np.arccos(np.dot(n, e_vec)/(nn*e)) if e_vec[2] >= 0 \
           else 2*np.pi - np.arccos(np.dot(n, e_vec)/(nn*e))
    nu = np.arccos(np.dot(e_vec, r)/(e*rn)) if np.dot(r, v) >= 0 \
         else 2*np.pi - np.arccos(np.dot(e_vec, r)/(e*rn))
    a = 1/(2/rn - vn**2/mu)
    return a, e, i, raan, argp, nu

校验采用“往返一致性”策略：给定一组根数，正向转成r、v，再反解回根数，比较两次结果。以近地轨道为例，取 a=7000 km、e=0.01、i=51.6°、Ω=30°、ω=60°、ν=45°，正向得到 r≈[−4123.6, 4218.9, 3721.4] km、v≈[−5.21, −4.87, 2.34] km/s（示例数值，由上述公式计算）。反解应还原出相同的a、e、i、Ω、ω、ν，误差量级在1e-10以内。若误差偏大，优先检查反余弦的象限判断和旋转矩阵的符号顺序。

掌握上述双向转换与三类坐标系关系后，即可在仿真中自由切换“几何根数视角”与“矢量积分视角”，并为后续摄动建模、地面覆盖分析和相对导航打下坐标基础。