你知道吗,Kalman Filter 的公式至少有四种解释方式。相应的,也就有至少四种不同理解它的办法。 明明最终推导出来都是同一套公式,有必要像孔乙己“回“字的四种写法那样掉书袋吗? 不,提供多种角度是为了理解 KF 的本质,给工程化参数调优提供一份理论支撑,避免出现“调无可调“的局面。知其然,知其所以然。

从四个角度推导 Kalman Filter:Time Update 与 Measurement Update

Kalman Filter 的公式并不长,七八个式子,抄一遍就能跑起来。麻烦的是它的推导有太多个版本,相信你也有过类似的感受:概率教材从贝叶斯公式讲起;导航和控制的工程手册直接定义估计误差,估计是否无偏、是否最优;数理统计的书两句话带过,说这不过是联合高斯分布的条件化;而学过 SLAM 技术的话,因子图优化、光束平差、平方根信息滤波等概念又把 KF 解释为带动态先验的递推加权最小二乘。

这些说法都对,而且指向同一组公式。但如果只学过其中一种,遇到另一种时往往认不出来,更糟的是会在需要改动滤波器时不知道该动哪里——为什么改了增益就必须换 Joseph 形式,为什么量测明明够多协方差却降不下去,为什么 UKF 里没有 $H$ 矩阵却还能算增益,这些问题的答案分散在不同的视角里。

我会把四条路线放在同一套符号和假设下各走一遍:

  1. 贝叶斯递推:把滤波看成条件分布的递推,解释为什么 KF 天然分成预测和校正两步。
  2. 估计误差传播与最优增益:只用一阶和二阶矩,解释每个协方差公式怎么来的、增益在什么意义上最优。
  3. 仿射变换与联合高斯条件化:只调用高斯分布的两条封闭性结论,是最短的一条路径。
  4. MAP 与加权最小二乘:把每个历元看成一次二次优化,连接信息滤波和因子图。

第 5 章把四条路线的中间量逐一对应起来,说明它们为什么必然汇合;第 6 章讨论假设被打破时公式该怎么改,以及各条路线分别通向哪一类扩展算法。

0. 问题定义与统一符号

四条推导路线能互相印证或者解释的前提,是它们从同一个问题出发,为便于对比,我们定义同一套记号。后面每一章直接引用,不再重复定义。读到任何一条路线出现陌生符号时,都可以翻回这里对照。

0.1 线性状态空间模型

Kalman Filter 处理的是离散时间线性高斯系统,它由两个方程描述。状态方程说明状态如何从上一历元演化到当前历元:

\[x_k = F_{k-1}x_{k-1}+B_{k-1}u_{k-1}+w_{k-1}\]

量测方程说明当前状态如何被观测到:

\[z_k = H_k x_k + v_k\]

各个符号的含义和维度如下:

  • $x_k\in\mathbb{R}^n$:历元 $k$ 的状态向量。GNSS 里通常是位置、速度、接收机钟差与钟漂,有时再加对流层延迟、频间偏差或载波相位模糊度。
  • $z_k\in\mathbb{R}^m$:历元 $k$ 的量测向量,例如各卫星的伪距、多普勒或载波相位(在线性化滤波中一般是它们相对预测值的残差)。
  • $F_{k-1}\in\mathbb{R}^{n\times n}$:状态转移矩阵,把 $k-1$ 历元的状态推到 $k$ 历元。
  • $B_{k-1}u_{k-1}$:已知的确定性输入项。纯 GNSS 滤波通常没有这一项,组合导航里它可以是 IMU 比力和角速度积分带来的已知增量。
  • $H_k\in\mathbb{R}^{m\times n}$:量测矩阵,GNSS 里每一行往往是卫星视线方向的单位矢量加上钟差系数。
  • $w_{k-1}\sim\mathcal{N}(0,Q_{k-1})$:过程噪声,描述状态演化模型本身的不确定性。
  • $v_k\sim\mathcal{N}(0,R_k)$:量测噪声,描述观测本身的不确定性。

有三个约定需要特别强调。

第一,噪声下标跟随区间而不是历元。$w_{k-1}$ 表示作用在 $k-1\to k$ 这一段时间上的过程噪声,所以它和 $F_{k-1}$、$Q_{k-1}$ 用同一个下标;而 $v_k$ 是历元 $k$ 当场产生的量测噪声。后面推导误差传播时,交叉项是否为零和这个时间对应关系直接相关。

第二,$F$、$B$、$Q$ 全部带下标,是因为采样间隔 $\Delta t$ 未必恒定。GNSS 场景里丢历元、变频输出都很常见,把矩阵写成时变形式可以避免推导时把 $\Delta t$ 当常数处理。$H_k$、$R_k$ 同理:每个历元可见卫星数不同,$m$ 本身就是变化的。

第三,$B_{k-1}u_{k-1}$ 是完全已知的确定性项,不含随机成分。它只会平移状态的均值,不会改变任何协方差。四条推导路线里它都只出现在均值公式中,所以后面讨论协方差时会直接把它略去,这不是遗漏。

0.2 需要明确的基本假设

标准 Kalman Filter 的公式之所以是闭式的,靠的是下面这几条假设。它们在推导里会被反复用到,每次用到时都会点明是哪一条。

  • 初始后验高斯:上一历元的后验分布是高斯分布,即 $p(x_{k-1}\mid z_{1:k-1})=\mathcal{N}(\hat x_{k-1}^+,P_{k-1}^+)$。线性高斯系统下这一条会自我保持:只要初始时刻成立,之后每个历元都成立。
  • 噪声零均值且白:$E[w_k]=0$,$E[v_k]=0$,并且 $E[w_iw_j^T]=Q_i\delta_{ij}$,$E[v_iv_j^T]=R_i\delta_{ij}$。也就是不同历元之间的噪声互不相关,噪声序列没有时间相关性。
  • 过程噪声与历史独立:$w_{k-1}$ 与 $x_{k-1}$、与上一历元的估计误差 $\tilde x_{k-1}^+$、与全部历史量测 $z_{1:k-1}$ 都不相关。Time Update 里交叉项消失就靠这一条。
  • 量测噪声与状态独立:$v_k$ 与 $x_k$、与预测误差 $\tilde x_k^-$ 不相关,且 $E[w_iv_j^T]=0$。Measurement Update 里 $R_k$ 能干净地加到新息协方差上,就靠这一条。
  • Markov 性与条件独立:给定 $x_{k-1}$ 后 $x_k$ 与更早的历史无关,给定 $x_k$ 后 $z_k$ 与其他历元的量测无关。这一条是贝叶斯递推能写成两步的结构性前提。
  • 协方差正定性:$P_{k-1}^+\succ 0$、$R_k\succ 0$、$Q_{k-1}\succeq 0$。$R_k$ 正定保证新息协方差 $S_k$ 可逆,增益有定义。

这里还有一个容易被忽略的层次差别:高斯假设并不是所有路线都必需。第 1 章的贝叶斯递推和第 4 章的 MAP 需要完整的分布假设,因为它们直接操作概率密度;第 2 章的误差传播只用到一阶和二阶矩,即便噪声不是高斯的,得到的结果仍然是最优线性无偏估计,只是不再是全局最优估计。第 3 章介于两者之间:仿射变换部分只需要矩,条件化部分才需要联合高斯。这个差别在第 5 章汇总时还会再提。

本文主体按上面这组标准独立噪声假设推导。工程中经常被违反的情形——过程噪声与量测噪声互相关、多接收机钟差共享晶振导致 $Q$ 出现非对角块、量测之间存在相关性——留到第 6 章统一讨论,说明公式需要补哪些交叉协方差项。

0.3 估计量与误差的记号

Kalman Filter 的每个历元都有“量测前”和“量测后”两个状态,符号上必须区分清楚,否则推导中很容易串。本文统一采用上标 $-$ 表示先验(Time Update 之后、Measurement Update 之前),上标 $+$ 表示后验(吸收了当前历元量测之后):

\[\hat x_k^-=E[x_k\mid z_{1:k-1}],\qquad \hat x_k^+=E[x_k\mid z_{1:k}]\]

其中 $z_{1:k}={z_1,z_2,\dots,z_k}$ 表示从第一个历元到第 $k$ 个历元的全部量测。相应的估计误差定义为真值减估计值:

\[\tilde x_k^-=x_k-\hat x_k^-,\qquad \tilde x_k^+=x_k-\hat x_k^+\]

协方差矩阵则是估计误差的二阶矩:

\[P_k^-=E\left[\tilde x_k^-(\tilde x_k^-)^T\right],\qquad P_k^+=E\left[\tilde x_k^+(\tilde x_k^+)^T\right]\]

注意 $P$ 描述的是估计误差的不确定性,不是状态本身的不确定性。这两者在无偏估计下数值相同,但含义不同,讨论滤波器一致性时区别很重要。此外,$\hat{\cdot}$ 表示估计值,$\tilde{\cdot}$ 表示误差,$I$ 表示与上下文维度匹配的单位阵,$|a|_W^2=a^TWa$ 表示以 $W$ 为权矩阵的加权平方范数(第 4 章会大量使用)。

0.4 两步递推及目标公式

有了上面的记号,Kalman Filter 就是下面这样一个两步递推:每个历元先用动力学模型把上一历元的后验推到当前历元,得到先验;再用当前历元的量测把先验校正成后验。

Time Update 传播均值和协方差:

\[\hat x_k^- = F_{k-1}\hat x_{k-1}^+ + B_{k-1}u_{k-1}\] \[P_k^- = F_{k-1}P_{k-1}^+F_{k-1}^T+Q_{k-1}\]

Measurement Update 先算新息及其协方差,再算增益,最后更新状态和协方差:

\[\nu_k=z_k-H_k\hat x_k^-\] \[S_k=H_kP_k^-H_k^T+R_k\] \[K_k=P_k^-H_k^TS_k^{-1}\] \[\hat x_k^+=\hat x_k^-+K_k\nu_k\] \[P_k^+=(I-K_kH_k)P_k^-\]

其中 $\nu_k$ 称为新息(innovation),是实际量测与量测预测之差;$S_k$ 是新息协方差;$K_k$ 是 Kalman 增益。最后一式通常称为协方差更新的简式,它只在 $K_k$ 取最优值时成立。数值上更稳健的是 Joseph 形式:

\[P_k^+=(I-K_kH_k)P_k^-(I-K_kH_k)^T+K_kR_kK_k^T\]

Joseph 形式对任意增益 $K_k$ 都成立,而且右侧是两个对称半正定项之和,浮点运算下更容易维持 $P_k^+$ 的对称性和正半定性。第 2 章会说明这两种形式如何互相转化,第 6 章再补充实现层面的取舍。

这七个式子就是全文的终点。接下来四章分别从贝叶斯递推、估计误差传播、仿射变换与联合高斯条件化、MAP 与加权最小二乘出发,每一章都独立走到这同一组公式;第 5 章再把四条路线上的中间量一一对应起来。


1. 角度一:贝叶斯递推

这一章回答:KF 为什么天然分成“先预测、再校正”两步?

贝叶斯视角不把状态看成一个待求的数,而看成一个随机变量;滤波器维护的也不是一个点估计,而是这个随机变量在“已知全部历史量测”条件下的概率分布 $p(x_k\mid z_{1:k})$。一旦这样定义问题,两步结构就不是人为设计出来的实现技巧,而是两件本质不同的事情在概率上的必然分工:时间往前走一步,和一批新量测到达。前者引入一个新的随机变量并丢弃旧的,对应边缘化,也就是积分;后者对同一个随机变量施加新的约束,对应似然相乘。KF 的 Time Update 和 Measurement Update 分别就是这两个操作在线性高斯假设下的闭式结果。

1.1 递推的起点:上一历元的后验分布

递推必须有起点。按假设 0.2 的第一条,上一历元结束时我们持有的是一个高斯后验:

\[p(x_{k-1}\mid z_{1:k-1}) = \mathcal{N}\!\left(x_{k-1};\,\hat x_{k-1}^+,\,P_{k-1}^+\right)\]

这里 $\mathcal{N}(x;\mu,\Sigma)$ 表示以 $x$ 为自变量、均值 $\mu$、协方差 $\Sigma$ 的高斯密度。滤波开始时($k=0$)这个分布由初始化给出,即 $p(x_0)=\mathcal{N}(x_0;\hat x_0,P_0)$,GNSS 里通常用一次单点定位的最小二乘解和它的协方差来填。

这一步看着像废话,其实是整套推导的关键:接下来两小节要证明的,本质上是一个不变量。只要输入是高斯分布,Time Update 和 Measurement Update 的输出仍然是高斯分布。有了这个封闭性,滤波器就不需要存储一个函数,只需要存储 $(\hat x,P)$ 这两个参数——它们是后验分布的充分统计量。一般的贝叶斯滤波做不到这一点:任意分布经过一次非线性传播和一次似然加权后,形状会越来越复杂,理论上需要无穷多个参数才能精确描述。KF 能做成 $O(n^2)$ 内存的递推算法,全靠高斯族对下面这些操作封闭。

1.2 Time Update:Chapman–Kolmogorov 方程

Time Update 要求的是 $p(x_k\mid z_{1:k-1})$:条件里只有到 $k-1$ 为止的量测,$z_k$ 还没用上。$x_k$ 并不直接出现在已有的后验里,它是通过 $x_{k-1}$ 间接相关的,所以先把 $x_{k-1}$ 引进来再积掉。

用全概率公式:

\[p(x_k\mid z_{1:k-1}) = \int p(x_k, x_{k-1}\mid z_{1:k-1}) \,dx_{k-1} = \int p(x_k\mid x_{k-1},z_{1:k-1})\, p(x_{k-1}\mid z_{1:k-1}) \,dx_{k-1}\]

这里 $x_{k−1}$ 的引入可以从边缘化的角度去理解:$x_k$ 的条件先验分布,是 ($x_k, x_{k−1}$) 条件联合分布关于 $x_{k−1}$ 的边缘分布。即按照上一时刻状态的所有可能情况,对当前状态事件进行分解。

再用 Markov 性:给定 $x_{k-1}$ 之后,$x_k$ 与更早的历史无关,因此 $p(x_k\mid x_{k-1},z_{1:k-1})=p(x_k\mid x_{k-1})$,于是

\[p(x_k\mid z_{1:k-1}) = \int p(x_k\mid x_{k-1})\, p(x_{k-1}\mid z_{1:k-1}) \,dx_{k-1}\]

这就是 Chapman–Kolmogorov 方程。它对任何 Markov 过程都成立,和高斯与否无关——高斯假设只决定这个积分能不能算出闭式解。

现在代入具体模型。状态方程 $x_k=F_{k-1}x_{k-1}+B_{k-1}u_{k-1}+w_{k-1}$ 给出转移密度:

\[p(x_k\mid x_{k-1}) = \mathcal{N}\!\left(x_k;\,F_{k-1}x_{k-1}+B_{k-1}u_{k-1},\,Q_{k-1}\right)\]

于是积分号里是两个高斯密度的乘积,形式上是一个卷积。硬算的话,要把指数上两个二次型合并、对 $x_{k-1}$ 配方,把与 $x_{k-1}$ 有关的部分凑成一个完整的高斯核(对 $x_{k-1}$ 积分后是常数),剩下的部分就是关于 $x_k$ 的二次型。标量情形下这个结果很熟悉:两个高斯卷积后仍是高斯,均值相加、方差相加。

不过有一条更省事的路。被积函数是 $(x_{k-1},x_k)$ 的联合高斯密度,而积分掉一个分量就是求边缘分布;联合高斯的边缘分布必定是高斯,所以我们只需要算出它的均值和协方差,不必真的做积分。均值直接对状态方程取条件期望,注意 $w_{k-1}$ 零均值且与 $z_{1:k-1}$ 无关:

\[\hat x_k^- = E\!\left[x_k\mid z_{1:k-1}\right] = F_{k-1}E\!\left[x_{k-1}\mid z_{1:k-1}\right]+B_{k-1}u_{k-1} = F_{k-1}\hat x_{k-1}^++B_{k-1}u_{k-1}\]

协方差则把状态方程减去上式,确定性项 $B_{k-1}u_{k-1}$ 正好抵消:

\[x_k-\hat x_k^-=F_{k-1}\left(x_{k-1}-\hat x_{k-1}^+\right)+w_{k-1}\]

预测协方差定义为条件二阶矩 $P_k^-=E!\left[(x_k-\hat x_k^-)(x_k-\hat x_k^-)^T\mid z_{1:k-1}\right]$。把上式代入并展开:

\[\begin{aligned} P_k^- &= E\!\left[ \bigl(F_{k-1}(x_{k-1}-\hat x_{k-1}^+)+w_{k-1}\bigr) \bigl(F_{k-1}(x_{k-1}-\hat x_{k-1}^+)+w_{k-1}\bigr)^T \mid z_{1:k-1} \right] \\ &= F_{k-1}\, E\!\left[ (x_{k-1}-\hat x_{k-1}^+)(x_{k-1}-\hat x_{k-1}^+)^T \mid z_{1:k-1} \right] F_{k-1}^T \\ &\quad + F_{k-1}\, E\!\left[ (x_{k-1}-\hat x_{k-1}^+)w_{k-1}^T \mid z_{1:k-1} \right] \\ &\quad + E\!\left[ w_{k-1}(x_{k-1}-\hat x_{k-1}^+)^T \mid z_{1:k-1} \right] F_{k-1}^T \\ &\quad + E\!\left[ w_{k-1}w_{k-1}^T \mid z_{1:k-1} \right] \end{aligned}\]

第一项正是 $F_{k-1}P_{k-1}^+F_{k-1}^T$,第四项由过程噪声定义就是 $Q_{k-1}$。中间两项是交叉项。由假设,$w_{k-1}$ 与历史独立——它既独立于 $z_{1:k-1}$,也独立于由历史决定的状态 $x_{k-1}$(因而也独立于误差 $x_{k-1}-\hat x_{k-1}^+$)。再加零均值 $E[w_{k-1}]=0$,交叉项期望为零,于是只剩对角两项:

\[P_k^- = F_{k-1}P_{k-1}^+F_{k-1}^T+Q_{k-1}\]

合起来就是

\[p(x_k\mid z_{1:k-1}) = \mathcal{N}\!\left(x_k;\,\hat x_k^-,\,P_k^-\right)\]

概率上的图像很清楚:Time Update 做的是把概率质量沿着动力学搬运,同时抹开。$F_{k-1}$ 负责搬运和形变,$Q_{k-1}$ 负责抹开。因为 $Q_{k-1}\succeq 0$,这一步只会让分布变得更宽——没有新信息进来,不确定性不可能减少。这也是为什么长时间失锁后 GNSS 滤波的协方差会一路涨上去。

1.3 Measurement Update:贝叶斯公式

新的量测 $z_k$ 到达后,要求的是 $p(x_k\mid z_{1:k})$。把 $z_{1:k}$ 拆成 $z_k$ 和 $z_{1:k-1}$,对 $z_k$ 用贝叶斯公式:

\[p(x_k\mid z_{1:k}) = \frac{p(z_k\mid x_k,z_{1:k-1})\,p(x_k\mid z_{1:k-1})}{p(z_k\mid z_{1:k-1})}\]

再用条件独立性:给定 $x_k$ 之后,$z_k$ 只由量测方程决定,与历史量测无关,即 $p(z_k\mid x_k,z_{1:k-1})=p(z_k\mid x_k)$。分母不含 $x_k$,只是归一化常数,所以

\[p(x_k\mid z_{1:k}) \propto p(z_k\mid x_k)\, p(x_k\mid z_{1:k-1})\]

这就是“先验乘似然”。注意这里的先验不是滤波初值,而是上一小节刚算出来的预测分布:相对于本历元的量测 $z_k$ 而言,它就是先验。似然由量测方程给出:

\[p(z_k\mid x_k)=\mathcal{N}\!\left(z_k;\,H_kx_k,\,R_k\right)\]

两个高斯密度相乘,指数相加。丢掉与 $x_k$ 无关的常数,取负对数得到

\[J(x_k) = \frac{1}{2}\left(x_k-\hat x_k^-\right)^T\left(P_k^-\right)^{-1}\left(x_k-\hat x_k^-\right) + \frac{1}{2}\left(z_k-H_kx_k\right)^TR_k^{-1}\left(z_k-H_kx_k\right)\]

关键观察是:$J$ 是 $x_k$ 的二次函数,所以 $p(x_k\mid z_{1:k})\propto e^{-J(x_k)}$ 仍然是高斯密度。既然是高斯,它的协方差就是 $J$ 的 Hessian 的逆,均值就是 $J$ 的极小点。直接求二阶导:

\[\left(P_k^+\right)^{-1} = \nabla^2J = \left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k\]

令一阶导为零:

\[\left[\left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k\right]\hat x_k^+ = \left(P_k^-\right)^{-1}\hat x_k^-+H_k^TR_k^{-1}z_k\]

这一对式子叫信息形式,读起来非常直观:后验精度等于先验精度加上量测带来的精度,信息是可加的;后验均值则是先验和量测按各自精度加权后的折中。第 4 章会从 MAP 准则出发再次得到同一组式子,因为负对数后验的极小点正是 MAP 估计——在高斯情形下它和后验均值重合。

剩下的只是代数变形。对协方差用 Woodbury 矩阵恒等式(详细展开见 4.4 节):

\[P_k^+ = \left[\left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k\right]^{-1} = P_k^--P_k^-H_k^T\left(H_kP_k^-H_k^T+R_k\right)^{-1}H_kP_k^-\]

括号里正是新息协方差 $S_k=H_kP_k^-H_k^T+R_k$,把 $K_k=P_k^-H_k^TS_k^{-1}$ 代进去就得到 $P_k^+=(I-K_kH_k)P_k^-$。均值同理:把 $\hat x_k^+=P_k^+\left[(P_k^-)^{-1}\hat x_k^-+H_k^TR_k^{-1}z_k\right]$ 展开整理,可以化成

\[\hat x_k^+=\hat x_k^-+K_k\left(z_k-H_k\hat x_k^-\right)\]

至此第 0 章列出的全部公式都从概率密度递推中出来了。

还有一个副产品值得留意:刚才被当成归一化常数丢掉的分母,本身是有意义的。对预测分布积分可得

\[p(z_k\mid z_{1:k-1}) = \int p(z_k\mid x_k)p(x_k\mid z_{1:k-1})\,dx_k = \mathcal{N}\!\left(z_k;\,H_k\hat x_k^-,\,S_k\right)\]

也就是说,新息 $\nu_k=z_k-H_k\hat x_k^-$ 服从零均值、协方差为 $S_k$ 的高斯分布。这一条是很多工程手段的理论依据:新息的归一化平方 $\nu_k^TS_k^{-1}\nu_k$ 服从自由度为 $m$ 的卡方分布,可用于 NIS 一致性检验和抗差剔除;把各历元的 $p(z_k\mid z_{1:k-1})$ 连乘则得到整段数据的边缘似然,是噪声参数辨识、IMM 模型概率更新的计算基础。

1.4 本角度重点解释的问题

为什么 Time Update 是积分,而 Measurement Update 是相乘? 因为两者处理的事情不同。时间前进意味着状态从 $x_{k-1}$ 换成了 $x_k$,我们必须引入一个新变量并把旧变量丢掉,丢掉的操作就是边缘化,也就是积分。量测到达则没有引入新的状态变量,只是对已有的 $x_k$ 增加了一个约束,约束的强度由似然函数描述,叠加约束的操作就是相乘。一句话:积分对应“换变量”,相乘对应“加信息”。

为什么预测分布叫先验,更新结果叫后验? 先验和后验是相对于本历元的量测说的,不是相对整段数据说的。$p(x_k\mid z_{1:k-1})$ 已经吸收了全部历史量测,从整段数据看它绝不是无信息的;但它还没看到 $z_k$,所以在处理 $z_k$ 这一步里它扮演先验的角色。滤波器的每个历元都是一次完整的贝叶斯推断,上一次的后验经动力学传播后成为下一次的先验,这正是“递推”二字的含义。

为什么线性、高斯、Markov 三个假设合起来才有闭式解? Markov 性负责把递推拆成两步,让每一步只依赖上一步的结果,否则条件分布会挂着越来越长的历史。高斯族对边缘化和条件化封闭,保证两步之后分布形状不变。线性则保证仿射变换后的高斯还是高斯,且参数可以显式算出。三者缺一:丢掉 Markov 性,递推结构不成立;丢掉线性,Chapman–Kolmogorov 积分一般算不出闭式;丢掉高斯,即使积分能算,也无法用有限个参数描述结果。

KF 和一般贝叶斯滤波、粒子滤波是什么关系? 它们用的是同一套方程——Chapman–Kolmogorov 加贝叶斯公式,区别只在于用什么表示概率分布,以及两步操作如何近似执行。KF 用均值和协方差表示分布,两步都解析完成,在线性高斯下是精确的。EKF 和 UKF 仍然用高斯表示,但通过一阶线性化或 sigma 点来近似传播矩,此时高斯只是近似而非精确。粒子滤波用一组带权样本表示分布:Time Update 变成按转移密度对每个粒子采样,Measurement Update 变成按似然给粒子重新加权再重采样,代价是采样噪声和维数灾难,换来的是能处理任意非线性、非高斯、多峰的情形。网格滤波则干脆把状态空间离散化,用直方图表示分布。所以 KF 不是另一种滤波器,而是一般贝叶斯滤波在线性高斯这个特例下退化出来的精确解。


2. 角度二:估计误差传播与最优增益

这一章回答:每个协方差公式怎样从估计误差本身推出来,Kalman 增益为什么最优?

这条路线里不会再出现概率密度,全程只操作估计误差这个随机向量和它的二阶矩。代价是要多假设一件事——估计器的结构必须先给定;收益是多得到一个结论——增益不是抄来的,而是解一个优化问题得到的。第 1 章里 $K_k$ 是从高斯代数中自动掉出来的,这里它是最优化的产物,我们能明确回答“最优”指的是在什么集合里、按什么准则最优。

另一个值得先说的差别是假设强度。下面全部推导只用到零均值和不相关这两类二阶矩性质,从头到尾没有用高斯。所以这条路线的结论比第 1 章更宽:即便噪声不是高斯的,只要一二阶矩对,得到的仍是所有线性无偏估计中均方误差最小的那个,也就是 BLUE。工程上写 KF 代码时脑子里跑的通常就是这套逻辑,$Q$、$R$、$K$ 的调参直觉也从这里来。

2.1 Time Update:预测误差传播

先定义两个误差量。上一历元的后验误差:

\[\tilde x_{k-1}^+=x_{k-1}-\hat x_{k-1}^+\]

当前历元的预测误差:

\[\tilde x_k^-=x_k-\hat x_k^-\]

预测器为什么取 $\hat x_k^-=F_{k-1}\hat x_{k-1}^++B_{k-1}u_{k-1}$,其实也有理由,不是随手写的。设上一历元估计无偏,即 $E[\tilde x_{k-1}^+]=0$,那么对状态方程取期望可知,只有这样取才能保证 $E[\tilde x_k^-]=0$,即无偏性能一路传下去。换句话说,Time Update 的均值公式是“保持无偏”这个要求唯一确定的。

把状态方程和预测式相减,确定性项 $B_{k-1}u_{k-1}$ 完全抵消:

\[\tilde x_k^- = \left(F_{k-1}x_{k-1}+B_{k-1}u_{k-1}+w_{k-1}\right)-\left(F_{k-1}\hat x_{k-1}^++B_{k-1}u_{k-1}\right) = F_{k-1}\tilde x_{k-1}^++w_{k-1}\]

这个式子说明预测误差由两部分组成:旧误差被 $F_{k-1}$ 搬运过来的部分,加上这段区间内新注入的过程噪声。对误差外积取期望,展开成四项:

\[\begin{aligned} P_k^- &=E\!\left[\tilde x_k^-(\tilde x_k^-)^T\right]\\ &=F_{k-1}E\!\left[\tilde x_{k-1}^+(\tilde x_{k-1}^+)^T\right]F_{k-1}^T +F_{k-1}E\!\left[\tilde x_{k-1}^+w_{k-1}^T\right]\\ &\quad+E\!\left[w_{k-1}(\tilde x_{k-1}^+)^T\right]F_{k-1}^T +E\!\left[w_{k-1}w_{k-1}^T\right] \end{aligned}\]

四项逐个看。第一项里 $F_{k-1}$ 和 $F_{k-1}^T$ 分别来自误差外积的左右两侧,这是协方差在线性变换下的标准行为,不是什么特殊构造。第四项就是 $Q_{k-1}$,来自过程噪声自身的协方差。中间两项互为转置,它们为零靠的是假设“过程噪声与上一历元估计误差不相关”,即 $E[\tilde x_{k-1}^+w_{k-1}^T]=0$。于是

\[P_k^-=F_{k-1}P_{k-1}^+F_{k-1}^T+Q_{k-1}\]

如果这个交叉项不为零,记 $M_{k-1}=E[\tilde x_{k-1}^+w_{k-1}^T]$,正确的预测协方差应该是

\[P_k^-=F_{k-1}P_{k-1}^+F_{k-1}^T+F_{k-1}M_{k-1}+M_{k-1}^TF_{k-1}^T+Q_{k-1}\]

这不是纯理论的担心。连续系统离散化时如果过程噪声和量测噪声由同一个物理源驱动、EKF 里 $F$ 依赖于本身含误差的估计值、或者状态中人为塞进了噪声的滤波形式,都会让 $M_{k-1}\neq 0$。第 6 章会具体讨论。

工程含义也就摆在这里:$F_{k-1}P_{k-1}^+F_{k-1}^T$ 是旧的不确定性被动力学搬运和放大的结果,$Q_{k-1}$ 是这一段时间里模型本身不可信而新注入的不确定性。把 $Q$ 调大,$P_k^-$ 变大,滤波器就更愿意听量测的;调小则更相信模型的外推。

2.2 Measurement Update:从线性校正器出发

量测更新这一步,我们不假设任何分布,只限定估计器的结构:它必须是量测的仿射函数,

\[\hat x_k^+=a_k+K_k z_k\]

其中 $a_k$ 和 $K_k$ 待定。再要求它无偏,即对任意真值都有 $E[\tilde x_k^+]=0$。把 $z_k=H_kx_k+v_k$ 代进去取期望,利用 $E[v_k]=0$ 和 $E[x_k]=E[\hat x_k^-]$(预测无偏),可解出 $a_k=(I-K_kH_k)\hat x_k^-$。代回即得

\[\hat x_k^+ = \hat x_k^-+K_k\left(z_k-H_k\hat x_k^-\right)\]

所以“预测值加上增益乘残差”这个熟悉的结构并不是拍脑袋设计的,它是“仿射”加“无偏”两个要求联合决定的唯一形式。定义新息

\[\nu_k=z_k-H_k\hat x_k^-\]

它有一个在后面反复要用的展开:把 $z_k=H_kx_k+v_k$ 代入,得

\[\nu_k=H_k\left(x_k-\hat x_k^-\right)+v_k=H_k\tilde x_k^-+v_k\]

也就是说新息由两部分构成——状态预测错的部分被 $H_k$ 投影到量测空间,加上量测本身的噪声。滤波器无法区分这两者,这正是增益必须做折中的根源。顺带得到新息协方差:

\[S_k=E\!\left[\nu_k\nu_k^T\right]=H_kP_k^-H_k^T+R_k\]

接着推更新后的误差。用真值减去估计式:

\[\tilde x_k^+ = x_k-\hat x_k^--K_k\left(H_k\tilde x_k^-+v_k\right) = \left(I-K_kH_k\right)\tilde x_k^--K_kv_k\]

再取外积期望,四项中两个交叉项因 $E[\tilde x_k^-v_k^T]=0$ 而消失,剩下

\[P_k^+ = \left(I-K_kH_k\right)P_k^-\left(I-K_kH_k\right)^T+K_kR_kK_k^T\]

这就是 Joseph 形式。有一点必须强调:整个推导过程从没有用到 $K_k$ 的具体取值,所以 Joseph 形式对任意增益都成立。取 $K_k=0$ 它退化成 $P_k^+=P_k^-$,符合“不用量测就什么都没学到”;取一个很差的 $K_k$,它照样给出正确的(只是更大的)协方差。这个性质在 2.4 节和工程实现里都很关键。

2.3 用最小均方误差求最优增益

现在才轮到“最优”。我们要在所有线性无偏估计器中挑一个使估计误差最小的 $K_k$。有两条等价路线。

路线一:最小化误差协方差。 先把 Joseph 形式展开,并代入 $S_k=H_kP_k^-H_k^T+R_k$:

\[P_k^+(K_k) = P_k^--K_kH_kP_k^--P_k^-H_k^TK_k^T+K_kS_kK_k^T\]

标量情形下大家会直接对它求导,矩阵情形也一样。用两个常见的矩阵求导恒等式 $\partial\,\mathrm{tr}(AB)/\partial A=B^T$ 和 $\partial\,\mathrm{tr}(ACA^T)/\partial A=2AC$($C$ 对称),对 $\mathrm{tr}(P_k^+)$ 求导并令其为零:

\[\frac{\partial\,\mathrm{tr}(P_k^+)}{\partial K_k} = -2P_k^-H_k^T+2K_kS_k=0 \quad\Longrightarrow\quad K_k=P_k^-H_k^T S_k^{-1}\]

由于 $S_k\succ 0$,二阶条件满足,这是极小值而不是极大值或鞍点。

不过对 $K_k$ 配方能得到比“迹最小”更强的结论。把上面的 $P_k^+(K_k)$ 关于 $K_k$ 配方:

\[P_k^+(K_k) = \underbrace{P_k^--P_k^-H_k^TS_k^{-1}H_kP_k^-}_{\text{与 }K_k\text{ 无关}} + \left(K_k-P_k^-H_k^TS_k^{-1}\right)S_k\left(K_k-P_k^-H_k^TS_k^{-1}\right)^T\]

第二项是关于 $S_k$ 的合同型,恒为半正定,当且仅当 $K_k=P_k^-H_k^TS_k^{-1}$ 时为零。这说明最优增益不只让迹最小,而是让 $P_k^+$ 在矩阵半序意义下最小:对任意其他增益 $K$,都有 $P_k^+(K)\succeq P_k^+(K_k^{\text{opt}})$。推论是,无论你关心的是位置方差、某个线性组合 $a^Tx$ 的方差、协方差的行列式还是任意单调的标量指标,最优的都是同一个 $K_k$。这个增益不依赖于你选哪个代价函数,是相当强的性质。

路线二:正交性原理。 最优估计的一个标志是:残余误差里再也榨不出可用信息,也就是它与所用数据不相关。写成条件即 $E[\tilde x_k^+\nu_k^T]=0$。直接代入 2.2 的两个展开:

\[E\!\left[\tilde x_k^+\nu_k^T\right] = E\!\left[\left(\left(I-K_kH_k\right)\tilde x_k^--K_kv_k\right)\left(H_k\tilde x_k^-+v_k\right)^T\right] = \left(I-K_kH_k\right)P_k^-H_k^T-K_kR_k\]

令其为零,整理得 $P_k^-H_k^T-K_k\left(H_kP_k^-H_k^T+R_k\right)=0$,同样给出

\[K_k=P_k^-H_k^T\left(H_kP_k^-H_k^T+R_k\right)^{-1}\]

两条路线殊途同归。第二条更有几何味道:把随机变量看成内积空间中的向量,最优估计就是把 $x_k$ 投影到由可用数据张成的子空间上,而投影的特征就是残差垂直于该子空间。第 3 章的联合高斯条件化,实际上是同一件事的分布语言版本。

顺便可以看出增益的结构含义:$K_k=P_k^-H_k^TS_k^{-1}$ 中,$P_k^-H_k^T$ 是状态与量测的交叉协方差,$S_k^{-1}$ 是量测自身不确定性的倒数。量测越准($R_k$ 小)增益越大,预测越准($P_k^-$ 小)增益越小,与直觉一致。

2.4 从 Joseph 形式化简到简式

把最优增益代回 Joseph 形式。上一节配方结果的第一项就是答案:

\[P_k^+=P_k^--P_k^-H_k^TS_k^{-1}H_kP_k^-=P_k^--K_kH_kP_k^-=\left(I-K_kH_k\right)P_k^-\]

也可以直接做代数变形,看得更清楚化简发生在哪里。把 Joseph 形式拆开重组:

\[P_k^+ = \left(I-K_kH_k\right)P_k^- - \left[\left(I-K_kH_k\right)P_k^-H_k^T-K_kR_k\right]K_k^T\]

方括号里正是 2.3 中的正交性条件 $E[\tilde x_k^+\nu_k^T]$,在最优增益下为零,于是整项消失,剩下简式。这也就回答了一个常被忽略的问题:简式只在增益取最优值时成立,Joseph 形式对任意增益都成立。

这个区别在工程上有实际后果。凡是动过增益的实现——抗差估计里对可疑量测降权、故障检测后部分屏蔽、衰减记忆滤波里人为放大增益、定点实现里对 $K$ 做了截断或量化——增益都不再是最优值,此时继续用简式算出的 $P_k^+$ 会偏小。协方差偏小意味着滤波器过度自信,后续增益进一步减小,量测越来越推不动状态,最终表现为估计发散但协方差却很漂亮,是很难查的一类问题。

即便增益确实是最优的,浮点实现下两种形式的行为也不同。简式 $(I-K_kH_k)P_k^-$ 是两个一般矩阵相乘,结果在数学上对称、在浮点上不对称,误差会随历元累积,长时间运行后可能出现负的对角元和 Cholesky 分解失败。Joseph 形式右侧是两个合同变换之和,$APA^T$ 和 $KRK^T$ 都天然对称半正定,结构上更难被舍入误差破坏。另外,由于 $P_k^+$ 在最优点处对 $K_k$ 的一阶导为零,$K_k$ 上的小误差只会带来 $P_k^+$ 上的二阶误差,Joseph 形式对增益计算的不精确也更不敏感。代价是计算量明显更大,嵌入式实现里要在稳健性和算力之间权衡,或者干脆走 UDKF、SRKF 这类平方根形式。

最后看协方差是怎么收缩的。利用 $K_kS_kK_k^T=P_k^-H_k^TS_k^{-1}H_kP_k^-$,简式还可以写成

\[P_k^+=P_k^--K_kS_kK_k^T\]

减掉的这一项是半正定的,且秩不超过 $m$,也就是量测的维数。这说明一次量测更新只能在最多 $m$ 个方向上压缩不确定性,而且压缩发生在 $H_k^T$ 的列空间所张成的方向上;落在 $H_k$ 零空间里的状态方向,协方差一点都不会降。GNSS 里这件事随处可见:只有伪距时,接收机高程和钟差沿着某个组合方向几乎不可分离,那个方向上的协方差长期居高不下;只有多普勒观测时,位置状态基本靠 Time Update 撑着。协方差矩阵能不能降下来,取决于量测有没有真正约束到那个方向,而不取决于你喂了多少条观测。

2.5 本角度重点解释的问题

为什么这条路线不需要高斯假设? 因为它只处理一阶矩和二阶矩。整个推导中用到的全部性质是零均值、不相关和线性运算,没有任何地方用到密度函数的具体形状。结论相应地也要弱一点:得到的是所有线性无偏估计中均方误差最小的那个。如果噪声恰好是高斯的,可以进一步证明最优的非线性估计也是线性的,此时 KF 就是全局最优;如果噪声不是高斯的,可能存在更好的非线性估计器,但 KF 仍然是最好的线性估计器。这个区别在有粗差、多路径这类重尾噪声的 GNSS 场景里很实际——它解释了为什么抗差方法有时能明显超过标准 KF。

$Q$、$R$、$K$ 的调参直觉从哪来? 从 $K_k=P_k^-H_k^TS_k^{-1}$ 和 $S_k=H_kP_k^-H_k^T+R_k$ 这两式就能读出来。增大 $Q$ 使 $P_k^-$ 变大,增益变大,滤波器更跟随量测,响应快但噪声大;增大 $R$ 使 $S_k$ 变大,增益变小,输出更平滑但对真实机动反应迟钝。两个极端是:$R\to 0$ 时 $K_k\to$ 量测的伪逆方向,滤波器几乎完全信任量测;$P_k^-\to 0$ 时 $K_k\to 0$,滤波器完全不理量测。真正决定行为的其实是 $Q$ 与 $R$ 的比值,而不是它们的绝对大小——这也是为什么调参时同时放大两者往往看不出变化。

为什么说“最优”必须限定在某个集合里? 因为最优是相对候选集合而言的。这一章把候选集合限定成“$z_k$ 的仿射函数且无偏”,在这个集合里 $K_k$ 是唯一最优解,而且在矩阵半序意义下最优。跳出这个集合,比如允许非线性估计器,或者放弃无偏性接受一点偏差换更小的均方误差(岭回归式的做法),结论就可能不同。明确这一点,才能理解为什么在非高斯、非线性或存在粗差的场景里,KF 会被别的方法超过——它没有失效,只是那些方法不在它的候选集合里。


3. 角度三:随机变量仿射变换与联合高斯条件化

这一章回答:如何只利用高斯随机变量的两个标准结论,最紧凑地推出完整 KF?

前两章都在“推导”,这一章几乎不推导,只是查表。高斯分布族有两条封闭性结论,它们在任何一本概率教材里都有:

  1. 仿射变换:高斯向量的仿射函数仍是高斯向量,均值和协方差按 $E[Ax+b]=AE[x]+b$、$\mathrm{Cov}(Ax+b)=A\,\mathrm{Cov}(x)\,A^T$ 变换。
  2. 条件化:联合高斯向量中,一部分分量在给定另一部分取值后的条件分布仍是高斯,均值和协方差有闭式表达。

Time Update 是第一条的直接应用,Measurement Update 是第二条的直接应用。整章不需要积分、不需要配方、不需要优化,是四条路线里最短的一条。它的另一个价值在后面才显现:因为全程只用到均值、协方差和交叉协方差这几个量,一旦把它们换成用 sigma 点或样本估计出来的版本,同一套公式立刻变成 UKF 或集合卡尔曼滤波,连雅可比矩阵都不需要。

3.1 Time Update:随机变量的仿射变换

把 $x_{k-1}$ 和 $w_{k-1}$ 叠成一个增广向量。以 $z_{1:k-1}$ 为背景条件(下面所有分布都隐含这个条件,为了式子干净就不写了),它服从

\[\begin{bmatrix} x_{k-1}\\ w_{k-1} \end{bmatrix} \sim \mathcal{N}\!\left( \begin{bmatrix} \hat x_{k-1}^+\\ 0 \end{bmatrix}, \begin{bmatrix} P_{k-1}^+ & M_{k-1}\\ M_{k-1}^T & Q_{k-1} \end{bmatrix} \right)\]

这里先不急着把交叉块 $M_{k-1}=\mathrm{Cov}(x_{k-1},w_{k-1})$ 设为零,留着它可以顺手把互相关情形一起推完。状态方程恰好是这个增广向量的仿射函数:

\[x_k = \begin{bmatrix}F_{k-1} & I\end{bmatrix} \begin{bmatrix} x_{k-1}\\ w_{k-1} \end{bmatrix} +B_{k-1}u_{k-1}\]

套用第一条结论。均值:

\[\hat x_k^-=F_{k-1}\hat x_{k-1}^++I\cdot 0+B_{k-1}u_{k-1}=F_{k-1}\hat x_{k-1}^++B_{k-1}u_{k-1}\]

协方差(常数项 $B_{k-1}u_{k-1}$ 不参与):

\[P_k^- = \begin{bmatrix}F_{k-1} & I\end{bmatrix} \begin{bmatrix} P_{k-1}^+ & M_{k-1}\\ M_{k-1}^T & Q_{k-1} \end{bmatrix} \begin{bmatrix}F_{k-1}^T\\ I\end{bmatrix} = F_{k-1}P_{k-1}^+F_{k-1}^T+F_{k-1}M_{k-1}+M_{k-1}^TF_{k-1}^T+Q_{k-1}\]

标准假设下 $M_{k-1}=0$,中间两项消失,回到 $P_k^-=F_{k-1}P_{k-1}^+F_{k-1}^T+Q_{k-1}$。和 2.1 节的结论完全一样,但这里连误差量都不用定义——交叉项不是“被证明为零”,而是联合协方差矩阵里本来就摆在那里的非对角块,你决定它是不是零。这种写法在处理复杂相关结构时会省很多事。

值得留意的是,这一步全程没有用到条件化。Time Update 里没有任何新信息进来,我们只是在计算一个已知随机变量的确定性函数(外加一个独立噪声)的分布,这纯粹是“换变量”。

3.2 Measurement Update:构造状态与量测的联合分布

要用第二条结论,先得有一个联合分布。同样的技巧:把 $x_k$ 和 $v_k$ 叠起来,注意 $(x_k,z_k)$ 是它的仿射函数,

\[\begin{bmatrix} x_k\\ z_k \end{bmatrix} = \begin{bmatrix} I & 0\\ H_k & I \end{bmatrix} \begin{bmatrix} x_k\\ v_k \end{bmatrix}\]

在标准假设下 $\mathrm{Cov}(x_k,v_k)=0$,$x_k$ 的(先验)分布是 $\mathcal{N}(\hat x_k^-,P_k^-)$,$v_k\sim\mathcal{N}(0,R_k)$。再套一次仿射变换法则,就得到

\[\begin{bmatrix} x_k\\ z_k \end{bmatrix} \sim \mathcal{N}\!\left( \begin{bmatrix} \hat x_k^-\\ H_k\hat x_k^- \end{bmatrix}, \begin{bmatrix} P_k^- & P_k^-H_k^T\\ H_kP_k^- & H_kP_k^-H_k^T+R_k \end{bmatrix} \right)\]

三个块逐个核对一下。左上角是状态自身的先验协方差 $P_k^-$。右下角是量测自身的协方差,$\mathrm{Cov}(z_k)=H_kP_k^-H_k^T+R_k=S_k$,也就是新息协方差——它在这里的身份是量测的边缘分布协方差,第 1 章里是归一化常数的协方差,第 2 章里是新息的二阶矩,三者是同一个东西。右上角是交叉协方差:

\[\mathrm{Cov}(x_k,z_k) = E\!\left[\tilde x_k^-\left(H_k\tilde x_k^-+v_k\right)^T\right] = P_k^-H_k^T\]

写成 $P_{xz}=P_k^-H_k^T$、$P_{zz}=S_k$,Kalman 增益就是

\[K_k=P_{xz}P_{zz}^{-1}\]

这个写法比 $P_k^-H_k^TS_k^{-1}$ 更本质:增益的含义是“状态与量测的相关性”除以“量测自身的不确定性”。$H_k$ 只是在线性模型下计算 $P_{xz}$ 的一个中间产物,不是增益定义的一部分。UKF 里没有 $H_k$,但用 sigma 点照样能算出 $P_{xz}$ 和 $P_{zz}$,$K_k=P_{xz}P_{zz}^{-1}$ 原样成立;6.2 节还会回到这一点。

3.3 对联合高斯分布做条件化

现在调用第二条结论。若

\[\begin{bmatrix}a\\b\end{bmatrix} \sim \mathcal{N}\!\left( \begin{bmatrix}\mu_a\\ \mu_b\end{bmatrix}, \begin{bmatrix}\Sigma_{aa}&\Sigma_{ab}\\ \Sigma_{ba}&\Sigma_{bb}\end{bmatrix} \right)\]

则给定 $b$ 后 $a$ 的条件分布是高斯,且

\[E[a\mid b]=\mu_a+\Sigma_{ab}\Sigma_{bb}^{-1}\left(b-\mu_b\right), \qquad \mathrm{Cov}(a\mid b)=\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}\]

这个结论的证明思路值得看一眼,因为它和第 2 章的正交性原理是同一件事。把 $a$ 拆成两部分:

\[a=\underbrace{\left(a-\Sigma_{ab}\Sigma_{bb}^{-1}b\right)}_{\text{与 }b\text{ 不相关}}+\Sigma_{ab}\Sigma_{bb}^{-1}b\]

第一部分与 $b$ 的协方差恰好是 $\Sigma_{ab}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{bb}=0$。在高斯情形下不相关等价于独立,所以给定 $b$ 时第一部分的分布不变,第二部分则成了已知常数。条件均值和条件协方差随之而来。$\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$ 就是 $\Sigma_{bb}$ 在 $\Sigma$ 中的 Schur 补,第 4 章消元时还会再遇到它。

把 3.2 的三个块代进去:

\[E[x_k\mid z_k] = \hat x_k^-+P_k^-H_k^TS_k^{-1}\left(z_k-H_k\hat x_k^-\right) = \hat x_k^-+K_k\nu_k\] \[\mathrm{Cov}(x_k\mid z_k) = P_k^--P_k^-H_k^TS_k^{-1}H_kP_k^- = P_k^--K_kS_kK_k^T = \left(I-K_kH_k\right)P_k^-\]

这就是标准的 Measurement Update。三小节,两条查来的结论,没有一步真正的推导。

有一处容易被绕进去,顺带说清楚:上面条件化的是 $z_k$,而我们真正要的是以 $z_{1:k}$ 为条件的分布。二者一致,是因为所有分布从一开始就隐含地以 $z_{1:k-1}$ 为条件——先验 $\mathcal{N}(\hat x_k^-,P_k^-)$ 已经把历史量测吸收进去了,而 $v_k$ 与历史独立,所以在这个背景条件下再对 $z_k$ 条件化,得到的正是 $p(x_k\mid z_{1:k})$。

3.4 本角度重点解释的问题

为什么“仿射变换”就足以推出 Time Update? 因为 Time Update 根本不是一个推断问题。$x_k$ 由 $x_{k-1}$ 和 $w_{k-1}$ 通过一个确定的仿射映射生成,我们要做的只是把输入的分布推到输出,这是概率论里最基本的操作。没有观测到任何量,也就没有任何需要“反推”的成分。

Measurement Update 为什么必须用条件化? 因为这里做的事情反过来了:我们观测到的是 $z_k$,想知道的是 $x_k$,而 $x_k$ 并不是 $z_k$ 的函数。这属于“已知联合分布和其中一部分的取值,求另一部分的分布”,正是条件化解决的问题。第 1 章的贝叶斯公式、第 2 章的正交投影、这里的条件化,说的是同一件事的三种语言。

$S_k$ 和 $P_k^-H_k^T$ 到底从哪来? 它们不是推导出来的中间量,而是联合协方差矩阵的两个块——在你写下联合分布的那一刻就已经确定了。第 2 章需要绕一圈误差外积才能算出它们,这里则是构造分布时的副产品。这也是为什么这条路线最紧凑。

为什么量测只能沿“状态—量测相关”的方向修正状态? 从 $K_k=P_{xz}P_{zz}^{-1}$ 看得很直接:如果某个状态分量与全部量测的交叉协方差为零,$P_{xz}$ 中对应的那一行就是零,该状态的修正量恒为零,不管量测残差有多大。GNSS 里典型例子是只有伪距观测时的载波相位模糊度:它与伪距量测在模型上就没有关联,$P_{xz}$ 对应行为零,喂再多伪距也不会动它一下。反过来,如果两个状态在 $P_k^-$ 中高度相关,那么修正其中一个必然会带动另一个——接收机高程和钟差就是这样被绑在一起的,这也是为什么单靠 GNSS 很难把高程做准。协方差矩阵的非对角结构,决定了信息在状态之间怎么流动。


4. 角度四:MAP 准则与加权最小二乘

这一章回答:KF 如何等价为递推的二次优化,以及它为何能连接信息滤波、最小二乘和因子图?

前三章都把滤波看成“传播一个分布”。这一章换个立场:把每个历元看成一次优化,滤波器输出的是某个代价函数的极小点。代价函数由一堆加权残差平方相加而成,每个残差对应一条我们相信的信息——上一历元的估计、动力学模型、本历元的每条观测。这个视角在数值上不一定最省,但它是通往 SRIF、滑动窗口、因子图和非线性优化的入口,也是把 KF 和测量平差、Bundle Adjustment 这些东西接上的地方。

先约定记号:$|a|_W^2=a^TWa$。高斯密度的负对数(去掉常数)恰好是这种形式,$W$ 取协方差的逆。

4.1 Time Update:动态约束与上一历元后验

Time Update 在这个视角下是一次变量消元。此刻我们手上有两条信息。一条是上一历元的后验,写成关于 $x_{k-1}$ 的二次代价:

\[\left\|x_{k-1}-\hat x_{k-1}^+\right\|_{\left(P_{k-1}^+\right)^{-1}}^2\]

另一条是动力学模型,它约束的是 $x_{k-1}$ 和 $x_k$ 的关系,残差就是状态方程两边之差:

\[\left\|x_k-F_{k-1}x_{k-1}-B_{k-1}u_{k-1}\right\|_{Q_{k-1}^{-1}}^2\]

两项相加,得到关于 $(x_{k-1},x_k)$ 的联合代价 $J(x_{k-1},x_k)$,它正是这两个变量联合高斯分布的负对数。我们想要的是只关于 $x_k$ 的那部分,所以要把 $x_{k-1}$ 消掉。

概率上消掉一个变量是积分(边缘化),优化上消掉一个变量是取最小(profile)。高斯情形下这两件事给出同一个结果,只差一个与 $x_k$ 无关的常数——因为对高斯做边缘化和做部分最小化,留下的都是同一个二次型。这个等价性只在高斯(或者说二次代价)下成立,非线性、非高斯时两者会分开,滑动窗口边缘化的很多坑就出在这里。

先算均值,有个偷懒办法:$J\geq 0$,而在 $x_{k-1}=\hat x_{k-1}^+$、$x_k=F_{k-1}\hat x_{k-1}^++B_{k-1}u_{k-1}$ 处两项残差同时为零,即 $J=0$。既然全局最小值在这里取到,消元后的代价必然在

\[\hat x_k^-=F_{k-1}\hat x_{k-1}^++B_{k-1}u_{k-1}\]

处取到极小,这就是预测均值。

再算协方差。把联合代价的二次项部分写成信息矩阵(即 Hessian 的一半,按 $J=|\cdot|^2$ 的约定这里直接取 $\Lambda$):

\[\Lambda = \begin{bmatrix} \left(P_{k-1}^+\right)^{-1}+F_{k-1}^TQ_{k-1}^{-1}F_{k-1} & -F_{k-1}^TQ_{k-1}^{-1}\\ -Q_{k-1}^{-1}F_{k-1} & Q_{k-1}^{-1} \end{bmatrix}\]

消去 $x_{k-1}$ 对应在信息矩阵上取 Schur 补,剩下的 $x_k$ 块信息矩阵是

\[\left(P_k^-\right)^{-1} = Q_{k-1}^{-1}-Q_{k-1}^{-1}F_{k-1}\left[\left(P_{k-1}^+\right)^{-1}+F_{k-1}^TQ_{k-1}^{-1}F_{k-1}\right]^{-1}F_{k-1}^TQ_{k-1}^{-1}\]

右边正好是 Woodbury 恒等式的右侧,把它反过来用即得

\[P_k^-=F_{k-1}P_{k-1}^+F_{k-1}^T+Q_{k-1}\]

于是消元后的预测代价就是

\[\left\|x_k-\hat x_k^-\right\|_{\left(P_k^-\right)^{-1}}^2\]

均值和协方差都和前三章一致。

这段推导在因子图语言里有个更好记的说法:Time Update 就是把变量节点 $x_{k-1}$ 从图中消元,消元的代数操作是 Schur 补,代价是原本与 $x_{k-1}$ 相连的所有节点之间会产生新的边(fill-in)。这里只有 $x_k$ 一个邻居,所以填充是平凡的;但在滑动窗口和 MSCKF 里边缘化一个位姿时,它连接的所有特征点之间都会稠密化,这正是边缘化策略必须小心设计的原因。反过来说,如果消元,把所有历元的变量都留在图里,得到的就是批处理平滑或因子图优化——KF 只是它的一个逐步消元的特例。

4.2 Measurement Update:最大后验估计

量测更新在这个视角下是求最大后验估计:

\[\hat x_k^{\mathrm{MAP}} = \arg\max_{x_k} p\left(x_k\mid z_{1:k}\right)\]

用第 1 章的贝叶斯公式 $p(x_k\mid z_{1:k})\propto p(z_k\mid x_k)p(x_k\mid z_{1:k-1})$,两个因子都是高斯,取负对数并丢掉常数($\arg\min$ 不受正比例系数影响,所以 $\tfrac12$ 也一并去掉):

\[\hat x_k^+ = \arg\min_x \left[ \left\|x-\hat x_k^-\right\|_{\left(P_k^-\right)^{-1}}^2 + \left\|z_k-H_kx\right\|_{R_k^{-1}}^2 \right]\]

结构非常干净:第一项要求解不要离预测太远,权是先验精度;第二项要求解能解释观测,权是量测精度。两项都是加权残差平方,地位完全平等——先验在数学上就是一条“伪观测”,只不过它约束的是状态本身而不是状态的函数。GNSS 平差里给某个参数加约束方程,做的就是同一件事。

三个概念需要分清。MAP 是估计准则,说的是“取后验密度最大的点”;最小化负对数后验是求解方法,因为对数单调,取极大和取负对数的极小是一回事;而在标准 KF 里,这个目标函数是严格凸的二次函数,有唯一闭式解,不需要迭代。这一点常被误解:把 KF 写成优化问题不意味着它是非线性优化,只有到了 EKF/IEKF 才真的要迭代。

4.3 法方程与信息形式

目标函数是二次的,求梯度令其为零即可。对上式关于 $x$ 求导:

\[\left(P_k^-\right)^{-1}\left(x-\hat x_k^-\right)-H_k^TR_k^{-1}\left(z_k-H_kx\right)=0\]

整理成法方程:

\[\left[\left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k\right]\hat x_k^+ = \left(P_k^-\right)^{-1}\hat x_k^-+H_k^TR_k^{-1}z_k\]

左边括号里就是后验信息矩阵,也就是目标函数的 Hessian:

\[\left(P_k^+\right)^{-1}=\left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k\]

定义信息向量 $\eta=P^{-1}\hat x$,法方程可以写成更对称的一对:

\[\left(P_k^+\right)^{-1}=\left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k, \qquad \eta_k^+=\eta_k^-+H_k^TR_k^{-1}z_k\]

这就是信息滤波(Information Filter)的 Measurement Update。它的优点一眼可见:量测更新是纯粹的加法,多条独立量测可以任意顺序累加,甚至可以分布式地在不同节点上算完再求和。代价是状态更新变复杂了——Time Update 在信息形式里反而要做矩阵求逆,正好和协方差形式相反。

从这里再往前一步就是 SRIF。把两个加权范数写成堆叠的最小二乘问题:若 $(P_k^-)^{-1}=L_P^TL_P$、$R_k^{-1}=L_R^TL_R$,则目标等价于

\[\min_x \left\| \begin{bmatrix} L_P\left(x-\hat x_k^-\right)\\ L_R\left(z_k-H_kx\right) \end{bmatrix} \right\|^2\]

对这个堆叠矩阵做 QR 分解,就直接得到后验的信息平方根因子,全程不碰 $P$ 也不碰 $\Lambda$。这正是站内那篇讲 SRIF 的文章的出发点,也说明为什么 SRIF 更像“递推形式的加权最小二乘”而不是“换了个变量的 KF”。

选哪种形式,看维度。状态维 $n$、量测维 $m$:增益形式要求解一个 $m\times m$ 的 $S_k$,信息形式要求解一个 $n\times n$ 的信息矩阵。$m\gg n$ 时(比如多系统多频点,几十条观测更新十来个状态)信息形式更划算;$m\ll n$ 时(比如 PPP 里几十个状态、每次只来几条观测)增益形式更划算。信息形式还有一个协方差形式做不到的能力:完全无先验可以直接写成 $\Lambda=0$,而 $P^-=\infty$ 在浮点里没法表示,这也是初始化阶段信息形式更自然的原因。

4.4 从信息形式恢复 Kalman 增益形式

最后把信息形式变回大家熟悉的增益形式。用 Woodbury 矩阵恒等式

\[\left(A+UCV\right)^{-1}=A^{-1}-A^{-1}U\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}\]

取 $A=(P_k^-)^{-1}$、$U=H_k^T$、$C=R_k^{-1}$、$V=H_k$,得到

\[P_k^+ = \left[\left(P_k^-\right)^{-1}+H_k^TR_k^{-1}H_k\right]^{-1} = P_k^--P_k^-H_k^T\left(R_k+H_kP_k^-H_k^T\right)^{-1}H_kP_k^-\]

括号里是 $S_k$,于是 $P_k^+=P_k^–K_kH_kP_k^-=(I-K_kH_k)P_k^-$,且

\[K_k=P_k^-H_k^T\left(H_kP_k^-H_k^T+R_k\right)^{-1}\]

均值那一半要多一步代数,这里补上(第 1 章引用的就是这一段)。由法方程,

\[\hat x_k^+=P_k^+\left[\left(P_k^-\right)^{-1}\hat x_k^-+H_k^TR_k^{-1}z_k\right] = \left(I-K_kH_k\right)\hat x_k^-+\left(I-K_kH_k\right)P_k^-H_k^TR_k^{-1}z_k\]

关键是化简 $z_k$ 前面那一坨。注意 $H_kP_k^-H_k^T=S_k-R_k$,

\[\left(I-K_kH_k\right)P_k^-H_k^T = P_k^-H_k^T-K_k\left(S_k-R_k\right) = P_k^-H_k^T-P_k^-H_k^T+K_kR_k = K_kR_k\]

所以 $(I-K_kH_k)P_k^-H_k^TR_k^{-1}=K_k$,代回即得

\[\hat x_k^+=\hat x_k^--K_kH_k\hat x_k^-+K_kz_k=\hat x_k^-+K_k\left(z_k-H_k\hat x_k^-\right)\]

四条路线到这里全部走到了同一组公式。顺带一提,$(I-K_kH_k)P_k^-H_k^T=K_kR_k$ 这个恒等式本身也有意思:它说明后验状态与量测残差之间的协方差结构,正好让 Joseph 形式里的交叉项抵消,和 2.4 节看到的是同一件事。

4.5 本角度重点解释的问题

为什么先验和量测都表现为加权残差? 因为高斯密度的负对数就是加权残差平方。$-\log\mathcal{N}(a;\mu,\Sigma)=\tfrac12|a-\mu|_{\Sigma^{-1}}^2+\text{const}$,这个形式对先验和似然一视同仁。所以“把先验当成一条伪观测”不是类比,而是字面意义上成立的:你完全可以把 $\hat x_k^-$ 当作一组直接观测状态的虚拟量测,量测矩阵为 $I$,噪声协方差为 $P_k^-$,塞进一次普通的加权最小二乘里,结果与 KF 完全一致。

权矩阵为什么是协方差的逆? 同样来自指数上的二次型。直观上也说得通:方差越大的量越不可信,权重越低,而“越不可信”的量化正是方差的倒数。工程上这里有两个常见错误——把标准差而不是方差填进 $R$,以及给 $R$ 只填对角而忽略量测间的相关性。前者直接把权算错,后者会重复计数信息,两种都会让滤波器过度自信。

KF、递推最小二乘、信息滤波、因子图是什么关系? 它们是同一个二次优化问题的不同切法。递推最小二乘(RLS)是 KF 在 $F=I$、$Q=0$ 时的特例,也就是被估计量不随时间变化的情形。信息滤波和 KF 是同一组方程在 $(\Lambda,\eta)$ 与 $(P,\hat x)$ 两套坐标下的表示,SRIF 和 SRKF 则各自再取一次平方根。因子图把每个残差项画成一个因子、每个状态画成一个节点,Measurement Update 是往图里加因子,Time Update 是从图里消元;如果不急着消元而是保留整个窗口,就是滑动窗口估计;全部保留就是批处理平滑。KF 之所以能做成 $O(1)$ 存储的递推,就是因为它在每个历元都立刻把旧变量消元了。

为什么线性高斯下 MAP、后验均值、MMSE 三者重合? 高斯分布单峰对称,众数就是均值,所以 MAP 等于后验均值;而后验均值本身就是最小均方误差估计(这是 MMSE 的一般结论,与是否高斯无关);再加上线性高斯下后验均值恰好是量测的线性函数,它同时落进了第 2 章那个“线性无偏估计”的候选集合里,于是又和 BLUE 重合。四个准则指向同一个解,这就是四条路线能汇合的深层原因。反过来,只要跳出线性高斯,它们立刻分家:重尾分布的众数和均值不在一处,非线性模型下 MMSE 估计也不再是线性的。


5. 四种推导如何汇合到同一组公式

四条路线走的是完全不同的数学,落点却是同一组式子。先用一张表对齐它们的核心操作:

推导角度 Time Update 的核心操作 Measurement Update 的核心操作 最适合解释
贝叶斯递推 概率边缘化 先验乘似然并归一化 KF 的概率本质
误差传播与最优增益 传播估计误差协方差 最小化后验误差协方差 工程公式、$Q/R/K$ 的意义
仿射变换与联合高斯条件化 高斯变量仿射变换 联合高斯条件分布 最紧凑的数学推导
MAP 与加权最小二乘 动态因子消元 二次代价闭式求解 信息滤波、最小二乘和因子图

5.1 中间量的一一对应

$P_k^-$ 在四条路线里分别是:预测分布 $p(x_k\mid z_{1:k-1})$ 的协方差参数;预测误差 $\tilde x_k^-$ 的二阶矩;增广向量 $(x_{k-1},w_{k-1})$ 经仿射变换后的协方差;联合代价消去 $x_{k-1}$ 后所剩二次型的信息矩阵之逆。四种说法背后是同一个矩阵,但强调的东西不同——第 1 章强调它描述的是一个分布的宽度,第 2 章强调它描述的是估计误差的大小,第 4 章强调它是 Schur 补的产物。

$S_k$ 分别是:贝叶斯归一化常数 $p(z_k\mid z_{1:k-1})$ 的协方差;新息 $\nu_k=H_k\tilde x_k^-+v_k$ 的二阶矩;联合分布 $(x_k,z_k)$ 中量测的边缘协方差;法方程经 Woodbury 变换后出现在括号里的那个矩阵。第三种说法最能解释它的名字——它就是量测自身的不确定性,只不过量测的不确定性同时来自状态预测不准和量测噪声两部分。

$K_k$ 分别是:高斯乘积配方后自动出现的系数;使后验误差协方差在矩阵半序下最小的最优线性校正系数;联合高斯条件化公式里的 $\Sigma_{ab}\Sigma_{bb}^{-1}$;MAP 法方程经 Woodbury 恒等式改写后的解系数。最有推广价值的是第三种写法 $K_k=P_{xz}P_{zz}^{-1}$,因为它不含 $H_k$,在 UKF、EnKF 乃至存在噪声互相关的情形下都原样成立。

Joseph 形式与简式 的关系,只有第 2 章讲得清楚——因为只有那条路线里增益是个自由变量。其余三条路线一上来就在最优增益处工作,自然直接得到简式,看不到“非最优增益下会怎样”这个维度。

5.2 四条路线的分工

如果只能记一条,记第 3 章:三小节、两条查表结论,是最快能重建出全部公式的路径。但另外三条各有它们不可替代的地方。

想理解 KF 到底在做什么、它和粒子滤波之类的东西是什么关系,看第 1 章。只有概率视角能说清 KF 不是一个孤立的算法,而是通用贝叶斯滤波在线性高斯特例下的精确解。

要写代码、调参数、排查滤波器发散,看第 2 章。$Q$ 和 $R$ 的比值如何决定跟随性、为什么改过增益就必须用 Joseph 形式、协方差为什么在某些方向上死活降不下来,这些问题的答案都在误差传播那一章。

要做大规模问题——PPP、多站网解、滑动窗口、GNSS/视觉/惯导融合,看第 4 章。信息形式、稀疏结构、变量消元、鲁棒核函数,都是从二次优化视角自然长出来的,而从协方差视角看过去它们几乎不可见。

5.3 为什么它们必然汇合

四条路线用了四个不同的最优性准则:后验分布(贝叶斯)、最小均方误差(BLUE)、条件期望(投影)、最大后验(MAP)。在一般情形下这四个准则给出四个不同的答案。它们在这里重合,是因为线性高斯这个组合太特殊:高斯单峰对称使众数等于均值,于是 MAP 等于后验均值;后验均值天然是 MMSE 最优解;线性模型使后验均值恰好是量测的线性函数,于是它又落进了线性无偏估计的候选集合,与 BLUE 重合。

这也意味着一旦离开线性高斯,四条路线就会分家,而且分家的方式决定了各类扩展算法的形态:坚持贝叶斯准则、放弃闭式解,得到粒子滤波;坚持高斯表示、用矩近似,得到 UKF;坚持 MAP 准则、用迭代求解,得到 IEKF 和因子图优化。下一章就沿着这个思路,看看模型变复杂之后各条路线分别通向哪里。


6. 推导之后需要补充的边界条件

前面五章都在标准假设下工作。真实系统很少这么听话,这一章把几类常见的偏离拿出来,说明公式该怎么改、四条路线分别给出什么样的推广。

6.1 噪声不独立时

处理相关性有一个统一的办法,就是第 3 章那条最本质的增益公式:

\[K_k=P_{xz}P_{zz}^{-1}, \qquad P_k^+=P_k^--K_kP_{zz}K_k^T\]

只要重新算对交叉协方差 $P_{xz}=\mathrm{Cov}(x_k,\nu_k)$ 和新息协方差 $P_{zz}=\mathrm{Cov}(\nu_k)$,公式的形状就不用变。下面几种情形都可以照这个套路处理。

状态误差与过程噪声相关。 记 $M_{k-1}=E[\tilde x_{k-1}^+w_{k-1}^T]\neq 0$,3.1 节已经推过:

\[P_k^-=F_{k-1}P_{k-1}^+F_{k-1}^T+F_{k-1}M_{k-1}+M_{k-1}^TF_{k-1}^T+Q_{k-1}\]

什么时候会出现?连续系统离散化时若过程噪声在区间内同时影响了上一历元的估计和本区间的转移;EKF 中 $F$ 由含误差的估计值线性化而来;或者把有色噪声的滤波器状态并入状态向量而没有正确处理耦合。

过程噪声与量测噪声相关。 设 $E[w_{k-1}v_k^T]=C_k\neq 0$。由于 $\tilde x_k^-$ 含 $w_{k-1}$,它与 $v_k$ 不再无关,于是

\[P_{zz}=H_kP_k^-H_k^T+H_kC_k+C_k^TH_k^T+R_k, \qquad P_{xz}=P_k^-H_k^T+C_k\] \[K_k=\left(P_k^-H_k^T+C_k\right)\left(H_kP_k^-H_k^T+H_kC_k+C_k^TH_k^T+R_k\right)^{-1}\]

忽略 $C_k$ 的后果是双向的:$P_{zz}$ 算错会让新息检验的门限失准,$P_{xz}$ 算错会让增益偏离最优。GNSS 里这种相关常来自同一个物理源同时进入动力学和观测,例如把接收机钟建模成随机游走的同时又用同一段钟差改正观测。

量测噪声在时间上相关(有色噪声)。 标准 KF 要求 $v_k$ 是白的,但 GNSS 里多路径、对流层残差、星历改正误差、SBAS 差分改正的残余误差全都是有色的,相关时间从几十秒到几十分钟不等。用白噪声近似会严重高估信息量——同一个偏差被当成 $N$ 次独立观测重复计数,协方差降得比实际快得多,典型现象是位置输出很平滑但真实误差远大于 $3\sigma$ 包络。正确做法是状态增广:把误差建成一阶 Gauss–Markov 过程并放进状态向量,让滤波器自己估计它;或者用量测差分法消掉相关部分(代价是差分后的噪声反而变成相关的,需要配合处理)。

不同量测之间相关($R$ 非对角)。 常见来源包括:双差观测因共用参考星和参考站而强相关;同一历元不同频点共用同一套电离层、对流层模型误差;同一接收机的所有观测共用一个钟差改正误差。把这些当成独立观测,同样是重复计数信息,结果是过度自信。相应地,$Q$ 也会有同样的问题——多接收机或多系统钟差共享同一块晶振时,$Q$ 的钟差块必然带互相关项,而且随手补的相关系数很容易让 $Q$ 失去正半定性,这个话题在站内那篇讲多钟差过程噪声互相关的文章里有详细讨论。

6.2 非线性模型时

真实模型通常是

\[x_k=f\left(x_{k-1},u_{k-1}\right)+w_{k-1}\] \[z_k=h\left(x_k\right)+v_k\]

此时四条路线各自延伸出一类方法,而且延伸的方向和它们各自的“瓶颈”直接相关。

误差传播视角 → EKF。 在当前估计点做一阶泰勒展开,取雅可比 $F_{k-1}=\partial f/\partial x$、$H_k=\partial h/\partial x$,其余公式照抄。它的问题也来自这里:协方差是按线性化后的模型传播的,线性化误差不进入 $P$,所以 EKF 系统性地低估不确定性,强非线性或初值偏差大时容易发散。均值传播也有偏——$E[f(x)]\neq f(E[x])$。

仿射变换与条件化视角 → UKF、EnKF。 这条路线只需要 $P_{xz}$ 和 $P_{zz}$,不需要 $H_k$。用一组确定性采样的 sigma 点过非线性函数,再由样本统计出这两个矩,就得到 UKF;用蒙特卡洛样本代替 sigma 点,就是集合卡尔曼滤波。均值和协方差的传播精度到二阶(EKF 只到一阶),而且不用求雅可比,模型里有查表、分段函数这类不可导环节时特别有用。

MAP 视角 → IEKF、因子图优化。 EKF 只在预测点线性化一次,相当于对非线性最小二乘做了一步高斯牛顿。既然目标函数摆在那里,为什么不多迭代几次?每次在最新估计点重新线性化,就是迭代 EKF——它本质上是用高斯牛顿求解 MAP 问题。再往前,加 Levenberg–Marquardt 阻尼、加 Huber 或 Cauchy 鲁棒核抑制粗差、保留多个历元变量做滑动窗口,就进入了因子图优化的领域。这条路线的优势是能明确写出目标函数,从而知道自己在优化什么、收敛到哪里;EKF 则说不清它在优化什么。

贝叶斯视角 → 粒子滤波、高斯和滤波。 当后验根本不是单峰高斯时(多路径导致的多假设、整周模糊度的离散性、城市峡谷里的多模态位置分布),前三条路线的高斯表示从根上就不成立。此时只能回到最一般的贝叶斯递推,用样本或高斯混合来表示分布。代价是计算量和维数灾难,所以实践中多用于低维子问题,或者与 KF 混合(Rao-Blackwellized 粒子滤波:非线性的少数几维用粒子,条件线性的其余维度用 KF)。

6.3 数值实现时

不要显式求逆。 $K_k=P_k^-H_k^TS_k^{-1}$ 在代码里不应该写成求逆再相乘,而应该解线性方程组 $S_kX=H_kP_k^-$,用 Cholesky 或 LDLT 分解。这样既快又稳,还顺带得到 $S_k$ 是否正定的判据——如果 Cholesky 分解失败,说明 $P$ 或 $R$ 已经坏了,此时应该报警而不是继续算。

量测维数高时用顺序更新。 若 $R_k$ 是对角的,$m$ 条量测可以一条一条地更新,每次只需要一个标量除法,完全避免 $m\times m$ 求逆。$R_k$ 非对角时,先对 $R_k$ 做 Cholesky 分解并对量测方程做白化变换,再顺序更新。GNSS 里可见卫星数动辄二三十颗,这个技巧很实用,而且天然适合边解码边更新的流式处理。

协方差的对称性和正定性要主动维护。 每次更新后做一次 $P\leftarrow(P+P^T)/2$ 几乎没有代价,能挡掉大部分舍入误差累积。再加上对角元非负检查、条件数监控,以及必要时的特征值截断。更彻底的做法是走平方根形式:UDKF 维护 $P=UDU^T$,SRKF 维护 $P=SS^T$,SRIF 维护 $\Lambda=R^TR$,它们在结构上就不可能产生负方差,代价是实现复杂一些。

注意状态量纲带来的条件数问题。 GNSS 状态向量里位置是米、钟差是米、钟漂是米每秒、对流层是毫米、模糊度是周,数值尺度可能差好几个数量级。$P$ 的条件数因此很大,直接算容易丢精度。做法是对状态做无量纲化缩放,或者至少在数值上把量级拉近。这个问题在双精度下常被掩盖,一旦移植到单精度或定点的嵌入式平台就会暴露。

用一致性检验监控滤波器。 第 1 章证明了新息服从 $\mathcal{N}(0,S_k)$,这给了一组现成的在线检查手段:归一化新息平方 $\mathrm{NIS}_k=\nu_k^TS_k^{-1}\nu_k$ 应服从自由度为 $m$ 的卡方分布,长期均值应接近 $m$;新息序列应当是白的,自相关明显不为零说明模型有色误差没建对;有真值时还可以算 NEES 检查状态误差与 $P$ 是否匹配。NIS 系统性偏大说明 $Q$ 或 $R$ 给小了(过度自信),偏小说明给大了(过于保守)。这套检验既是调参依据,也是运行时的健康监测,比只盯着定位残差有用得多。

返回文章列表 返回首页