简介：面向惯性天文导航与非线性滤波研究的MATLAB算法资源，提供鲁棒UKF/AUKF及其对比算法的完整实现。压缩包共55个文件，以52个.m脚本/函数为主，涵盖无迹变换、预测更新、数值雅可比、马氏距离等核心模块，并附3个txt说明文档；整体仅36KB，轻量易读，便于二次开发。内容既包含EKF、UKF、粒子滤波等多个演示程序，也提供KF的多种更新形式与采样工具，可系统性对比不同滤波在非线性系统中的表现；从无迹变换、数据处理到观测融合均有可直接运行的脚本，尤其适合探索仅用方位信息等受限观测条件下的导航方法。这些代码可用于惯性测量与天文观测数据的融合试验，抑制噪声与模型误差，提升导航解算精度和抗干扰能力。已有349人学习浏览，适合具备一定MATLAB基础、希望深入理解UKF原理并在惯性天文导航场景中实践的研究者与工程师。
1. 惯性天文组合导航里的 UKF 工具：噪声一变，普通滤波器就翻车
惯导短时精度很高，但误差随时间积分增长；星敏感器能提供绝对姿态，却输出频率低、偶尔还带野值。把两者用 UKF 融合，是惯性天文组合导航的常见路线，也是 matlab_ukf_utilities 这类工具包的核心场景。普通 UKF 对量测噪声统计太敏感，星敏感器噪声一变或出现一两个野值，滤波输出就剧烈跳动。AUKF 就是在 UKF 上加噪声自适应、强跟踪和野值拒斥，让滤波更鲁棒。这篇笔记写给做 SINS/CNS 组合预研、算法验证或毕业设计的工程师：在 MATLAB 里搭出能跑的 UKF/AUKF 循环，并知道每个参数动了会有什么后果。
2. 惯性天文组合导航里为什么非用 UKF：先搞清状态和量测长什么样
2.1 惯性递推为什么发散，天文观测为什么能兜底
捷联惯导递推的是姿态、速度和位置。陀螺和加速度计的零偏误差不会消失，而是被积分放大：陀螺漂移 0.01°/h，一小时后姿态误差就能到几十角分量级，位置误差按时间的三次方增长。这不是滤波能解决的，因为惯导系统本身只有相对测量能力，没有绝对参考。
天文导航提供的是绝对参考。星敏感器拍摄星图后匹配星库，直接输出载体相对惯性空间的姿态四元数，测量误差通常在角秒到角分级且不随积分放大。两者特性正好互补：SINS 短期精度高、更新快，CNS 长期稳定、更新慢。组合导航就是让高速的惯导递推在慢速的天文量测辅助下不漂移，滤波器在这里扮演的是数据融合器和误差估计器。
用状态空间语言描述这个场景：状态量取 SINS 误差方程中的姿态误差角 φ、速度误差 δv、位置误差 δp，再加上陀螺漂移 ε 和加速度计零偏 ∇，共 15 维。状态方程就是惯导误差传播方程，它是强非线性的，因为里面有方向余弦矩阵、三角函数和地球自转角速度的耦合。量测量是星敏感器给出的姿态与惯导推算姿态之差，这个差在小角度假设下可以近似线性，但状态方程的非线性没法回避。
2.2 UKF 替代 EKF 的两个理由：不线性化、不抄雅可比
EKF 在这个场景的麻烦不在精度，而在工程实现。状态方程里每个变量对另一个变量的偏导都要手推，姿态误差方程里方向余弦矩阵和三角函数混在一起，写雅可比矩阵时很容易漏项。漏一项在仿真里可能看不出问题，上真实数据就会在某个工况下发散，且发散前毫无征兆。
UKF 的核心是无迹变换（UT）：取 2n+1 个 sigma 点，让这些点穿过非线性方程，再从传播后的点集里重建均值和协方差。整个过程不需要任何解析导数，对非线性的逼近精度可以达到二阶甚至三阶（对高斯分布），而 EKF 只有一阶。对惯导误差传播这种"本身是小量但必须精细传播"的问题，UKF 的样本式传播更符合直觉，也更容易维护。
需要强调的是，UT 不是蒙特卡洛采样，sigma 点的数量和位置是确定性的，只与状态维数 n 和三个参数（alpha、beta、kappa）有关。这也是为什么工具包要把 sigma 点生成单独抽成一个函数：参数要实验、要统一改，写死在主循环里会变成玄学调参。
2.3 sigma 点参数表：alpha、beta、kappa 各自的职责
sigma 点的散布程度由比例 UT 的三个参数决定，它们不是随意设的。
| 参数 | 作用 | 常用值 | 影响 |
|---|---|---|---|
| alpha | 控制 sigma 点离均值的距离 | 1e-2 ~ 1 | 越小点越靠近均值，对强非线性更稳，但过小会损失统计精度 |
| beta | 高阶信息的先验修正 | 2（高斯最优） | 影响协方差权值中高阶项，改它不如改 alpha 来得直观 |
| kappa | 额外的散布尺度 | 3 - n 或 0 | 与 alpha 共同决定 lambda，15 维时常用 0 |
| lambda | 由 alpha、kappa 推导的中间量 | 计算得到 | lambda = alpha^2 * (n+kappa) - n |
对 15 维的 SINS/CNS 组合，我会配 alpha=1e-2、beta=2、kappa=0。alpha 取 1e-2 是因为姿态误差方程在误差较大时非线性明显，sigma 点散布太远会把异常值传播进均值和协方差，造成更新步修正过度。如果状态全是对准良好的小偏差，alpha 取大一点（0.1~1）计算更高效。参数对结果的影响不是线性的，最好离线扫一遍。
2.4 utilities 工具包的函数边界与数据流
一个可复用的 MATLAB 工具包不需要是一个大而全的类库，把函数拆成单一职责就够。下面是我维护的一个最小结构，对应标题里的 utilities 那一层：
| 文件 | 职责 | 关键输入输出 |
|---|---|---|
| ukf_sigma_points.m | 生成 sigma 点与权值 | 输入 x、P，输出 X、Wm、Wc |
| ukf_predict.m | sigma 点经状态方程传播 | 输入 X、Wm、Wc、f、Q、dt，输出 x_pred、P_pred |
| ukf_update.m | 量测更新 | 输入 X、Wm、Wc、z、h、R，输出 x_upd、P_upd、nu、S |
| estimate_R_from_innov.m | 开窗估计量测噪声 | 输入新息缓冲，输出 R_hat |
| fading_factor.m | 强跟踪渐消因子 | 输入 nu、S，输出 lambda |
| outlier_test.m | 野值拒斥 | 输入 nu、S，输出拒斥标记 |
| quat2dcm.m / dcm2quat.m | 姿态转换 | 四元数与 DCM 互转 |
拆开的好处是每个函数都能单独测试。sigma 点函数可以拿一个已知高斯分布验证样本矩；野值拒斥函数可以喂一个刻意做大的新息看它是否触发。我在早期把 sigma 点生成写进主循环，结果每次调参数都要在 200 行代码里找位置，改成独立函数后，调试时间至少省一半。
数据流是这样的：主循环里先调用 ukf_sigma_points 生成点集，传给 ukf_predict 做时间更新；收到星敏量测时，用预测后的分布重新采样点集，传给 ukf_update 做量测更新，同时把新息喂给 adaptive 模块，用来在线调整 R 或触发野值拒斥。自适应模块不直接改状态，只修改主循环下一次迭代用到的参数，这样能保证每个模块可回滚。
3. 在 MATLAB 里搭可复现的 UKF 导航滤波：最小闭环与代码骨架
3.1 15 维 SINS/CNS 状态方程与量测方程怎么写
状态向量定义如下（15 维）：
x = [φ; δv; δp; ε; ∇]
其中 φ 是姿态误差角（3 维），δv 是速度误差（3 维），δp 是位置误差（3 维），ε 是陀螺常值漂移（3 维），∇ 是加速度计零偏（3 维）。
状态方程用 MATLAB 函数给出，输入是当前状态和步长，输出是下个时刻的状态近似。这里我直接用一阶欧拉离散，适合算法验证；上实物数据时建议换成四阶 Runge-Kutta 或解析离散。
function x_next = f_sins_err(x, dt) % 捷联惯导误差传播（15 维，一阶欧拉离散） % x: [姿态误差角(3); 速度误差(3); 位置误差(3); 陀螺漂移(3); 加计零偏(3)] phi = x(1:3); dv = x(4:6); dp = x(7:9); eps = x(10:12); nabla = x(13:15); % 地球自转角速度在导航系的分量（简化，忽略位置变化引起的项） wie_n = [0; 0; 7.292115e-5]; % 姿态误差方程：phi_dot = - cross(phi, wie_n) + eps % 这里用叉乘近似的简化形式，完整版要代入载体运动角速度 phi_dot = - cross(phi, wie_n) + eps; % 速度误差方程：dv_dot = -cross(phi, f_n) + nabla f_n = [0; 0; 9.8]; % 比力在导航系投影，这里取静止近似 dv_dot = - cross(phi, f_n) + nabla; % 位置误差方程 dp_dot = dv; % 陀螺漂移和加计零偏建模为随机游走 eps_dot = zeros(3, 1); nabla_dot = zeros(3, 1); x_next = x + dt * [phi_dot; dv_dot; dp_dot; eps_dot; nabla_dot]; end这段代码最关键的是叉乘符号。姿态误差角方程里，误差传播矩阵是反对称阵，叉乘方向错了，滤波会直接发散射出去。我见过不止一个人把 cross(phi, wie_n) 写成 cross(wie_n, phi)，两者正好差一个负号。另外这里把比力 f_n 近似成重力，是因为算法骨架阶段没有惯导解算输出，验证滤波逻辑够用；真跑组合时，f_n 要由惯导解算实时传入，否则速度误差通道的激励是错的。
量测方程相对简单。星敏感器给出姿态四元数 q_cns，惯导给出 q_ins，两者之间的误差角就是量测：
function z = h_cns(x) % 从状态中的姿态误差角构造量测预测值 % 在误差状态框架下，量测预测就是当前姿态误差角的估计值 z = x(1:3); end这里量测方程等于从状态里直接取前三维，是因为我们把误差状态定义为"惯导输出减去真实值"的差值，而星敏感器提供的正是这个差值的观测。实际工程里还要考虑量测时间戳对齐和星敏感器安装角，算法层面这个量测方程是能直接用的。
3.2 sigma 点生成与权值计算：这段代码别抄错
UT 变换的第一步是生成点集和权值，代码不长，但两处容易出错：协方差阵的 Cholesky 分解要求 P 对称正定；权值公式里 Wc(1) 带 beta 项，Wm(1) 不带，抄串就会导致协方差重建不准确。
function [X, Wm, Wc] = ukf_sigma_points(x, P, alpha, beta, kappa) % 比例无迹变换：生成 2n+1 个 sigma 点及对应权值 % x: n×1 状态均值；P: n×n 状态协方差 % alpha: sigma 点散布因子；beta: 高阶修正；kappa: 额外尺度 n = numel(x); lambda = alpha^2 * (n + kappa) - n; % 对 (n+lambda)*P 做 Cholesky 分解，取下三角 S = chol((n + lambda) * P, 'lower'); X = zeros(n, 2*n + 1); X(:, 1) = x; for i = 1:n X(:, i+1) = x + S(:, i); X(:, i+n+1) = x - S(:, i); end Wm = zeros(1, 2*n + 1); Wc = zeros(1, 2*n + 1); Wm(1) = lambda / (n + lambda); Wc(1) = lambda / (n + lambda) + (1 - alpha^2 + beta); Wm(2:end) = 1 / (2 * (n + lambda)); Wc(2:end) = 1 / (2 * (n + lambda)); end注意 chol 的输入因子是 (n+lambda)P 而不是 P。如果 lambda 是负的（alpha 很小且 n 很大时可能出现），这个系数可能为负，chol 直接报错。我的经验是先把 P 做一次对称化 P = (P+P')/2，再给主对角加 1e-12 的微小正则，能挡掉绝大多数数值问题。kappa 在 n=15 时我取 0，这样 lambda = alpha^2n - n，alpha=0.01 时 lambda 是负数，所以正则项不是可有可无。
权值验证有个简单办法：把 Wm 和 Wc 分别求和，Wm 必须等于 1，Wc 求和理论上接近 1（高斯分布下 beta=2 时恰好等于 1）。每次改参数表后跑这个检测，能掐掉一堆手误。
3.3 UKF 主循环：预测、更新与协方差修正
预测步的核心是让每个 sigma 点穿过状态方程，再加权重建统计量：
function [x_pred, P_pred] = ukf_predict(X, Wm, Wc, f, Q, dt) % 时间更新：sigma 点经状态方程传播 n = size(X, 1); npts = size(X, 2); Xn = zeros(n, npts); for i = 1:npts Xn(:, i) = f(X(:, i), dt); end x_pred = Xn * Wm'; P_pred = zeros(n, n); for i = 1:npts dx = Xn(:, i) - x_pred; P_pred = P_pred + Wc(i) * (dx * dx'); end P_pred = P_pred + Q; % 叠加过程噪声 end这里的过程噪声 Q 直接加在协方差上，没有做任何矩阵分解，所以 Q 必须是正定的。Q 设为零矩阵在某些教科书结论里可行，但工程上会让长时间无量测更新时的 P_pred 越来越小，滤波越来越"自信"，一旦来一个带野值的量测就直接炸。Q 至少要给对角线一个表示随机游走强度的下限。
更新步要算新息协方差 S 和互协方差 Pxz，然后解卡尔曼增益：
function [x_upd, P_upd, nu, S] = ukf_update(X, Wm, Wc, x_pred, P_pred, z, h, R) % 量测更新：把 sigma 点映射到量测空间，做标准卡尔曼更新 n = size(X, 1); nz = numel(z); npts = size(X, 2); Zn = zeros(nz, npts); for i = 1:npts Zn(:, i) = h(X(:, i)); end z_pred = Zn * Wm'; S = zeros(nz, nz); Pxz = zeros(n, nz); for i = 1:npts dz = Zn(:, i) - z_pred; dx = X(:, i) - x_pred; S = S + Wc(i) * (dz * dz'); Pxz = Pxz + Wc(i) * (dx * dz'); end S = S + R; K = Pxz / S; % 卡尔曼增益，相当于 Pxz * inv(S) nu = z - z_pred; % 新息，用于自适应和野值检测 x_upd = x_pred + K * nu; P_upd = P_pred - K * S * K'; endK = Pxz / S 在 MATLAB 里是右除，等价于 Pxz * inv(S)，但数值上比显式求逆稳定。P_upd 用 P_pred - KSK' 形式算完后可能出现小幅不对称，我通常在之后补一行 P_upd = (P_upd + P_upd')/2。
3.4 一个能直接跑通的 SINS/CNS 组合滤波骨架
把上面函数串起来，最小闭环长这样。状态初始化为零向量，协方差对角线根据对准精度和器件指标给粗值，然后进入滤波循环：
% 组合滤波主循环骨架：惯导 100 Hz 递推，星敏 1 Hz 量测更新 clear; rng(1); % 参数初始化 n = 15; alpha = 1e-2; beta = 2; kappa = 0; dt = 0.01; % 惯导/滤波 100 Hz N_cns = 100; % 星敏 1 Hz，每 100 拍更新一次 N_total = 10000; % 仿真 100 秒 % 状态与协方差初值 x = zeros(n, 1); P = blkdiag((0.1*pi/180)^2*eye(3), ... (0.01)^2*eye(3), ... (10)^2*eye(3), ... (0.01*pi/180/3600)^2*eye(3), ... (1e-4)^2*eye(3)); Q = blkdiag((0.005*pi/180)^2*eye(3)*dt, ... (0.001)^2*eye(3)*dt, ... (0.1)^2*eye(3)*dt, ... (1e-6*pi/180/3600)^2*eye(3)*dt, ... (1e-6)^2*eye(3)*dt); R = (3/3600*pi/180)^2 * eye(3); % 星敏 3 arcsec 1σ % 量测真值（假设真实姿态误差角已知，用于生成量测） z_true = zeros(3, N_total / N_cns); for k = 1:N_total % 1) 预测步 [X, Wm, Wc] = ukf_sigma_points(x, P, alpha, beta, kappa); [x_pred, P_pred] = ukf_predict(X, Wm, Wc, @f_sins_err, Q, dt); % 2) 有星敏量测时量测更新 if mod(k, N_cns) == 0 idx = k / N_cns; z = z_true(:, idx) + sqrtm(R) * randn(3, 1); % 模拟量测 % 注意：这里要用预测分布重新生成 sigma 点 [X2, Wm2, Wc2] = ukf_sigma_points(x_pred, P_pred, alpha, beta, kappa); [x_upd, P_upd, nu, S] = ukf_update(X2, Wm2, Wc2, x_pred, P_pred, z, @h_cns, R); x = x_upd; P = P_upd; else x = x_pred; P = P_pred; end % 记录状态用于事后画图（略） end这个骨架里最容易忽略的细节是：量测更新前必须用预测后的 x_pred、P_pred 重新生成 sigma 点，而不是复用预测步的点集。原因在于预测步的 sigma 点已经穿过状态方程，分布已经改变；直接用它们做量测传播，相当于把时间更新和量测更新强耦合，协方差会低估。这个错误在 EKF 里不存在对应物（EKF 直接线性化），所以从 EKF 转过来的人最容易犯。
量测更新频率和预测频率不一致时，滤波曲线会出现周期性的"锯齿"：惯导误差在两次星敏更新之间逐渐累积，量测到达时又被拉回。这是正常现象，不是滤波发散。如果锯齿幅度过大，说明过程噪声 Q 取值偏小或量测噪声 R 偏大，滤波器对量测信任不足。
4. 从 UKF 到 AUKF：把鲁棒性做进去的三个手段
4.1 量测噪声 R 的动态估计：新息序列开窗
普通 UKF 把 R 当常数，但星敏感器的噪声在不同星图、不同亮度、不同动态下并不一致。R 估小了，滤波器过度信任量测，野值一次就能把姿态打偏；R 估大了，天文校正作用被稀释，惯导漂移在长航时里压不住。AUKF 的第一个手段就是用新息序列在线估计 R。
新息定义是 nu = z - z_pred。理论上新息的样本协方差应等于理论新息协方差 S = Pzz + R。因此可以从最近 N 个新息样本的均值估计出 R：
function R_hat = estimate_R_from_innov(nu_buf, S_buf, R_prev, N, R_min, R_max) % 开窗估计量测噪声协方差 % nu_buf : nz×N 最近 N 个新息向量 % S_buf : 每次更新得到的理论新息协方差，存放在 nz×nz×N 数组中 nz = size(nu_buf, 1); % 样本协方差 R_sample = zeros(nz, nz); for i = 1:N R_sample = R_sample + nu_buf(:, i) * nu_buf(:, i)'; end R_sample = R_sample / N; % 减去状态不确定性贡献的部分（用历史 S 减去历史 R 近似） Pzz_mean = zeros(nz, nz); for i = 1:N Pzz_mean = Pzz_mean + (S_buf(:, :, i) - R_prev); end Pzz_mean = Pzz_mean / N; R_hat = R_sample - Pzz_mean; R_hat = (R_hat + R_hat') / 2; % 强制对称 R_hat = R_hat + 1e-9 * eye(nz); % 主对角正则，防止奇异 R_hat = min(max(R_hat, R_min), R_max); % 限幅，防止自适应漂移 end两个参数最影响效果：窗口长度 N 和下限 R_min。N 太小（小于 20），R 的估计抖动大，滤波会在"极其信任量测"和"几乎不信量测"之间摆动；N 太大，噪声统计突变时响应太慢，比如星敏感器从强光遮挡恢复后，R 需要几百拍才回调。我一般取 N=30~50，R_min 取标称 R 的 1/10 到 1/5。
这里有个工程近似要说明：Pzz_mean 用 S_buf - R_prev 来近似，严格推导里应该用 H P H^T 的理论值，但 UKF 没有显式 H 矩阵，直接从每次更新的 S 里剥离 R 是最省事的做法。代价是当滤波器状态收敛后 Pzz 很小，R_hat 近似等于样本新息协方差，这在实际日志回放里表现是够用的。追求严格推导的话就是 Sage-Husa 族的递推式，本质差别不大。
4.2 强跟踪渐消因子：新息异常变大时的"悔过"机制
R 自适应解决的是噪声统计慢变问题，但组合导航里还有一类更棘手的场景：星敏感器数据本身没问题，是载体机动或模型简化导致状态方程失配。这时候新息也会异常大，但问题不在 R 而在预测协方差被低估。直接调 R 是治标不治本。
强跟踪滤波（Fading UKF）的思路是：检测到新息能量异常时，在预测协方差上乘一个渐消因子 lambda，人为放大模型不确定性，让滤波器对新的量测更开放：
function lambda = fading_factor(nu, S, gamma) % 强跟踪渐消因子计算（简化版） % gamma: 新息白噪声阈值，通常取 chi2inv(0.99, nz) nz = numel(nu); md = nu' / S * nu; % 马氏距离 if md > gamma lambda = md / nz; % 偏离越大，预测协方差放大倍数越大 lambda = min(lambda, 5.0); % 限幅，防止一步就崩 else lambda = 1.0; end end使用时机：先做野值拒斥，再做渐消因子判定。如果新息是因为野值产生的，拒斥后根本不会走到更新步；能走到这里的异常新息，大概率来自状态预测不准，这时放大预测协方差是合理的后悔药。lambda 限幅我放到 5 以内，超过这个值意味着系统模型已经不能解释当前量测，再放大也没有意义，应该去查状态方程或时间同步。
严格强跟踪里渐消因子由两级滤波递推计算，表达式里包含 V(k) 平滑估计。我给的 md/nz 是工程简化，因为它复用已经算好的马氏距离，不需要额外状态。代价是这个因子对突变响应更激进，所以限幅必须跟上。两种我都试过，简化版在 SINS/CNS 场景下结果差异很小，但代码量少一个量级。
4.3 野值拒斥：星敏感器被遮挡时别让它带偏滤波器
星敏感器最典型的野值来自两个地方：一是飞行器机动导致星图拖尾，二是太阳或地气光辉进入视场，星点提取失败或匹配出错。野值的特点是新息马氏距离远大于正常分布。卡方检验是最直接的拒斥手段：
function flag = outlier_test(nu, S, chi2_threshold) % 野值检测：新息马氏距离服从自由度为 nz 的卡方分布 % nu: 新息向量；S: 新息协方差；chi2_threshold: 卡方分位数 md = nu' / S * nu; flag = md > chi2_threshold; end自由度是量测维数，姿态量测为 3 时，chi2inv(0.997, 3) 约 14.16。阈值设多少要结合星敏感器的真实故障概率：故障率高于千分之三时固定阈值会产生大量误拒斥，我习惯把它做成参数，用一段不含野值的真实日志先标定正常新息的分布，再反推阈值。
拒斥后的策略有三种：直接跳过更新保持预测；保持预测但把 R 暂调大 10 倍再更新；连续 N 拍拒斥后切纯惯导模式并报警。我一般推荐第一种加第三种：单拍野值直接丢弃，连续 5 拍以上拒斥说明不是野值而是系统性故障，此时还硬融合只会在姿态里引入持续偏置。切换纯惯导后，要靠地面注入或重新捕获星图来恢复，这个逻辑要在组合导航状态机里单独写。
4.4 AUKF 参数表与初始化建议
把三块自适应机制合成 AUKF 后，参数比普通 UKF 多了不少，下面是我经过多轮回放后固定下来的默认值，可以直接当起点：
| 参数 | 推荐值 | 含义与调整方向 |
|---|---|---|
| alpha / beta / kappa | 1e-2 / 2 / 0 | sigma 点分布，不要轻易动 |
| R 开窗长度 N | 30 ~ 50 | 越长越稳，越短响应越快 |
| R 下限 | 标称 R 的 1/10 | 防止自适应把 R 压到接近零 |
| R 上限 | 标称 R 的 10 倍 | 防止把模型失配误判为噪声变大 |
| 渐消因子限幅 | 1 ~ 5 | 超过 5 直接报警查模型 |
| 野值卡方阈值 | chi2inv(0.997, nz) | 故障率越高阈值应越小 |
| 连续拒斥切换门限 | 5 拍 | 连续 5 拍拒斥转纯惯导 |
初始化顺序有讲究：先把 R 自适应关掉，用固定 R 跑通基线；再单独开野值拒斥，确认拒斥率符合预期；最后开 R 自适应和渐消因子。一次全部打开，出问题时分不清是谁引起的。每开一个模块就要留一段对应的记录，比如把新息马氏距离、R_hat 估计值、lambda 变化都存下来，回放时一眼能定位到是哪块在起作用。
5. 惯性天文 UKF 调试避坑：五个让我返工过的实际问题
调试 UKF 的判断依据只有一个核心：新息序列是否和白噪声一致。以下问题在现象层最终都表现为新息统计异常，但触发原因各不相同。前四条在 MATLAB 2023b 和 2026b 两个版本上都遇到过，最后一条是换平台移植时踩出来的。每条按现象、原因、解决的顺序写，方便你拿着现象对照自己的曲线。
5.1 协方差阵非正定，chol 报错后滤波直接中断
现象：跑一段时间后 ukf_sigma_points 里 chol((n+lambda)*P) 报错，提示矩阵不是正定；把断点设在出错前看 P，对角线全是正常量级，但矩阵已经不对称。
原因：两条路径叠加。一是数值误差，P 在反复做 P - KSK' 后失去对称性，负特征值被累积放大；二是长时间没有量测更新时，Q 对 P 的驱动不足以抵消初始协方差的收缩，P 接近奇异。
解决：三层防护。第一层，每次更新后强制 P = (P+P')/2，把对称性修好；第二层，chol 前给 P 加 1e-12 * eye(n)，正则量级远小于状态方差，不影响滤波精度；第三层，把 Q 的下限设成非零，保证任何时刻 P 都有正定增量。这三层我都加在工具函数内部，主循环里就不会因为这些数值问题返工。我在 MATLAB 2023b 和 2026b 上都复现过同样的现象，结论一致：这是数值问题，不是滤波逻辑问题。
5.2 星敏更新频率太低，滤波曲线出现周期性大幅锯齿
现象：仿真曲线在两次星敏更新之间缓慢漂移，量测更新瞬间拉回，拉回幅度在后半段越来越大，看起来像收敛不彻底。
原因：SINS 递推 100 Hz，星敏 1 Hz，两拍之间滤波只有预测没有量测。惯性误差在 1 秒内累积的量如果和 R 不匹配，就会表现出"快到更新点时误差已经很大"的锯齿。锯齿大不是发散，是预测协方差与真实误差不匹配的信号。
解决：不要靠调 Q 硬压锯齿，先检查时间同步。星敏输出是曝光时刻还是读出时刻，MATLAB 里模拟时如果不对齐，等于给量测人为加了延迟，滤波永远收敛不到静止点。其次把锯齿幅度和理论预测协方差对角线对比：如果实际误差超过 3 倍 sigma，说明过程噪声 Q 低估了陀螺噪声；如果远小于 sigma，说明 Q 高估，滤波偏保守。按这个口径调 Q，比肉眼调曲线靠谱得多。
5.3 野值没拒斥，姿态跳变一次后几条曲线都不回来
现象：无自适应版本在第 400 拍附近出现一次很大的尖峰后，姿态误差曲线没有回到零附近，而是收敛到一个偏置，直到下一次大机动才被纠正。
原因：星敏感器受阳光影响输出一个明显错误姿态，新息是正常值的 10 倍以上。普通 UKF 没有拒斥机制，直接把野值当成可信量测，把状态和协方差一起拉偏。协方差还因此收缩，后面正确量测也拉不回来，滤波器失去了对自身不确定度的正确认知。
解决：在更新步之前加 outlier_test，拒斥后不执行更新，保持预测结果。这里要说明一个容易混淆的点：野值拒斥和渐消因子都看马氏距离，但作用相反。野值拒斥是丢弃量测不更新，渐消因子是保留量测但放大预测协方差。先判野值，后判渐消，顺序不能反，反了会把野值当成模型失配，反而加剧污染。
5.4 自适应 R 没限幅，滤波器把 R 估计到接近零
现象：开窗估计的 R_hat 在一段新息连续偏小的阶段持续下降，最后趋近于零。此时滤波器几乎完全信任星敏量测，下一拍一个普通扰动就能让姿态估计剧烈震荡。
原因：新息序列在滤波收敛后存在时间相关性（实际系统里通常是有色噪声），不满足自适应估计的"新息白噪声"假设。开窗样本协方差因此被低估，R_hat 也被带偏。
解决：给 R_hat 加硬限幅是最省事的兜底，R_min 设成标称 R 的十分之一，R_max 设成标称 R 的十倍。更接近本质的做法是先把新息序列做白化检验（比如看自相关函数是否落在置信区间内），如果有色性明显，就只能降低自适应权重或改用遗忘因子更小的递推形式。工程上我两个都做：限幅兜底，白化检验做在线监控，两者独立，互不干扰。
5.5 初始协方差拍脑袋设置，前几十秒剧烈震荡甚至发散
现象：仿真前 30 秒姿态误差曲线振幅逐渐增大，最大到好几度，之后才被量测拉回；有些随机种子下直接发散，后面的结果全部不可用。
原因：初始 P 给得过大，比如姿态误差初值实际有 10 度却给了 0.1 度量级的初始方差，滤波器对自己的初值过度自信，早期量测无法快速修正；反之初始 P 过小，收敛速度差，前期误差主导曲线。二者都会让滤波初期进入非线性强的区域，UT 的局部近似失效。
解决：初始 P 应该按实际对准工艺给，常见做法是静态对准后姿态误差大概角分级，就给角分级的方差；速度误差按载体运动状态给 0.01~0.1 m/s。不要因为"不知道就给大一点"而盲目放大，协方差大不是后悔药，是让滤波器把真实值当异常值的开门证件。我的一般习惯是先用一段实测静止数据做离线标定：把初始 P 扫 3 个数量级，看哪组能在最短时间内使新息进入白噪声区间，就用那组上系统。
6. 把 AUKF 的结论跑成可判读的数字：蒙特卡洛与日志回放
验证 AUKF 是否真的比 UKF 鲁棒，不能只跑一条仿真曲线。固定噪声的一次仿真对参数有偶然性，给同一份代码换随机种子，表现可能差好几倍。我习惯跑 50 次蒙特卡洛，每次重新生成量测噪声和过程噪声序列，统计姿态误差的 RMS 与最大偏差：
% 蒙特卡洛验证骨架：50 次仿真，统计姿态误差 RMS 和最大偏差 n_mc = 50; rms_phi = zeros(n_mc, 3); max_phi = zeros(n_mc, 3); for mc = 1:n_mc rng(mc * 100); % 固定种子，保证可复现 % 运行第 3.4 节的主循环，每次重新生成噪声序列 % 记录姿态误差时间序列 phi_err rms_phi(mc, :) = sqrt(mean(phi_err.^2, 1)); max_phi(mc, :) = max(abs(phi_err), [], 1); end % 判读标准： % 姿态 RMS 低于 10 arcsec，最大偏差低于 30 arcsec，说明滤波统计上可信判读时把 AUKF 和普通 UKF 各跑一遍，重点看两类指标：平均 RMS 是否改善，最大偏差是否被压住。鲁棒性体现在后者：正常 UKF 在 50 次里可能有一两次发散或出轨，AUKF 把最大偏差从度级压到角分级，这才是它存在的意义。平均 RMS 反而可能比手动调好的 UKF 略差，因为自适应在噪声平稳时牺牲了一点点最优性换取稳定性，这个取舍要提前向使用者讲清楚。
回放实际日志时，我还会盯一个新指标：新息的白噪声程度。把新息序列做自相关，正常情况下只有零延迟显著；如果出现周期性相关，通常对应时间同步误差或载体周期运动建模缺失。整条链跑通后，我的习惯是把每次调参后的 R_hat、lambda、拒斥计数一起落盘，下次排故障不用重跑数据。这个习惯帮我省过不少返工，希望也能帮到你。
本文还有配套的精品资源，点击获取
