我在做火箭仿真系列的稳定性与动力学分析模块时，遇到最多的一个误解是：仿真里已经把轨迹画出来了，甚至姿态也在动，就默认动力学没问题。但真正把模块拆开看，轨迹只是"结果"，而稳定性与动态响应特性才是判断一枚火箭设计是否可靠的关键。这个模块的核心作用是回答两件事：给定初始扰动后，火箭的姿态会不会回到参考状态？如果会，收敛速度、振荡幅度、阻尼大小是否在设计包线内。
这篇文章就围绕我自己实现这套模块的过程展开，包括动力学建模、稳定性判据、数值线性化、时域响应分析与实际踩坑记录。如果你在做飞行器仿真、小型探空火箭项目，或者想把一个"能出轨迹"的仿真脚本升级成"能指导设计"的分析工具，那这部分内容应该对你有直接帮助。
1. 为什么火箭仿真里必须单独拆出一个稳定性与动力学模块
1.1 只有弹道曲线的仿真，说服力远远不够
我最早做的火箭仿真非常简单：给一个推力曲线，加上简化的阻力公式，积分出一条高度时间曲线。看起来像模像样，但换个尾翼尺寸、改一下质心位置，轨迹图的差别往往很小，你根本看不出哪个方案更稳。真正决定一枚火箭是"稳定飞行"还是"边飞边翻"的，是姿态扰动导致的力矩响应，而这些在只关注坐标位置的仿真里是隐性问题。
稳定性与动力学分析模块就是把这些隐性问题显式化。它不是替代轨迹计算，而是把力和力矩的平衡关系、姿态扰动的收敛性单独抽出来，形成一套可量化、可对比的指标。比如等效静稳定裕度、俯仰振荡阻尼比、短周期模态频率等。有了这些指标，你就可以在试飞之前判断"这发弹加长了翼展之后，抗风能力是否变差"，而不是等落地后看残骸。
1.2 模块在仿真系列里的定位：承上启下、可插拔
在我们这个仿真系列里，稳定性与动力学分析模块的位置很明确：上面接气动参数模块，下面接六自由度运动仿真内核，旁边还能挂控制系统模块用于闭环分析。这么做的好处是，气动数据、质量特性、推力曲线都可以独立替换，不需要改动稳定性核心代码，就可以对比"不同尾翼方案下的稳定性差异"或"推进剂消耗对动态特性的影响"。
模块的输入通常包括：
- 气动系数表：升力系数斜率、阻力系数、俯仰力矩系数斜率、俯仰阻尼导数；
- 质量特性：随时间变化的质心位置、转动惯量；
- 推力曲线与推力偏心参数；
- 初始状态偏差范围，用于设置扰动条件。
输出则是：
- 静稳定裕度随时间的变化曲线；
- 扰动运动特征值，即模态频率与阻尼比；
- 攻角、角速度的时域响应序列；
- 判定结论，是否满足预设的稳定边界。
把输入输出边界定清楚之后，这个模块就可以独立运行，也可以嵌入整个仿真系列做批量化方案对比。我在项目里通常还会让模块自动生成一份稳定性报告，包含图像和指标表格，方便团队评审时直接看结果。
2. 动力学模型：从力与力矩到六自由度运动的推导
2.1 坐标系约定：先统一坐标系，再谈公式
动力学建模的第一件事不是写方程，而是定义坐标系。火箭仿真里常用三套坐标系：地面惯性坐标系，用于描述位置和重力方向；火箭体轴坐标系，原点在质心，x轴沿弹体轴线指向头部，用来表达推力方向和力矩；速度坐标系，x轴沿速度方向，用来定义攻角和侧滑角。
我习惯采用3-2-1欧拉角顺序描述姿态，即先偏航、再俯仰、最后滚转。这样在攻角不大的情况下，姿态角与气动角的关系最直观。这里容易忽略的一个细节是：气动力和力矩的计算应该在速度坐标系或风轴系里进行，然后通过转换矩阵投影到体轴系。如果直接在体轴系里套用气动系数，侧滑角交叉耦合会被掩盖，尤其在滚转速率较大时误差很明显。
姿态动力学用四元数积分角运动比欧拉角更稳，不会出现万向节锁死。我最终的代码里，位置速度在前半部分用惯性系，姿态用四元数，角速度用体轴系，转换环节全部显式写成函数并做了单测。
2.2 力与力矩模型的细化程度
动力学模块里的力主要考虑三部分：
- 重力，简化取标准重力常值，按惯性系y轴方向施力即可。对短程火箭来说，重力随高度变化的影响很小，我在第一阶段没有加入；
- 推力，默认沿体轴线的轴向方向，同时增加一个推力偏斜角参数，用来模拟发动机装配误差。推力偏心会直接激励横向角运动，是稳定性分析里重要的扰动来源之一；
- 气动力，包括阻力、升力和侧向力。表达式统一写成 (F = C \cdot q \cdot S_{ref})，其中 (q) 是动压，(S_{ref}) 是参考面积。升力系数在小攻角范围内按 (C_L = C_{L\alpha} \cdot \alpha) 线性化。
力矩部分比力更关键。俯仰力矩写作：
[ M_z = C_{m\alpha} \cdot \alpha + C_{mq} \cdot q + C_{m\delta} \cdot \delta ]
其中 (C_{m\alpha}) 是静稳定力矩导数，(C_{mq}) 是俯仰阻尼导数，(C_{m\delta}) 是控制舵偏或燃气舵产生的控制力矩导数。在无控火箭分析里，第三项直接置零，但模块里还是预留了接口，方便以后接控制系统做闭环分析。
2.3 转动惯量变化：不能当常值处理的坑
固体火箭推进剂燃烧导致质量、质心和转动惯量都在变化，这一项最容易在大规模仿真中被忽略。把转动惯量固定成初始值，短期内轨迹差别不大，但稳定性特征值会明显飘移，因为转动惯量直接出现在姿态动力学方程的分母上，惯量变化会影响振荡频率和阻尼比。
我的处理方式是做一张质量表格，包含时间、质量、质心位置、俯仰转动惯量、偏航转动惯量和滚转动惯量。仿真积分每走一步，根据当前时间查表插值更新这些参数。发射初期推进剂刚点燃时质量变化最快，时间步长需要适当加密，否则插值出来的惯量值会有轻微跳变，造成响应曲线上的毛刺。
3. 稳定性判据：不是"翻不翻"，而是"扰动后回不回来"
3.1 静稳定性与等效静稳定裕度
静稳定性描述的是火箭在受到攻角扰动时，气动力矩是否有趋势把它推回平衡姿态。判据是俯仰力矩系数对攻角的导数 (C_{m\alpha} < 0)。飞行力学里更常说的是等效静稳定裕度：
[ S_m = \frac{x_{cp} - x_{cg}}{L_{ref}} ]
其中 (x_{cp}) 是压心位置，(x_{cg}) 是质心位置，(L_{ref}) 是参考长度。静稳定裕度越大，火箭越"僵硬"，抗扰动能力越强，但操纵性也越差。对无控火箭而言，一般希望它在整个飞行段内都保持正裕度，工程上常见取值在 0.1 到 0.3 之间。
火箭和飞机的最大区别在于：飞机重心基本不变，而火箭发射中质心持续前移，压心也会随马赫数变化。所以我不会只看某一时刻的静稳定裕度，而是画出整个主动段剩余时间的裕度曲线，检查是否存在"最小静稳定裕度时间点"，这个点往往就在推进剂快烧尽、质心后移幅值最大的时候。
3.2 动态稳定性：扰动模态、阻尼比与频率
静稳定只保证"有力往回推"，但推过头就会振荡，振荡能不能衰减要看阻尼。把扰动运动方程在平衡点附近线性化，可以得到典型的二阶系统形式：
[ \ddot{\theta} + 2\zeta\omega_n\dot{\theta} + \omega_n^2\theta = 0 ]
其中 (\omega_n) 是短周期模态自然频率，(\zeta) 是阻尼比。阻尼比为负时，扰动不仅不衰减，反而增大，称为动不稳定。阻尼比在 0 到 1 之间时，振荡逐渐收敛，但会伴随超调。阻尼比偏小是很多火箭试飞中姿态振荡衰减慢、发散风险高的根源。
对无控旋转稳定火箭，动稳定性分析还要考虑滚转速率与俯仰偏航模态的耦合，会出现一个章动模态和一个进动模态。这个模块的通用版本里我暂时只做刚体的短周期模态提取，把章动分析留到后续转动动力学版本。分析的逻辑不变：线性化、提取特征值、看实部符号和虚部频率。
3.3 稳定性分析的三种实现方式
我在模块里同时保留了三种稳定性检查方法，因为不同阶段的仿真用得不一样：
- 解析近似法：直接用简化公式估算静稳定裕度和短周期模态频率，适合设计前期快速筛选方案；
- 数值线性化法：对六自由度状态方程施加微小扰动，用有限差分构造雅可比矩阵，再求特征值。适合中段详细分析；
- 时域仿真法：给初始状态设置一个脉冲扰动，积分数十秒，观察攻角响应包络。这最接近真实试飞场景，适合最终验证。
三种方法的结果应该有相当的一致性。如果数值线性化得到的阻尼比和时域响应曲线提取到的对数衰减率对不上，那说明要么线性化步长不合适，要么时域数据滤波破坏了模态信息。这两类问题我后面会专门展开讲。
4. 模块实现：状态方程、线性化与响应分析的代码落法
4.1 状态向量设计与积分器选择
状态向量我用了 13 维：惯性系位置 3 个、惯性系速度 3 个、四元数 4 个、体轴系角速度 3 个。之所以不用欧拉角，除了避免万向节锁死，另一个原因是数值线性化时，四元数约束下的扰动处理虽然麻烦一点，但数值行为更干净。
积分器选择了四阶龙格库塔方法，固定步长 0.01 秒。固定步长的好处是线性化时状态轨迹更容易复现，自适应步长在比较响应曲线时会引入额外的插值误差。整个主动段也就几十秒，计算量不大，稳定性和效率都能兼顾。
状态方程的右端函数集中在一个函数里，接收当前状态和一组配置参数，返回状态导数。配置参数包括气动系数表、质量插值表、推力曲线、参考面积和参考长度。这样做的目的是让数值线性化只用调用同一个右端函数就能构造雅可比矩阵。
4.2 数值线性化：有限差分雅可比
数值线性化的思路不复杂：把状态方程的右端函数对每个状态变量求偏导。先找一个基准状态，通常取初始平衡状态或某一特征时刻的状态，然后对每个变量施加一个微小扰动 (\epsilon)，计算 ((f(x+\epsilon)-f(x-\epsilon))/(2\epsilon))，得到雅可比矩阵。
步长的选择很关键。我在调试中遇到过两个极端：步长取 1e-5 时，因为浮点舍入误差明显，特征值虚部出现抖动；步长取 1e-3 时，动力学函数的非线性结构又不够精确，阻尼比误差偏大。最终通过扫描 1e-4 到 1e-8 范围，选定 1e-6 作为默认步长。以下是一段用于求取雅可比矩阵的示意代码，核心逻辑可以直接复用。
import numpy as np def compute_jacobian(deriv_func, state, cfg, eps=1e-6): n = len(state) jac = np.zeros((n, n)) for i in range(n): step = np.zeros(n) step[i] = eps jac[:, i] = (deriv_func(state + step, cfg) - deriv_func(state - step, cfg)) / (2 * eps) return jac def extract_modes(jac): # 注意只保留姿态相关子空间更稳，我这里取完整矩阵的特征值 eigvals = np.linalg.eigvals(jac) return eigvals提取特征值之后，关注的是复共轭对。实部 (\sigma) 和虚部 (\omega) 换算成阻尼比 (\zeta = -\sigma / \sqrt{\sigma^2 + \omega^2})，自然频率 (\omega_n = \sqrt{\sigma^2 + \omega^2})。如果出现实部为正的特征值，说明这个运动模态不收敛，必须检查对应的扰动方向。
4.3 时域脉冲扰动与响应指标提取
时域响应分析更直观：在仿真初始时刻给俯仰角一个设定好的偏差，比如 3 度攻角，然后让模块自由飞行仿真，记录攻角随时间的变化序列。稳定的火箭会呈现出衰减振荡，我用攻角响应相邻峰值的比率估算对数衰减率，再换算阻尼比和时间常数。
提取峰值时不要直接把原始数据里的每个局部极大值都当真。仿真步长固定时，峰值点可能落在两个积分节点之间，直接取原始序列最大值的误差在采样率不足时会比较大。我用的办法是先对攻角时域序列做样条插值，再用插值函数找峰值，得到的阻尼比结果与特征值法之间的偏差能控制在 3% 以内。
5. 实测案例：一组典型仿真结果解读
5.1 案例参数与仿真设置
为了说清楚这个模块怎么用，我拿项目中做过的一发小型试验火箭参数举例，参数经过了脱敏处理，只保留量级作为参考。该弹采用尾翼稳定无控方案，总长约 2.2 米，起飞质量 48 公斤，推进剂占 16 公斤。
| 参数 | 数值 | 备注 |
|---|---|---|
| 参考长度 L_ref | 0.22 m | 取弹体直径 |
| 参考面积 S_ref | 0.038 m² | 按直径计算 |
| 初始质量 | 48 kg | 含推进剂 |
| 剩余质量 | 32 kg | 推进剂燃尽后 |
| 初始质心位置 | 1.34 m | 距头部 |
| 压心位置范围 | 1.28 m ~ 1.51 m | 随马赫数变化 |
| C_m_alpha | -0.075 /deg | 线性段近似 |
| C_m_q | -16.0 /rad/s | 俯仰阻尼导数 |
仿真从零时刻点火算起，初始攻角扰动设为 3 度，同时给俯仰角速度一个 0.5 deg/s 的初始偏差，用来观察阻尼项的作用。推力曲线按时间插值，取燃烧时间 8 秒，最大推力 4800 牛。
5.2 仿真结果：静稳定裕度与动态阻尼
首先看静稳定裕度，起飞瞬间约 0.21，随推进剂消耗和压心前移，在 6.8 秒附近裕度降到极小值 0.13，随后缓慢回升。这个谷值意味着扰动最不容易收敛的时段出现在接近助推末端，设计上是可以接受的，如果需要加固，尾部翼面积还能再增大。
动态模态的提取结果，短暂起飞的 1.0 秒时刻，数值线性化得到一对共轭特征值，实部约为 -0.52，虚部约 2.04，阻尼比约 0.25，自然频率约 2.1 rad/s。时域响应对数衰减率换算出的阻尼比是 0.24，两者一致。攻角响应曲线头两个振荡周期内，峰值从 3.0 度降到约 1.9 度，衰减趋势明显。
这个阻尼比在无控火箭里算中等水平，抗风扰动能力尚可，但余量不算大。如果蒙特卡洛仿真显示发射初段侧风 8 m/s 条件下最大攻角接近 6 度，就要考虑增大阻尼项或者提高静稳定裕度。从响应曲线上看，振荡周期约 3 秒，对应自然频率 2.1 rad/s 的物理意义很直观：姿态完成一次来回摆动大约 3 秒，这在试飞遥测数据里是可以直接看到的。
5.3 从结果反推设计改进
稳定性分析模块最有价值的用途是根据不满足要求的指标做设计微调。如果阻尼比偏低，最直接的手段是增大尾翼面积，等效压心后移，因而静稳定裕度增加，短周期频率升高，振荡更快。但代价是气动阻力增大，航程缩短，且滚转与偏航耦合变强。另一个措施是增加翼面展弦比来提升阻尼导数，这在保持压心位置变化不大的情况下改善阻尼比，效果更定向。
我在实际项目里通常先用模块做参数扫描：把尾翼面积从基准值乘 0.9 到 1.3 各跑一遍，输出每组参数对应的稳定裕度、阻尼比和超调量，形成类似设计曲线图。这样比凭感觉调参数可靠得多，因为火箭的压心随马赫数变化是非线性的，面积增大的收益在高速段可能明显缩小。
6. 踩坑记录：稳定性模块最容易出错的地方
6.1 转动惯量不更新，特征值漂到你怀疑人生
有一版仿真结果让我一度以为气动数据错了。特征值在某个时间窗口内自然频率突然从 2.1 掉到 1.6，阻尼比也忽高忽低。排了三天，最后发现是转动惯量没有按时间插值更新，积分器一直用初始状态的 (I_{yy}) 值。推进剂燃烧导致 (I_{yy}) 从 5.8 kg·m² 降到 4.2 kg·m²，变化超过 27%，模态频率自然显著改变。
从那以后，我在模块里加了一个断言：质量特性插值表的时间覆盖范围必须覆盖整个仿真时长，否则直接报错。这样虽然粗暴，但能避免静默错误，值得推荐给同样在做仿真的朋友。
6.2 数值线性化步长的选择不是小事
有限差分雅可比矩阵的步长选择，直接影响特征值提取。步长过大，非线性项被线性化到配置之外的区域；步长过小，浮点有限精度又把信号淹没在舍入噪声里。我建议的做法是先用一个典型的基准状态，扫描步长从 1e-3 到 1e-9，选出特征值实部、虚部都稳定不波动的区间。
我在模块里没有把扫描写进每个任务，而是把默认步长 1e-6 和"复算校验"做成标准流程。每次分析都额外跑一遍步长加倍和减半的特征值，三者差异超过 1% 就提示步长异常。这个方法简单有效，能挡掉很大一部分线性化引起的假发散结论。
6.3 数据滤波：别把物理振荡当成噪声滤掉
时域响应数据做平滑处理是常事，但一不小心就把真实的短周期振荡给滤没了。我最初用滑动平均窗口长度为 0.5 秒的输出做峰值提取，结果阻尼比被严重低估，因为每个振荡周期的波峰都被平滑抹掉了。
正确做法是平滑窗口取振荡周期的十分之一以内，或者干脆不对原始响应做平滑。相比平滑，更推荐先确认弹簧质量的响应序列没有异常毛刺，然后用样条插值精确找峰。如果遥测数据里有高频噪声，再用截止频率高于目标模态频率 5 倍以上的低通滤波器，才能保留模态信息。
6.4 模块接口与回归测试的组织心得
最后说一点代码组织经验。稳定性分析模块如果和轨迹仿真混在同一个文件里，后期换数据、加功能会非常痛苦。我最终把模块的输入全部配置化成字典或数据类，仿真内核只依赖配置对象，这样同一个分析代码可以在不同项目间复用，只需要换参数表。
同时，每次修改气动数据或质量表之后，我会跑一组固定基准案例的回归测试。基准案例的稳定性指标是历史版本计算过的，代码改动后如果指标对不上，说明改动引入了问题。这套回归流程说起来平淡无奇，但在实物飞行测试前的仿真阶段，它真的帮我避免过两次重大设计失误。
用这个稳定性与动力学分析模块跑完一发弹的完整评估流程后，我自己最深的体会是：火箭设计迭代里真正缺的不是"某一次仿真做得对不对"，而是一套能反复给出稳定指标的可靠工具。静稳定裕度、模态阻尼比、时域响应这几个可量化指标同时摆出来，设计方案之间的差异就变得一目了然。而后端模块化的写法也让后续接飞控、加蒙特卡洛分析变得顺理成章，这条工具链我们可以持续用很久。
