卫星"失踪"83 分钟——两位卡尔曼侦探的破案实录（EKF × UKF）
本文是《航天动力学中的数学方法》系列第九篇，对应第 8 章"状态估计与滤波（EKF/UKF）"。配套代码见
chapter8/。
导语： 凌晨两点，测控大厅。
一颗 400 公里高的低轨卫星正在过境，值班工程师盯着屏幕：雷达刚测回来的距离，和轨道数据库里的预测，对不上。更麻烦的是——数据库里的初始位置本身就是"错的"：位置偏了 2.7 公里，速度偏了 0.7 米/秒。
"卫星没有失联，"总师喝了一口咖啡，"失联的是我们对它位置的确信。"
接下来的 83 分钟（5000 秒），每 10 秒会有一帧测距数据传回来。数据只有一个数：卫星到地面站的距离——而且每次测量都带着 ±100 米的噪声，像 400 公里外一声带静电的口供。
任务：从这些"不太可靠的证词"里，把她找回来。
在测控大厅的这 83 分钟里，破案的其实是两位侦探：EKF 与 UKF。今天我们把这份卷宗完整摊开——不背公式，只讲破案。
案卷 01 · 案情分析：6 个未知数，1 条线索
先看看这案子为什么"难办"。
卫星此刻的状态，数学上就是 6 个数字：位置三分量 (x, y, z) + 速度三分量 (vx, vy, vz)。而雷达每次交给你的，只有 1 个数字：距离 ρ。
6 个未知数、1 条线索、单帧无解。 先验信息里，位置的偏差本身高达 2.7 公里，速度的偏差也有 0.7 米/秒——这些偏差会像复利一样随时间滚雪球。
那怎么破？靠两条腿走路：
- 证词（观测）：雷达测距，虽然有噪声，但它至少"看着"卫星；
- 常理（动力学模型）：卫星不会瞬移、不会无故拐弯。给她一个状态，用二体 + J2 摄动的运动方程往前推 10 秒，你能预测她"应该"在哪。
破案流程就是"两条腿交替迈步"，也就是滤波器的预测—更新循环：
▼text复制代码用动力学往前推一步 ──────→ "根据常理，她现在应该在 A 点" （预测） │ 雷达送来一帧测距 ───────→ "证词说，她离站 400.6 公里" （观测） │ 两者对照，按可信度折中 ─→ 锁定新位置，进入下一个 10 秒 （更新）
每 10 秒重复一次，501 帧之后，真相就会从噪声里浮出来。
侦破手记 ▎ 单帧欠定、噪声污染、初始线索还是错的——这案子听着耳熟？对，这就是轨道确定（Orbit Determination）。从阿波罗导航到现在的星历维护，航天工程每天都在破这种案。
案卷 02 · 侦探一号 EKF：把曲线掰直了推理
第一位侦探叫 EKF（扩展卡尔曼滤波），侦探界的老派名宿。他的破案哲学六个字：以直代曲，局部推演。
卫星的运动是弯的（非线性的），但如果你只关心她周围"一小步"范围内的行为，曲线可以用切线代替——就像地图上找不到一条弯曲的路，就用一小段直线描出它的走向。
图 1：两位侦探的手法对照。左边 EKF：在工作点画一条切线（橙色），用"直线逻辑"逼近曲线逻辑——代价是每次都要算一次导数（雅可比矩阵 J = df/dx）。右边 UKF：根本不掰直，直接在工作点周围"撒 13 个采样点"，让它们各自跑一遍动力学，回来统计。
EKF 的工作流拆开看只有两句话：
- 预测："用动力学推 10 秒——但我承认推演不完美，不确定性会变大一点。"（数学上：状态用 RK4 传播，协方差矩阵 P 同步膨胀、加上过程噪声 Q）
- 更新："证词来了。我该信自己多少、信雷达多少？"
第二条里的"该信谁"，就是卡尔曼滤波的灵魂——卡尔曼增益 K。它是一台"话事权分配器"：
- 如果自己的预测信心很足（P 小）、雷达很粗（R 大）→ K 小，以我为主，证词听听就好；
- 如果自己刚开案、啥都不知道（P 巨大）、雷达很准（R 小）→ K 大，先按证词把人拉过来再说。
用代码说话（chapter8/ekf.py，核心就三行）：
▼python复制代码K = self.P @ H.T @ np.linalg.inv(H @ self.P @ H.T + self.R) # 话事权分配 self.x = self.x + K @ (z - hx) # 按分配执行修正 self.P = (I - K @ H) @ self.P @ (I - K @ H).T + K @ self.R @ K.T # 更新信心（Joseph 形式）
侦破手记 ▎ 注意最后一行用的是 Joseph 形式——工程老手都知道，协方差矩阵"算崩了"（失去正定性）是卡尔曼滤波最常见的翻车方式，这种写法能把它牢牢摁在安全区。另外 EKF 的雅可比是数值差分算出来的（给每个状态分量各戳一下，看输出怎么变），教学友好，工程上通常换解析公式省算力。
代价也说清楚："以直代曲"在小步长下够用，但轨道一弯得厉害（强非线性），直线逻辑就会开始骗人。 所以我们需要第二位侦探。
案卷 03 · 侦探二号 UKF：撒一把点去试探
UKF（无迹卡尔曼滤波） 是新派侦探，手法完全不同——不掰直，直接采样。
他的逻辑很"笨"也很聪明：既然我不确定她精确在哪，那我围绕"最有可能的位置"撒下 2n+1 = 13 个采样点（n = 6 是状态维数，这些点叫 Sigma 点）：中心放一个，剩下 12 个沿 6 个方向对称成对撒开，撒多远由当前的不确定性 P 决定——P 越大，网撒得越开。
然后：13 个点全部扔进 RK4 动力学各自跑 10 秒，回来后统计它们的均值（新的"最可能位置"）和散布（新的不确定性）。
这套逻辑的卖点：
- 免求导：不用算雅可比，只要有"状态推演器"就能干活；
- 不怕弯：13 个点分布在曲线附近，动力学怎么弯，它们就怎么弯着走，不需要"掰直"这一步；
- 代价：算量变成 13 倍（每次要传播 13 条轨迹），但换来的是对强非线性场景的天然适应。
代码里最传神的一小段（chapter8/ukf.py）：
▼python复制代码sqrtP = self._cholesky_psd((self.n + self.lmbda) * self.P) # 用 Cholesky 分解决定撒点方向与距离 for i in range(n): sigma[i + 1] = self.x + sqrtP[:, i] # 沿第 i 个方向撒远点 sigma[i + 1 + n] = self.x - sqrtP[:, i] # 对称位置撒近点
侦破手记 ▎
_cholesky_psd是个有故事的函数：协方差矩阵在数值误差下偶尔会"轻微不老实"（失去正定），这时它先用特征值截断把矩阵修好再分解——工程代码的体面，就是所有崩溃点都提前垫好了垫子。
案卷 04 · 破案实录：四件物证
好，两位侦探同时进场，每 10 秒一轮"预测—更新"，全程 501 轮。测控中心的四条记录曲线，就是本案的四件物证：
物证一：位置误差——"74 公里的炸雷"与它的平息
图 2：位置估计误差全程记录（纵轴单位：米）。
读图要点：
- 开局第 1 帧：UKF 抛出一颗"炸雷"——误差跳到 7.4 万米（橙线尖峰）。别慌，下一帧就被观测拽回到 800 米。EKF（蓝线）则以 2.6 公里平稳起步；
- 一段"过山车"：从 800 米的谷底一路爬升到约 23 分钟处的主峰 4.3 万米，随后长下坡，途中在 40~50 分钟区间还停了一段 1.3 万米左右的"高原"——这是"带噪声证词 + 波浪式可观测性"的代价。单站测距只对视线方向敏感，切向分量要靠时间积累 + 动力学慢慢分辨，收敛因此呈台阶状，而不是一条平滑的下坡；
- 约 4500 秒（75 分钟）：误差首次跌破 1 公里；末段停在数百米（本次记录约 900 米）；
- 两条曲线：此后的 2500 秒里，EKF 与 UKF 几乎逐帧重合——手法完全不同，结论殊途同归。
诚实说明（技术侦探的操守）： 那颗 74 公里"炸雷"是教学版实现的偶发数值过冲——巨大的初始协方差、极端参数与数值差分在特定条件下叠加所致，下一帧即被修正。生产级实现会用解析导数、稳健的 UKF 参数与平方根滤波来规避。它不影响本案主结论；保留它，因为它恰好展示了滤波器"自我纠错"的真实一面。
物证二：速度误差——从 164 m/s 到散步的速度
图 3：速度估计误差全程记录（纵轴单位：米/秒）。
- 首帧更新后速度误差一度冲到 87 m/s，约 200 秒处爬到峰值 164 m/s——这是"开局 P 巨大、滤波器太敢动"的余波；
- 之后曲线一泻千里：15 分钟回到 40 m/s 量级，末段稳定在 ~1 m/s；
- 什么概念？这颗卫星以 7.6 km/s 飞行（约 2.7 万公里/小时），而滤波器最终把它的速度误差锁到了约 1 米/秒——相当于一架超音速巡航的飞机，仪表读数精确到"人散步的速度"。
物证三：3σ 包络——侦探的自信刻度
图 4：实线是真实误差，虚线是滤波器"自我申报"的 3σ 不确定度包络。
- 开场虚线约 2000 米，随后随误差一起涨到 1.5 万米，末段回落到 约 7000 米；
- 关键判据：虚线始终罩住实线——这叫道"一致性（consistency）"。翻译成人话：这位侦探从没吹牛，他说"误差不会超过 7 公里"，实际误差就老老实实待在包络里面；
- 有趣的是末段：申报 3σ ≈ 7 公里，实际误差只有 0.9 公里——滤波器比实际更"谦虚"了 8 倍。这不是笨，是工程纪律：宁可低估自己的精度，绝不吹牛——自信过头的滤波器，才是测控事故的起点。
物证四：新息——证词与推理的"落差"
图 5：观测新息 z − h(x)——"雷达说距离是 A，我推演的距离是 B"，两者之差。红线为 ±3σ 观测噪声带（±300 m）。
- 开局是一记 −1.3 万米的深谷：模型和"证词"彻底谈不拢——因为初始线索本来就错了 2.7 公里，第一帧落差自然巨大；
- 随后震荡收拢：约 15 分钟内基本贴合进 ±300 米（±3σ）的噪声带，此后偶尔越界、但姿态已稳——"证词与推理的口径"对上了；
- EKF 与 UKF 的新息曲线几乎完全重叠——两位侦探对"落差"的处置，殊途同归。
侦破手记 ▎ 工程上，新息序列是滤波器的"体检报告"：如果它持续超带（系统疯跑）、或像死水一样毫无波动（观测被无视了），就说明滤波器某一环出了问题。测控大厅里盯新息，就像医生盯心电图。
结案报告
本案卷宗摘要：
| 物证 | 开局 | 末段（83 分钟后） |
|---|---|---|
| 位置误差（EKF / UKF） | 2.6 km / 74 km 尖峰 | ~0.9 km / ~0.9 km |
| 速度误差 | 87 → 164 m/s（峰值） | ~1 m/s 量级 |
| 3σ 包络 | ~2000 m | ~7000 m（保守覆盖，实测误差 0.9 km） |
| 新息 | 一记 −1.3 万米深谷 | 收进 ±300 m 噪声带 |
两位侦探的选型对照：
| 维度 | EKF | UKF |
|---|---|---|
| 核心手法 | 线性化：以直代曲 | Sigma 点采样：撒网统计 |
| 需要的"情报" | 雅可比矩阵（导数） | 只要有状态推演器 |
| 首帧 | 平稳起步（2.6 km） | 一次 74 km 暴走后立即回落 |
| 强非线性 | 依赖线性假设的"脸色" | 天然适应 |
| 算量 | 基准 | ×13（每次传播 13 条轨迹） |
| 本案例终局 | 0.9 km / 1 m/s 量级 | 0.9 km / 1 m/s 量级 |
真相公式，全案浓缩成一句：
真相 = 模型递出的"预期" ⊕ 观测带来的"修正"，按双方的可信度加权。
信心足的一方说话大声（增益 K 大），信心虚的一方自动闭嘴——这套"按可信度加权"的哲学，就是卡尔曼滤波 1960 年论文里的全部秘密。
隐藏彩蛋 1： 阿波罗登月舱的导航计算机（AGC）跑的就是卡尔曼滤波的变体——1969 年人类第一次登月时，72KB 内存里装着和本文同一套哲学。
隐藏彩蛋 2： 你手机里的 GPS 定位、SpaceX 的卫星自主避碰、天问/EHT 深空网的轨道维持，背后站着的都是这两位侦探的徒子徒孙。
如何运行本文代码？
▼bash复制代码pip install numpy matplotlib scipy python chapter8/plots_ch8.py # 生成四张数据图（本文物证 1~4） python chapter8/plots_ch8_article.py # 生成两张概念插图（封面 + 方法论对比）
案件重演：五组动手实验
欢迎在 chapter8/ 里改参数，重开调查：
- 换一把更准的"证词枪"：把
plots_ch8.py里sigma_r = 100改成10——雷达精度提升 10 倍，收敛时间会缩短多少？ - 改开局线索：把初始误差
np.array([2000, -1500, 1000, ...])加大到 10 倍，看重启后首帧"炸雷"会不会更响、滤波器还能不能自己爬回来； - 调过程噪声 Q：把 Q 放大 10 倍，观察 3σ 包络怎么变——你会亲手体验"信心预算"的调节杆；
- 多请一位证人：
chapter8/observation_model.py里现成备着方位角+仰角观测obs_az_el()——把角度证词加进来，俯瞰"多源观测"如何让侦破提速； - 拆掉 UKF 的拐杖：把
alpha=1e-3改成1e-1，看撒点尺度变化对手表精度和数值稳定性的影响。
