变分贝叶斯自适应卡尔曼滤波：解决量测噪声未知的工程实践与MATLAB实现

简介：这份资源提供了一套变分贝叶斯自适应卡尔曼滤波的MATLAB实现，面向具有一定数学与编程基础的研究人员、工程师及高校师生，主要解决目标跟踪、导航与控制系统中非线性动态模型的滤波精度与稳定性问题。资源包内含17个文件，总大小约216KB，以MATLAB脚本为主，包含算法主程序、非线性建模、参数学习、迭代优化、高斯混合模型、最小均方误差等模块，兼有说明文档与论文文本，便于对照原理理解实现细节。目前已有73人学习使用。通过该资源可掌握变分贝叶斯推断与卡尔曼滤波的融合思路，了解如何在环境变化时自适应调整噪声参数，并通过示例代码验证算法在非线性场景下的滤波效果。压缩包目录层次清晰，适合作为算法研究与课程设计的参考资料，建议结合运行步骤动手实践，以深入理解参数学习与自适应机制的工程实现。

1. 变分贝叶斯自适应卡尔曼滤波：当卡尔曼的「最优」假设失效时

卡尔曼滤波在很多工程师手里像个黑匣子——给一组初始值，调好过程噪声协方差 Q 和量测噪声协方差 R，就以为输出是「最优估计」了。但「最优」两个字是有前提的：噪声统计特性必须精确已知，且在整个运行过程中保持不变。现实里这个前提往往不成立，雷达测距误差受目标姿态影响、GPS 量测方差随卫星几何分布变化、IMU 的噪声特性随温度漂移。R 一旦给错或者中途变了，标准卡尔曼滤波的表现会迅速恶化：R 给得偏大，滤波增益被压低，动态响应迟钝，目标转弯时直接跟丢；R 给得偏小，量测噪声毫无衰减地穿透滤波器，估计曲线毛刺丛生。

变分贝叶斯自适应卡尔曼滤波（VB-AKF）解决的就是这个问题。它在标准卡尔曼更新的基础上加了一层迭代：把量测噪声协方差 R 当成一个实时变化的随机变量，用变分贝叶斯方法在滤波的同时逼近它的后验分布。R 变了它能跟上，R 没变它也不乱动。这篇笔记从原理讲到 MATLAB 实现，再把参数调法和踩坑记录一次说清，新手工可以照着复现，熟手能直接拿去对比自己的实现。

2. 为什么是变分贝叶斯：标准卡尔曼滤波的盲区与其他自适应方案的局限

2.1 标准卡尔曼滤波的「五个已知量」：失配是怎么发生的

标准卡尔曼滤波对线性高斯系统是最小方差意义下的最优估计器，但这个最优性建立在一组苛刻的假设上：状态转移矩阵 A、量测矩阵 H、过程噪声协方差 Q、量测噪声协方差 R、初始状态分布全都精确已知。实际工程里，A 和 H 通常可以从运动模型或传感器几何中拿到，Q 和 R 却几乎都是拍脑袋定的。更麻烦的是，R 往往不是常数。

拿一个典型场景举例：某巡检机器人用里程计配合激光测距做定位，激光测距方差在室内、室外、强反光墙面三种环境下可以差一个数量级。固定 R 的卡尔曼滤波只能取一个折中值，结果在任何一种环境下都不是最优。R 给大（比如按保守值 10 给，实际是 0.1），滤波器认为量测很不可信，长期依赖运动模型外推，位置误差随时间慢慢累积；R 给小（按 0.1 给，实际是 10），滤波器把每次量测的噪声当成真实运动，状态估计方差很大，甚至可能出现滤波发散。

这个问题的本质是：标准卡尔曼滤波把 R 当成已知常量，而实际系统里 R 是一个随时间缓慢或突然变化的未知量。解决思路有两个方向，一是想办法在线估计 R 的值，二是把 R 的不确定性直接建模进滤波框架里。变分贝叶斯走的是第二条路，但理解它之前，先看看第一条路上那些方案为什么不够用。

2.2 创新序列法、多模型法与变分贝叶斯法的边界对比

创新序列法（Innovation-based Adaptation）是最直观的在线估计方案：卡尔曼滤波每一步都会产生残差（innovation），理论残差协方差是 H P⁻ Hᵀ + R，实际残差协方差统计出来以后，反推 R 的数值。这个思路在 R 基本不变、数据量足够的场景下能工作，但它有两个硬伤：一是需要一段滑动窗口内的残差样本来做统计，窗口内 R 不能突变，否则统计结果毫无意义；二是残差协方差和 R 之间隔着一个 H P⁻ Hᵀ 项，P 本身又在变化，反推出来的 R 方差很大，经常出现估计值为负这样的荒唐结果。

多模型法（IMM、MMAE）的思路是把 R 的可能取值离散成一组假设，每个假设跑一个独立的卡尔曼滤波器，按后验概率加权输出。这个方法在 R 恰好落在假设集合里的场景表现很好，但在集合之外表现很糟。想让集合覆盖所有可能的 R，假设数量指数增长，计算量直接失控。工程上一般只敢设三五个模型，R 真值落在模型间隙里时，加权结果比固定 R 好不到哪里去。

变分贝叶斯法走的是第三条路：不对 R 做离散假设，而是用一个连续分布去逼近 R 的后验。它不假设 R 在一段时间内恒定，每个时刻都可以更新对 R 的认知，R 慢漂移能跟，R 突变也能在十几步内追上。计算量是每次滤波更新里多几次迭代，量级可控。这是它在这三个方案里最适合工程落地的根本原因。

2.3 VB-AKF 的核心建模：Inv-Gamma 先验与定点迭代

变分贝叶斯自适应卡尔曼滤波的核心建模分三步。第一步，把量测噪声协方差 R 的对角元素建模为独立同分布的 Inverse-Gamma 分布。为什么选 Inv-Gamma 而不是高斯分布？因为协方差必须始终为正，高斯分布允许负值，而 Inv-Gamma 的支撑域就在正实数上。Inv-Gamma 有两个参数：形状参数 α 和尺度参数 β，其均值为 β/(α-1)，方差随 α 增大而减小。α 越大，说明对 R 的估计越确信。

第二步，用变分贝叶斯框架做后验逼近。真实后验 p(x, R | y) 很难直接求出，变分法的思路是找一个可分解的近似分布 q(x, R) ≈ q(x) q(R)，把依赖关系解耦，再用坐标上升法交替更新。固定 R 时更新 q(x)，走的就是标准卡尔曼更新的数学形式；固定 x 时更新 q(R)，可以得到 α 和 β 的解析更新公式。两步交替迭代若干次，状态估计和噪声统计一起收敛。

第三步是关键的时间更新，全靠遗忘因子 ρ 实现。每个时刻进入滤波之前，先把 α 和 β 都乘以 ρ（ρ 在 0.9 到 0.999 之间取值），相当于对过去所有时刻累积的噪声统计信息做指数衰减。ρ 越接近 1，历史信息保留得越多，R 的估计越稳但反应越慢；ρ 越小，旧信息丢得越快，R 突变时跟得越快但估计方差越大。这个参数是 VB-AKF 里最需要花心思调的，后面第 4 章单独展开。整体来看，VB-AKF 没有改变卡尔曼滤波的预测-更新结构，只是在更新步里嵌了一层定点迭代，工程上非常好落地。

3. 在 MATLAB 里跑通 VB-AKF：最小可复现例程与逐段说明

3.1 场景设计：一维匀速运动目标与量测方差突变

为了把 VB-AKF 和标准卡尔曼滤波的差距讲清楚，我构造一个最容易观察差异的仿真场景：一维匀速运动的物体，传感器只测位置，量测噪声方差在第 150 步从 0.1 突变到 10，相差 100 倍。这个场景刻意设计得极端，因为固定 R 的卡尔曼滤波在突变后会非常难看，VB-AKF 的自适应能力一眼就能看出来。

系统模型是常速运动模型，状态向量 x = [位置; 速度]，采样周期 dt = 0.1 秒，状态转移矩阵 A = [1 dt; 0 1]，量测矩阵 H = [1 0]，过程噪声协方差 Q 取 diag([0.01 0.01])，表示位置和速度上都有小幅随机扰动。总时长 300 步，突变点设在 150 步，留出足够的稳态段观察估计是否收敛。量测方差突变在传感器里对应的物理情形很多：激光测距仪从近距离转到远距离、雷达目标 RCS 突变的瞬间、GPS 从开阔地进入隧道口，都是同一类问题。

3.2 核心实现：仿真数据生成、标准 KF 与 VB-AKF 滤波循环

下面这段 MATLAB 脚本可以整段复制运行。为了让你看清差异，标准卡尔曼滤波和 VB-AKF 用同一个仿真数据，只在滤波环节分开写。

%% VB-AKF 与标准 KF 对比仿真 % 场景：一维匀速运动目标，位置量测，量测方差在第150步从0.1突变到10 clear; clc; rng(42); % ---------- 系统参数 ---------- dt = 0.1; % 采样周期 A = [1 dt; 0 1]; % 状态转移矩阵 H = [1 0]; % 量测矩阵，只测位置 Q = diag([0.01 0.01]); % 过程噪声协方差 T = 300; % 总步数 % ---------- 生成真实轨迹与量测 ---------- x_true = zeros(2, T); % 真实状态 [位置; 速度] y_meas = zeros(1, T); % 量测序列 R_true = 0.1 * ones(1, T); % 真实量测方差序列 R_true(150:end) = 10; % 第150步起方差突变为10 x_true(:,1) = [0; 1]; % 初始位置0，初始速度1 for k = 1:T if k > 1 w = mvnrnd([0 0], Q)'; % 过程噪声 x_true(:,k) = A * x_true(:,k-1) + w; end v = sqrt(R_true(k)) * randn; % 量测噪声 y_meas(k) = H * x_true(:,k) + v; end % ---------- 标准卡尔曼滤波（假设 R 恒为 0.1） ---------- x_kf = zeros(2, T); P_kf = eye(2); R_fixed = 0.1; % 错误但固定的 R 假设 for k = 1:T if k == 1 x_pred = x_true(:,1); P_pred = eye(2); else x_pred = A * x_kf(:,k-1); P_pred = A * P_kf * A' + Q; end S = H * P_pred * H' + R_fixed; K = P_pred * H' / S; % 用右除代替 inv，数值更稳 x_kf(:,k) = x_pred + K * (y_meas(k) - H * x_pred); P_kf = (eye(2) - K * H) * P_pred; end % ---------- VB-AKF 滤波 ---------- rho = 0.95; % 遗忘因子，需重点调参 Niter = 10; % 变分迭代次数 R0 = 0.1; % R 的先验均值初值 alpha0 = 1.5; % Inv-Gamma 形状参数初值 beta0 = R0 * (alpha0 - 1); % Inv-Gamma 尺度参数，保证均值恰为 R0 x_vb = zeros(2, T); P_vb = eye(2); alpha = alpha0; beta = beta0; R_hist = zeros(1, T); % 记录每个时刻的 R 估计 for k = 1:T if k == 1 x_pred = x_true(:,1); P_pred = eye(2); else x_pred = A * x_vb(:,k-1); P_pred = A * P_vb * A' + Q; end % ---- 变分时间更新：用遗忘因子衰减旧的噪声统计 ---- alpha = rho * alpha; beta = rho * beta; % ---- 变分定点迭代：交替更新状态与噪声参数 ---- for iter = 1:Niter R_hat = beta / (alpha - 1); % Inv-Gamma 均值 S = H * P_pred * H' + R_hat; K = P_pred * H' / S; x_upd = x_pred + K * (y_meas(k) - H * x_pred); P_upd = (eye(2) - K * H) * P_pred; innov = y_meas(k) - H * x_upd; % 更新后的残差 alpha = alpha + 0.5; beta = beta + 0.5 * (innov^2 + H * P_upd * H'); end x_vb(:,k) = x_upd; P_vb = P_upd; R_hist(k) = R_hat; end % ---------- 结果可视化 ---------- figure; subplot(2,1,1); plot(1:T, x_true(1,:), 'k', 'LineWidth', 2); hold on; plot(1:T, x_kf(1,:), 'b--', 'LineWidth', 1.2); plot(1:T, x_vb(1,:), 'r-', 'LineWidth', 1.2); legend('真实位置','标准KF','VB-AKF'); xlabel('时刻'); ylabel('位置'); subplot(2,1,2); semilogy(1:T, R_hist, 'r-', 'LineWidth', 1.5); hold on; semilogy(1:T, R_true, 'k--', 'LineWidth', 1.5); legend('VB-AKF估计的R','真实R'); xlabel('时刻'); ylabel('量测方差');

代码里的逻辑和参数说明值得逐段看。仿真数据生成部分，mvnrnd([0 0], Q)产生二维过程噪声，每次给状态加扰动；量测噪声的标准差是sqrt(R_true(k))，所以 R 突变后，量测点会立刻明显变散。标准卡尔曼滤波部分故意用固定R_fixed = 0.1，这个值在突变前是正确的，突变后就完全失效，用来充当对照组。注意这里用H * P_pred * H' / S而不是inv(S)，MATLAB 的右除运算内部走 QR 分解，数值稳定性比显式求逆好得多，这个习惯在滤波代码里值得贯彻到底。

VB-AKF 部分，alpha0 = 1.5和beta0 = R0 * (alpha0 - 1)的配合是有讲究的：Inv-Gamma 均值公式是 β/(α - 1)，这样设置让 R 的先验均值恰好落在 R0 = 0.1 上。如果直接把alpha0设成 1，均值公式分母为零，先验均值无定义，实现上就埋了一个隐患。迭代循环里，先用当前的 R_hat 做一次标准卡尔曼更新，得到 x_upd 和 P_upd，再算残差和投影协方差，把 α 和 β 各加 0.5，对应每个时刻量测带来的信息增量。Niter = 10 意味着这套动作重复十遍，状态和 R 的估计在迭代中互相校正、一起收敛。

3.3 结果怎么看：R 估计曲线与三段误差对比

脚本跑完后，重点看两个输出。第一个是 R_hat 的时间曲线：突变前它应该稳定在 0.1 附近，第 150 步之后开始爬升，在 ρ = 0.95 的衰减速度下，大约 20 到 30 步内逼近 10 附近，之后小幅波动。第二个是位置估计曲线：标准卡尔曼滤波在第 150 步之后明显跟丢，位置误差迅速增大且抖动剧烈；VB-AKF 在短暂的过渡段（大约 10-20 步）误差变大，之后重新收敛，稳态误差回到突变前量级。

如果 R_hat 曲线完全不跟着 R_true 走，先检查是不是 α 和 β 的初始化方向出了问题。很多人习惯直接beta0 = R0，但 α0 取 1 时均值无定义，R_hat 算出来是无穷大或 NaN，滤波直接翻车。另一个常见问题是 ρ 设得太接近 1（比如 0.999），第 150 步突变后 α 和 β 衰减极慢，R_hat 需要上百步才能爬到新值，看起来就像「没有自适应」。仿真里 ρ 取 0.95 就是为了在 300 步的时长内把突变响应看得清清楚楚。

4. 三个必调参数：遗忘因子 ρ、迭代次数 N 与初始先验 R0

4.1 遗忘因子 ρ：在「反应快」和「估计稳」之间找平衡

遗忘因子是 VB-AKF 里影响最大的参数，没有之一。它的物理含义是每个时刻对历史噪声统计信息的保留比例，ρ = 0.95 表示每个时刻旧信息打 95 折，10 步后旧信息只剩 0.95¹⁰ ≈ 0.6，半衰期大约 13 步；ρ = 0.99 的半衰期则是 ln(0.5)/ln(0.99) ≈ 69 步。这个半衰期直接决定了 R 突变后滤波器的追赶速度。R 真值在第 150 步突变 100 倍时，ρ = 0.95 的 VB-AKF 约 30 步能接近新值，ρ = 0.99 需要上百步，和固定 R 的表现差别不大。

ρ 不能为了「反应快」而无限调小。我见过有人把 ρ 调到 0.8，R_hat 曲线抖动得像噪声发生器，状态估计也跟着剧烈摆动。原因是 α 和 β 被最近一两步的残差主导，Inv-Gamma 分布的均值本身就有了很大的方差。实际选值规律是：量测噪声是突变型（进出隧道、目标 RCS 跳变），ρ 取 0.92 到 0.97，让滤波器在几十步内完成追赶；量测噪声是缓慢漂移型（热漂移、老化），ρ 取 0.98 到 0.999，优先保证稳态估计的平滑性。没有先验知识时，从 0.95 起步，看 R_hat 曲线的响应速度和波动幅度再微调。

4.2 变分迭代次数 N：5 次与 30 次的差别

VB 定点迭代的收敛速度比想象中快。第 1 次迭代把 R 从初始先验均值拉到接近当前残差水平，相当于完成粗估计；第 2 到第 5 次迭代用更新后的状态重新计算残差，修正 R 估计里的偏差；第 5 次之后每次迭代带来的变化已经很小，基本是数值层面的精细调整。我做过一个简单的对比实验：Niter 从 1 增加到 5，位置 RMSE 改善明显；从 5 增加到 15，改善幅度很小；从 15 增加到 30，差异基本淹没在随机噪声里。

从计算量角度看，VB-AKF 的单次迭代就是一次标量（或低维矩阵）求逆和几个矩阵乘法，10 次迭代的额外开销在 MATLAB 里通常连 1 毫秒都不到，在嵌入式平台上也就是多几个微秒量级。所以在计算资源允许的情况下，Niter 取 10 到 15 是稳的选择。你可以在仿真里把 Niter 设成 1、5、15 各跑一遍，对比 R_hat 曲线的收敛轨迹：Niter = 1 时 R_hat 有可见的振荡，Niter ≥ 5 时曲线基本光滑。

值得单独提醒的是，迭代次数和 ρ 有耦合关系。ρ 很小的时候，α 和 β 被快速衰减，即便 Niter 增加，R_hat 的稳态波动也不会明显改善；ρ 很大的时候，迭代次数对突变响应速度的影响也不大，因为瓶颈在遗忘速度。调参时先固定 Niter = 10，只调 ρ，等 ρ 定下来后再回头确认 Niter 是否够用。

4.3 初始先验 R0 的敏感性：别让 α 和 β 从离谱的地方起步

VB-AKF 对 R0 的敏感性比很多人想象的低，因为它有自修正能力——每步的残差会不断把 α 和 β 推向数据支持的方向。但初始值给得过于离谱会带来两个问题。第一个问题是收敛变慢：R0 比真值大两个数量级时，α 和 β 需要很长时间把均值拉回来，在这期间滤波器的增益被过大的 R_hat 压低，动态误差偏大。第二个问题是 ρ 会不断丢弃历史信息，如果收敛时间超过遗忘时间常数，R_hat 可能永远到不了真值附近，形成一个偏置的稳态。

实践中我一般这样给初值：如果对 R 完全没有概念，取系统启动后前 50 到 100 个量测样本，算一算残差的样本方差，把结果作为 R0。这个估计虽然粗糙，但数量级是对的，足够当 VB-AKF 的起点。α0 的取值也有讲究，前面代码里取了 1.5，实际在 1.01 到 10 之间都有人用。α0 越小，初始分布越宽，对 R 的初始约束越弱，自适应空间越大；α0 越大，初始分布越尖，R0 的分量越重。推荐从 α0 = 1.5 或 2.0 开始，这两个值能保证均值有限且分布足够宽，给数据留出充分的说话空间。

5. VB-AKF 避坑指南：5 个高频踩坑记录与排查方法

5.1 滤波中途发散，P 矩阵变成非正定

现象：程序跑了几十步后，状态估计突然跳变到巨大数值，P 矩阵出现负的特征值，后续出现 NaN。这个坑在 VB-AKF 里比标准卡尔曼滤波更容易触发，因为迭代循环里矩阵乘法次数更多，舍入误差累积更快。

原因：最直接的原因是 P 矩阵在反复迭代中失去了对称性，数值上出现微小不对称，在 Cholesky 或求逆操作里被放大；另一个推手是 R_hat 暂时性地被迭代推到极小的值（比如接近 0），导致 S 矩阵奇异或条件数过大。

解决：在每次迭代更新完 P_upd 后，加一句P_upd = (P_upd + P_upd') / 2，强制对称，成本极低但效果显著。同时给 R_hat 一个下限，比如R_hat = max(R_hat, 1e-6)，防止迭代过程中噪声协方差被压到零。最后检查一下 A 和 H 矩阵有没有病态，如果状态维度里同时存在位置和速度，且 dt 很小，P 的条件数会天然偏大，这时考虑用平方根滤波形式维护 P。

5.2 ρ 调小后 R_hat 剧烈抖动，状态估计跟着抖

现象：把 ρ 从 0.95 调到 0.85 想加快突变响应，结果 R_hat 曲线变成锯齿状，每步都在几个数量级之间跳，状态估计的误差反而比固定 R 还大。

原因：ρ 太小意味着每个时刻的 α 和 β 基本只由最近一两步的残差决定，而单步残差本身是随机量，它的平方波动很大。Inv-Gamma 均值对 β 的变化非常敏感，于是 R_hat 被单步噪音牵着走，估计方差极大。

解决：把 ρ 调回 0.92 以上，用半衰期公式对照你的场景需求。如果是突变型噪声，0.92 到 0.97 足够快；如果嫌慢，问题可能出在 Niter 太小而非 ρ 太小。另外，判断 R_hat 是否合理时要看趋势而不是单点值，画图时用semilogy看对数坐标下的整体走势，不要被个别尖峰误导。

5.3 R_hat 收敛到明显偏离真值的错误值

现象：R 真值是 10，VB-AKF 估计稳定在 3 到 4 附近，位置误差倒是没发散，但明显不是最优状态。更隐蔽的是，这种「稳定但错误」的估计比发散更难排查。

原因：这是 Q 与 R 的不可辨识性问题在作怪。VB-AKF 默认过程噪声协方差 Q 是已知且准确的，但如果你给的 Q 偏大，一部分量测误差会被模型解释成过程噪声，R 的估计就系统性地偏小。本质上，Q 和 R 都描述不确定性，量测误差到底归谁，光靠量测序列本身是分不清的。

解决：先确认 Q 的可信度。如果 Q 也是估计的，用一段真实数据做离线标定：固定 Q 跑 VB-AKF，看 R_hat 是否接近预期；把 Q 调小试一次，看 R_hat 是否上升。如果系统激励不充分，Q 和 R 的组合可能根本不唯一，这时任何自适应算法都有理论天花板，需要增加传感器或主动激励。如果 Q 和 R 真的同时未知，就得考虑扩展的变分方案，同时对两者建模，但参数更多、调参更难，工程上慎用。

5.4 多维量测下 R 估计发散，非对角项出现负值

现象：量测从一维扩展到二维或三维后，简单地把所有量测维度的 α、β 合并成向量来更新，结果 R_hat 对角线还好，非对角项直接发散，甚至有负值。

原因：VB-AKF 的经典形式是用独立的 Inv-Gamma 先验对 R 的每个对角元素建模，它从一开始就只能估计对角矩阵。如果直接把量测残差的外积（innov * innov'）整体当作更新量，非对角项的更新没有合法的先验约束，数值上就无法保证正定性。

解决：确认你的实现是逐维度独立更新 α 和 β，每个维度只用自己的标量残差和投影协方差，最终得到的 R_hat 是只含对角元素的对角矩阵。如果量测噪声确实存在明显相关性，需要换用 Inverse-Wishart 分布作为 R 的先验，那是另一套迭代公式和另一组参数更新规则，可以直接替换 VB-AKF 的「更新 αβ」步骤，但代码复杂度会上一个台阶。绝大多数导航和跟踪场景里，量测噪声关联性弱，对角假设够用。

5.5 量测维数升高后 inv(S) 数值不稳定

现象：量测向量维度从 1 升到 6 以后，inv(S)偶尔报奇异警告，滤波结果出现间歇性跳变，但明明 S 的理论值不奇异。

原因：S = H P⁻ Hᵀ + R_hat 里，H P⁻ Hᵀ 的秩受状态维度限制，如果量测维数高于状态维数，或者某些量测维度之间高度相关，S 的实际秩会小于名义维度，显式求逆时数值上就表现为奇异。

解决：把K = P_pred * H' * inv(S)改成K = P_pred * H' / S，MATLAB 的右除会走 QR 或 LU 分解而不是显式求逆，数值稳定性明显更好。如果问题依然存在，给 S 的对角线加一个微小的正则化项，比如S = S + 1e-6 * eye(d)，这在状态维数远低于量测维数时几乎是标配操作。另外，检查量测矩阵 H 是否满秩，如果两个量测维度其实是同一个信息的线性组合，先从量测设计层面消掉冗余，比在滤波层面硬扛要干净得多。

6. 进阶验证技巧：用「噪声突变」仿真检验自适应能力的完整方法

验证 VB-AKF 有没有真正在工作，单看一条 R_hat 曲线是不够的，因为单次仿真的随机性可能恰好产生好看的结果。我习惯用三段式指标加蒙特卡洛来下结论，这套方法也适用任何自适应滤波算法的评估。

三段式指标是把仿真时间按突变点切成三段：突变前稳态段、过渡段、突变后稳态段。R 突变发生在第 150 步，那么第 50 到 140 步算突变前，第 150 到 200 步算过渡段，第 200 到 290 步算突变后。分别统计标准 KF 和 VB-AKF 在这三段的位置 RMSE。只看整段平均 RMSE 会严重掩盖问题——突变前的优秀表现把突变后的恶化平均掉了，两个算法整段 RMSE 可能看起来差不多，分段后差距立刻清晰：标准 KF 在突变后 RMSE 可能扩大 5 到 10 倍，VB-AKF 通常只扩大 1.5 到 2 倍。

蒙特卡洛验证是把整个仿真包进一个循环里跑 50 到 100 次，每次换随机种子，最后统计所有实验中 R_hat 在每一时刻的中位数和 25% 到 75% 分位区间。R_hat 的单次轨迹可能抖动，但中位数曲线能反映算法的真实跟踪能力，分位区间宽度则暴露估计稳定性。个别发散样本会严重拉高均值，所以统计量一律取中位数而不是均值。做这个验证时，仿真场景里的突变幅度、突变时刻、Q 和 R0 的取值都应该覆盖一组典型工况，而不是只测一个参数组合，这样得出的结论才是算法层面的结论，而不是某组参数下的巧合。

最后说一个我自己的教训。早期在某个项目里用 VB-AKF，我盯着单次仿真的 R_hat 曲线，看到它跟随得又快又稳，就以为万事大吉。后来换了一个随机种子，结果直接发散，查了半天发现是那次仿真恰好碰上极端噪声序列。从此以后，新算法到手我先跑 100 次蒙特卡洛，看中位数和分位区间，再谈「有效」两个字。自适应滤波器最容易骗人的地方，就是单条曲线看起来很美，但你不确定它是靠算法收敛还是靠运气收敛。VB-AKF 把 R 从「拍脑袋的常量」变成了「在线估计的变量」，这是它值得投入的原因，但它的参数照样需要老老实实调、扎扎实实验证。希望帮到你。

本文还有配套的精品资源，点击获取