2. 坐标系与运动学：常用坐标系、坐标变换、火箭运动方程

搞火箭制导，第一件事不是写代码，而是把坐标系搞明白。

我见过太多新手，算法推得头头是道，一跑仿真就炸了。为什么？坐标系用错了。说白了，火箭在天上飞，你在地面看，你看到的「上下左右」和火箭感受到的「前后左右」完全是两码事。今天我们就来把这些坐标系掰扯清楚。

2.1 常用坐标系

火箭制导里常用的坐标系，其实就三个。我习惯把它们分成「地面看的」、「火箭想的」和「中间过渡的」。

2.1.1 发射坐标系

这个坐标系是固定在地面上的。原点通常选在发射点。X轴指向发射方位（正东或正北，看具体任务），Y轴垂直向上，Z轴按右手定则确定。

嗯，这里要注意：发射坐标系不是惯性系。它跟着地球自转在动。你想想看，火箭刚起飞时，它相对于发射坐标系的速度是0，但相对于地心，它已经有了一部分地球自转的速度。

我在项目中遇到过一个坑：某次仿真，我们把发射坐标系当成了惯性系来用，结果弹道偏差越来越大。后来排查才发现，地球自转带来的科里奥利力被忽略了。所以，发射坐标系只适合描述火箭相对于地面的位置和姿态，做动力学分析时一定要小心。

2.1.2 惯性坐标系

这才是真正的「不动」坐标系。通常我们用地心惯性系（ECI），原点在地心，X轴指向春分点，Z轴指向北极。

为什么需要它？因为牛顿定律只在惯性系里成立。你写运动方程，必须用惯性系。否则，那些虚拟力（离心力、科里奥利力）会让你头大。

我个人习惯：所有动力学计算都在惯性系里做，只在最后输出结果时，才转换到发射坐标系或速度坐标系。这样逻辑清晰，不容易出错。

2.1.3 速度坐标系

这个坐标系跟着火箭走。原点在火箭质心，X轴沿着速度方向，Y轴在纵向对称面内垂直于X轴，Z轴按右手定则。

为什么要搞这么个坐标系？因为气动力（升力、阻力）的定义和速度方向直接相关。阻力永远和速度方向相反，升力垂直于速度方向。你在速度坐标系里写气动力方程，简直不要太舒服。

2.2 坐标变换

坐标系之间怎么转？说白了就是旋转矩阵。我给大家总结一下最常用的几个变换。

2.2.1 发射坐标系 → 惯性坐标系

这个变换要考虑地球自转。假设发射时刻为t₀，地球自转角速度为ωₑ，那么从发射坐标系到惯性坐标系的变换矩阵为：

# Python示例：发射坐标系到惯性坐标系的变换
import numpy as np

def launch_to_inertial(lat, lon, t):
    """
    lat: 发射点纬度 (rad)
    lon: 发射点经度 (rad)
    t:   从发射开始的时间 (s)
    """
    omega_e = 7.2921150e-5  # 地球自转角速度 (rad/s)
    
    # 先绕Z轴转经度
    Rz_lon = np.array([
        [np.cos(lon), -np.sin(lon), 0],
        [np.sin(lon),  np.cos(lon), 0],
        [0,            0,           1]
    ])
    
    # 再绕Y轴转纬度
    Ry_lat = np.array([
        [np.cos(np.pi/2 - lat), 0, np.sin(np.pi/2 - lat)],
        [0,                     1, 0],
        [-np.sin(np.pi/2 - lat), 0, np.cos(np.pi/2 - lat)]
    ])
    
    # 最后考虑地球自转
    Rz_rot = np.array([
        [np.cos(omega_e * t), -np.sin(omega_e * t), 0],
        [np.sin(omega_e * t),  np.cos(omega_e * t), 0],
        [0,                   0,                    1]
    ])
    
    return Rz_rot @ Ry_lat @ Rz_lon
小技巧：我建议把坐标变换写成函数，每次调用时传入时间参数。这样仿真时，每个时间步都能自动更新变换矩阵，不会出错。

2.2.2 惯性坐标系 → 速度坐标系

这个变换需要知道火箭当前的速度方向。假设速度矢量在惯性系中的分量为(vx, vy, vz)，那么：

def inertial_to_velocity(v):
    """
    v: 速度矢量 (3x1)
    返回：从惯性系到速度系的变换矩阵
    """
    v_norm = np.linalg.norm(v)
    if v_norm < 1e-10:
        return np.eye(3)
    
    # 速度方向单位矢量
    v_hat = v / v_norm
    
    # 构造速度坐标系基矢量
    # X轴：速度方向
    x_v = v_hat
    
    # Y轴：在纵向对称面内，垂直于X轴
    # 这里假设纵向对称面由速度矢量和当地垂线决定
    # 实际工程中需要更精确的定义
    z_up = np.array([0, 0, 1])  # 简化处理
    y_v = np.cross(z_up, x_v)
    y_v = y_v / np.linalg.norm(y_v)
    
    # Z轴：右手定则
    z_v = np.cross(x_v, y_v)
    
    # 变换矩阵（从惯性系到速度系）
    C = np.vstack([x_v, y_v, z_v])
    return C
曾经踩过的坑：速度坐标系的定义在不同文献里可能不一样。有的把Y轴指向天，有的指向侧向。做项目前，一定要和团队统一约定。我曾经因为坐标系定义不一致，和同事对了两天的数据才找到问题。

2.3 火箭运动方程

好了，坐标系和变换都准备好了，现在可以写运动方程了。火箭的运动方程，说白了就是牛顿第二定律加上转动方程。

2.3.1 质心运动方程

在惯性坐标系下，火箭质心的运动方程为：

m * dv/dt = F_gravity + F_thrust + F_aero + F_control

其中：

  • F_gravity：地球引力。注意不是简单的mg，要考虑地球形状和高度变化。我一般用J2模型。
  • F_thrust：发动机推力。方向沿火箭纵轴，大小由发动机特性决定。
  • F_aero：气动力。包括阻力和升力，在速度坐标系里计算最方便。
  • F_control：控制力。来自舵面偏转或推力矢量。

2.3.2 绕质心转动方程

火箭的转动用欧拉方程描述：

I * dω/dt + ω × (I * ω) = M_aero + M_control

这里I是惯量张量，ω是角速度。M_aero是气动力矩，M_control是控制力矩。

核心要点：转动方程必须在箭体坐标系下写。因为惯量张量在箭体系里是常数（假设燃料消耗引起的质量变化已考虑），在其他坐标系里会随时间变化，那计算就复杂了。

2.3.3 完整的仿真流程

我给大家画个流程图，看看整个仿真是怎么串起来的：

火箭六自由度仿真流程 初始化状态 计算引力、推力、气动力/力矩 坐标变换（统一到惯性系） 数值积分（RK4） 输出状态（位置、速度、姿态）

2.4 代码实现示例

最后，给大家一个完整的运动方程积分示例。这个代码我用了很多年，结构清晰，适合做基础框架：

import numpy as np
from scipy.integrate import solve_ivp

def rocket_dynamics(t, state, params):
    """
    火箭六自由度动力学方程
    state: [x, y, z, vx, vy, vz, q0, q1, q2, q3, wx, wy, wz]
    params: 火箭参数（质量、惯量、推力等）
    """
    # 解包状态
    pos = state[0:3]
    vel = state[3:6]
    quat = state[6:10]  # 四元数表示姿态
    omega = state[10:13]  # 角速度
    
    # 计算引力（简化模型）
    r = np.linalg.norm(pos)
    g = -params['mu'] / r**3 * pos
    
    # 计算推力（假设沿箭体纵轴）
    thrust_body = np.array([params['thrust'], 0, 0])
    # 从箭体系到惯性系
    C_body_to_inertial = quat_to_dcm(quat)
    thrust_inertial = C_body_to_inertial @ thrust_body
    
    # 计算气动力（在速度系中）
    v_norm = np.linalg.norm(vel)
    if v_norm > 1e-6:
        # 动压
        q = 0.5 * params['rho'] * v_norm**2
        # 阻力（沿速度反方向）
        drag = -q * params['Cd'] * params['S'] * vel / v_norm
    else:
        drag = np.zeros(3)
    
    # 总力
    F_total = g * params['mass'] + thrust_inertial + drag
    
    # 加速度
    acc = F_total / params['mass']
    
    # 转动方程（在箭体系中）
    I = params['I']  # 惯量张量
    M_body = np.zeros(3)  # 力矩（简化）
    domega_body = np.linalg.inv(I) @ (M_body - np.cross(omega, I @ omega))
    
    # 四元数导数
    dquat = 0.5 * quaternion_multiply(quat, np.hstack([0, omega]))
    
    return np.hstack([vel, acc, dquat, domega_body])

# 使用示例
params = {
    'mass': 1000,    # kg
    'mu': 3.986e14,  # 地球引力常数
    'thrust': 5000,  # N
    'Cd': 0.3,       # 阻力系数
    'S': 0.5,        # 参考面积
    'rho': 1.2,      # 大气密度
    'I': np.diag([100, 200, 200])  # 惯量张量
}

# 初始状态
state0 = np.zeros(13)
state0[2] = 1000  # 初始高度1000m
state0[6] = 1     # 初始四元数（无旋转）

# 积分
sol = solve_ivp(rocket_dynamics, [0, 100], state0, 
                args=(params,), method='RK45', max_step=0.1)
个人建议：刚开始做仿真时，先用简单的质点模型（不考虑转动）跑通流程。等确认质心运动没问题了，再加入转动方程。一步到位容易出bug，而且很难排查。

好了，坐标系和运动学就讲到这里。这些是火箭制导的「地基」，地基打不牢，后面盖什么楼都会塌。下一节我们会进入真正的制导算法，到时候这些坐标系变换会反复用到，务必熟练掌握。

公众号:蓝海资料掘金营，微信deep3321

← 上一章 📖 返回目录 下一章 →