从贝叶斯估计到卡尔曼滤波(详细推导)(下)

承接上文 从贝叶斯估计到卡尔曼滤波(详细推导)(上)

---

二.高斯加权分布的多维积分运算

1.第一章总结

现在我们针对上一章中的重要结果做一个简单回顾与总结。

我们假设在一个动态过程，其形式满足下式
$$x_n=f_{n-1}(x_{n-1})+u_n+v_{n-1}(42)$$
并且观测过程满足
$$z_n=h_n(x_n)+w_n(43)$$
基于到$t_{n-1}$时刻为止的所有观测值，对$t_n$时刻状态向量的预测估计可以写为
$${\hat{x}}_{n|n-1}=\int \left[f_{n-1}(x_{n-1})+u_n+v_{n-1}\right]p(x_{n-1}|z_{1:n-1})d{x}_{n-1}(44)$$
对该状态向量的协方差矩阵的预测估计可以由下式得到
$$\begin{gathered} P_{n|n-1}^{xx}=\int \left[f_{n-1}(x_{n-1})+u_n-{\hat{x}}_{n|n-1}\right]{\left[f_{n-1}(x_{n-1})+u_n-{\hat{x}}_{n|n-1}\right]}^{T}p(x_{n-1}|z_{1:n-1})d{x}_{n-1} \\ +Q(45) \end{gathered}$$
另外，对观测向量$z_n$、$z_n$的协方差矩阵，$z_n$与$x_n$的互协方差矩阵的预测估计由下式得到

$$\begin{gathered} {\hat{z}}_{n|n-1}=\int h_n(x_n)p(x_n|z_{1:n-1})d{x}_n \\ +\int w_{n}p(x_n|z_{1:n-1})d{x}_n(46) \end{gathered}$$

$$\begin{gathered} P_{n|n-1}^{zz}=\int [h_n(x_n)-{\hat{z}}_{n|n-1}][h_n(x_n)-{\hat{z}}_{n|n-1}{]}^T \\ \times p(x_n|z_{1:n-1})d{x}_n+R(47) \end{gathered}$$

$$\begin{gathered} P_{n|n-1}^{xz}=\int [x_n-{\hat{x}}_{n|n-1}][h_n(x_n)-{\hat{z}}_{n|n-1}{]}^T \\ \times p(x_n|z_{1:n-1})d{x}_n(48) \end{gathered}$$

其中，
$$R:=\int w_n{w}_n^{T}p(x_n|z_{1:n-1})d{x}_n(49)$$
$$Q:=\int v_{n-1}{v}_{n-1}^{T}p(x_{n-1}|z_{1:n-1})d{x}_{n-1}(50)$$

2.修正卡尔曼滤波器更新方程的推导过程

在这一节，我们假设所有的条件密度都是高斯的，并在次假设下重新推导卡尔曼的更新方程。

根据贝叶斯准则，可以得到
$$p(x_n|z_{1:n})=\frac{p(x_n,z_n|z_{1:n-1})}{p(z_n)}(51)$$
现在定义一个联合向量
$$q:=\left[\begin{array}{c}{x}_n\\ z_n\end{array}\right](52)$$
则式(51)可被重写为
$$p(x_n|z_{1:n})=\frac{p(q_n|z_{1:n-1})}{p(z_n)}(53)$$
我们知道，多元高斯分布的一般形式可写为
$$N(t;s,\Sigma ):=\frac{1}{\sqrt{(2\pi {)}^n|\Sigma |}}exp\left\{-\frac{1}{2}(t-s{)}^T{\Sigma }^{-1}(t-s)\right\}(54)$$

则式(52)中的概率密度可以表示为
$$p(x_n|z_{1:n})=N(x_n;{\hat{x}}_{n|n},P_{n|n}^{xx})(55)$$
$$p(q_n|z_{1:n-1})=N(q_n;{\hat{q}}_{n|n-1},P_{n|n}^{qq})(56)$$
$$p(z_n)=N(z_n;{\hat{z}}_{n|n-1},P_{n|n-1}^{zz})(57)$$

忽略系数$-\frac{1}{2}$,则式(53)的左式中，后验概率密度的解析表达式中的指数项的指数可以写为
$$\begin{gathered} d:=(x_n-{\hat{x}}_{n|n}{)}^T[P_{n|n}^{xx}{]}^{-1}(x_n-{\hat{x}}_{n|n}) \\ =x_n^T[P_{n|n}^{xx}{]}^{-1}{x}_n-x_n^T[P_{n|n}^{xx}{]}^{-1}{\hat{x}}_{n|n} \\ -{\hat{x}}_{n|n}^T[P_{n|n}^{xx}{]}^{-1}{x}_n+{\hat{x}}_{n|n}^T[P_{n|n}^{xx}{]}^{-1}{\hat{x}}_{n|n}(58) \end{gathered}$$
同样，式(53)的右式中，后验概率密度的解析表达式中的指数项的指数可以写为
$$\begin{gathered} g:=(q_n-{\hat{q}}_{n|n-1}{)}^T[P_{n|n-1}^{qq}{]}^{-1}(q_n-{\hat{q}}_{n|n-1}) \\ -(z_n-{\hat{z}}_{n|n-1}{)}^T[P_{n|n-1}^{zz}{]}^{-1}(z_n-{\hat{z}}_{n|n-1})(59) \end{gathered}$$
使用式(52),则上式可以写为
$$\begin{gathered} g={\left[\begin{array}{c}{x}_n-{\hat{x}}_{n|n-1}\\ z_n-{\hat{z}}_{n|n-1}\end{array}\right]}^T{\left[\begin{array}{cc}{P}_{n|n-1}^{xx}& P_{n|n-1}^{xz}\\ P_{n|n-1}^{zx}& P_{n|n-1}^{zz}\end{array}\right]}^{-1}\left[\begin{array}{c}{x}_n-{\hat{x}}_{n|n-1}\\ z_n-{\hat{z}}_{n|n-1}\end{array}\right] \\ -(z_n-{\hat{z}}_{n|n-1}{)}^T[P_{n|n-1}^{zz}{]}^{-1}(z_n-{\hat{z}}_{n|n-1})(60) \end{gathered}$$
定义
$$\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]:={\left[\begin{array}{cc}{P}_{n|n-1}^{xx}& P_{n|n-1}^{xz}\\ P_{n|n-1}^{zx}& P_{n|n-1}^{zz}\end{array}\right]}^{-1}(61)$$
其中，
$$C_{11}=[P_{n|n-1}^{xx}-P_{n|n-1}^{xz}(P_{n|n-1}^{zz}{)}^{-1}{P}_{n|n-1}^{zx}{]}^{-1}(62)$$
$$C_{12}=-C_{11}{P}_{n|n-1}^{xz}(P_{n|n-1}^{zz}{)}^{-1}(63)$$
$$C_{21}=-C_{22}{P}_{n|n-1}^{zx}(P_{n|n-1}^{xx}{)}^{-1}(64)$$
$$C_{22}=P_{n|n-1}^{zz}-P_{n|n-1}^{zx}(P_{n|n-1}^{xx}{)}^{-1}{P}_{n|n-1}^{xz}(65)$$
将式(61)-式(65)带入式(60)，可得
$$g=x_n^T{C}_{11}{x}_n+x_n^T[-C_{11}{\hat{x}}_{n|n-1}+C_{12}(z_n-{\hat{z}}_{n|n-1})]+...(66)$$
由于式(66)与式(58)等价，比较两式的第一项，我们可以得到
$$P_{n|n}^{xx}=P_{n|n-1}^{xx}-P_{n|n-1}^{xz}(P_{n|n-1}^{zz}{)}^{-1}{P}_{n|n-1}^{zx}(67)$$
比较第二项，可以得到
$$\begin{gathered} [P_{n|n}^{xx}{]}^{-1}{\hat{x}}_{n|n}=C_{11}{\hat{x}}_{n|n-1}+C_{12}(z_n-{\hat{z}}_{n|n-1}) \\ =[P_{n|n-1}^{xx}-P_{n|n-1}^{xz}(P_{n|n-1}^{zz}{)}^{-1}{P}_{n|n-1}^{zx}{]}^{-1}{\hat{x}}_{n|n-1} \\ +[P_{n|n-1}^{xx}-P_{n|n-1}^{xz}(P_{n|n-1}^{zz}{)}^{-1}{P}_{n|n-1}^{zx}{]}^{-1}{P}_{n|n-1}^{zx} \\ \times (P_{n|n-1}^{zz}{)}^{-1}(z_n-{\hat{z}}_{n|n-1}) \\ =[P_{n|n}^{xx}{]}^{-1}{\hat{x}}_{n|n-1} \\ +[P_{n|n}^{xx}{]}^{-1}{P}_{n|n-1}^{zx}(P_{n|n-1}^{zz}{)}^{-1}(z_n-{\hat{z}}_{n|n-1})(68) \end{gathered}$$
进一步推导上式，$P_{n|n}^{xx}$左乘上式，可以得到
$${\hat{x}}_{n|n}={\hat{x}}_{n|n-1}+P_{n|n-1}^{zx}(P_{n|n-1}^{zz}{)}^{-1}(z_n-{\hat{z}}_{n|n-1})(69)$$
定义卡尔曼增益为
$$K_n:=P_{n|n-1}^{zx}(P_{n|n-1}^{zz}{)}^{-1}(70)$$
现在，我们通过式(69)和式(67)就得到了卡尔曼滤波器的更新公式！！
$$x_{n|n}={\hat{x}}_{n|n-1}+K_n(z_n-{\hat{z}}_{n|n-1})(71)$$
$$P_{n|n}^{xx}=P_{n|n-1}^{xx}-K_n(P_{n|n-1}^{zz})K_n^T(72)$$

注意，我们没有对观测方程和动态方程做出任何线性假设。总结来说，卡尔曼滤波器的更新公式，可以由高斯分布直接推导得到，而不需要做其他任何假设。

3.高斯密度下贝叶斯点预测的积分运算

在上一节我们在高斯假设下，重新推导了卡尔曼滤波器的更新公式。在本节，我们再在高斯假设下，进一步完善卡尔曼预测公式。

我们假设$v_n$是独立于$x_n$与$z_n$的零均值高斯随机噪声，服从高斯分布
$$v_{n-1}\sim p(x_n|x_{n-1})=N(x_n-[f_n(x_{n-1})+u_n];0,Q)=N(v_{n-1};0,Q)(73)$$
类似的，我们假设在观测方程中，$w_n$同样是独立于$x_n$与$z_n$的零均值高斯随机噪声，服从高斯分布
$$w_n\sim p(z_n|x_n)=N(z_n-h_n(x_n);0,R)=N(w_n;0.R)(74)$$
根据贝叶斯准则，我们知道
$$p(x_n|z_{1:n-1})\propto p(z_n|x_n)p(x_n|x_{n-1})(75)$$
由于两个高斯分布的概率密度函数的乘积仍然是服从高斯分布的，因此，$p(x_n|z_{1:n-1})$也是高斯的。假设所有的概率密度函数都是高斯分布的，则
$$p(x_n|z_{1:n-1})=N(x_n;{\hat{x}}_{n|n-1},P_{n|n-1}^{xx})(76)$$
$$p(x_{n-1}|z_{1:n-1})=N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})(77)$$

将式(77)带入式(44)、式(45)，并将式(76)带入式(46)、式(47)、式(48)，得到
$${\hat{x}}_{n|n-1}=\int f_{n-1}(x_{n-1})N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})d{x}_{n-1}+u(n)(78)$$

$$\begin{gathered} P_{n|n-1}^{xx}=\int \left[f_{n-1}(x_{n-1})+u_n-{\hat{x}}_{n|n-1}\right]{\left[f_{n-1}(x_{n-1})+u_n-{\hat{x}}_{n|n-1}\right]}^T \\ \times N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})d{x}_{n-1} \\ +Q(79) \end{gathered}$$

$${\hat{z}}_{n|n-1}=\int h_n(x_n)N(x_n;{\hat{x}}_{n|n-1},P_{n|n-1}^{xx})d{x}_n(80)$$

$$\begin{gathered} P_{n|n-1}^{zz}=\int [h_n(x_n)-{\hat{z}}_{n|n-1}][h_n(x_n)-{\hat{z}}_{n|n-1}{]}^T \\ \times N(x_n;{\hat{x}}_{n|n-1},P_{n|n-1}^{xx})d{x}_n+R(81) \end{gathered}$$

$$\begin{gathered} P_{n|n-1}^{xz}=\int [x_n-{\hat{x}}_{n|n-1}][h_n(x_n)-{\hat{z}}_{n|n-1}{]}^T \\ \times N(x_n;{\hat{x}}_{n|n-1},P_{n|n-1}^{xx})d{x}_n(82) \end{gathered}$$

其中，
$$R:=\int w_n{w}_n^{T}N(w_n;0,R)d{x}_n(83)$$
$$Q:=\int v_{n-1}{v}_{n-1}^{T}N(v_{n-1};0,Q)d{x}_{n-1}(84)$$

终于，我们得到了在高斯假设下的预测公式。并且，离 线性卡尔曼滤波器越来越近！！

总结回顾下，在本章，我们在第一章结果的基础上，重新推导了在概率密度函数服从高斯分布的前提下，预测与更新步骤中的方程。目前，我们还未得到最后可以在工程中实际可用的解析解。

三.线性卡尔曼滤波器

在这一章，我们将在前面几章的基础上，推导线性卡尔曼滤波器的预测与更新公式。

1.线性状态方程

回顾前面几章，我们为了最后得到解析解，已经做出了一系列假设。我们假设系统服从一阶马尔可夫过程；我们假设使用最小均方准则作为代价函数；我们假设概率密度函数服从高斯分布…现在，我们再做出假设：假设状态方程是线性的，则动态方程可以重写为
$$x_n=F{x}_{n-1}+u_n+v_{n-1}(85)$$
其中，$F$是与时间有关的状态转移矩阵。将式(85)带入式(78)，得到
$$\begin{gathered} {\hat{x}}_{n|n-1}=F\int x_{n-1}N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})d{x}_{n-1}+u(n) \\ =F{\hat{x}}_{n-1|n-1}+u_n(86) \end{gathered}$$
我们再将式(85)带入式(79)，得到

$$\begin{gathered} P_{n|n-1}^{xx}=\int F\left[x_{n-1}-{\hat{x}}_{n|n-1}\right]{\left[x_{n-1}-{\hat{x}}_{n|n-1}\right]}^T{F}^T \\ \times N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})d{x}_{n-1}+Q \\ =F\left\{\int \left[x_{n-1}-{\hat{x}}_{n|n-1}\right]{\left[x_{n-1}-{\hat{x}}_{n|n-1}\right]}^T \\ \times N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})d{x}_{n-1}\right\}{F}^T+Q(87) \end{gathered}$$
由于，
$$\begin{gathered} P_{n-1|n-1}^{xx}=\int [x_{n-1}-{\hat{x}}_{n|n-1}][x_{n-1}-{\hat{x}}_{n|n-1}{]}^T \\ \times N(x_{n-1};{\hat{x}}_{n-1|n-1},P_{n-1|n-1}^{xx})d{x}_{n-1}(88) \end{gathered}$$
因此，
$$P_{n|n-1}^{xx}=F{P}_{n-1|n-1}^{xx}{F}^T+Q(89)$$

2.线性观测方程

假设观测方程也是线性的，因此
$$z_n=H{x}_n+w_n(90)$$
其中，$H$是与时间和$x_n$无关的观测矩阵。
因此，我们可以写为
$$h(x_n)=H{x}_n(91)$$
假设，观测噪声$w_n$是零均值过程，则式(80),式(81)可以更新为
$${\hat{z}}_{n|n-1}=H{\hat{x}}_{n|n-1}(92)$$
$$P_{n|n-1}^{zz}=H{P}_{n|n-1}^{xx}{H}^T+R(93)$$
最后，将式(90),式(91)带入式(82)，得到
$$\begin{gathered} P_{n|n-1}^{xz}=\int [x_n-{\hat{x}}_{n|n-1}][z_n-{\hat{z}}_{n|n-1}{]}^T \\ \times N(x_n;{\hat{x}}_{n|n-1},P_{n|n-1}^{xx})d{x}_n \\ =\int [x_n-{\hat{x}}_{n|n-1}][H{x}_n+w_n-H{\hat{x}}_{n|n-1}{]}^T \\ \times N(x_n;{\hat{x}}_{n|n-1},P_{n|n-1}^{xx})d{x}_n \\ =P_{n|n-1}^{xx}{H}^T(94) \end{gathered}$$

3.总结

在观测方程与状态方程都满足严格线性的假设时，我们终于得到了可以使用的线性卡尔曼滤波器公式！！
针对卡尔曼滤波器的预测更新方程,现在做一个归纳总结。

Step1.状态向量的预测
$${\hat{x}}_{n|n-1}=F{\hat{x}}_{n-1|n-1}+u_n(95)$$
$$P_{n|n-1}^{xx}=F{P}_{n-1|n-1}^{xx}{F}^T+Q(96)$$

Step2.观测向量的预测
$${\hat{z}}_{n|n-1}=H{\hat{x}}_{n|n-1}(97)$$
$$P_{n|n-1}^{zz}=H{P}_{n|n-1}^{xx}{H}^T+R(98)$$
$$P_{n|n-1}^{xz}=P_{n|n-1}^{xx}{H}^T(99)$$

Step3.卡尔曼滤波器更新
$$K_n:=P_{n|n-1}^{zx}(P_{n|n-1}^{zz}{)}^{-1}(100)$$
$$x_{n|n}={\hat{x}}_{n|n-1}+K_n(z_n-{\hat{z}}_{n|n-1})(101)$$
$$P_{n|n}^{xx}=P_{n|n-1}^{xx}-K_n(P_{n|n-1}^{zz})K_n^T(102)$$

不知道有没有读者觉得上文中卡尔曼增益的推导公式有些晦涩难懂？别担心。在线性假设的前提下，我们用一种更加简单的方式推导一下卡尔曼增益公式。

我们已知状态向量的滤波公式为${\hat{x}}_{n|n}={\hat{x}}_{n|n-1}+K_n(z_n-{\hat{z}}_{n|n-1})$。记滤波误差为${\tilde{x}}_{n|n}={\hat{x}}_{n|n}-x_n$,其中$x_n$为$n$时刻的真值。将式(101)代入，可得
$$\begin{array}{rl}{\tilde{x}}_{n|n}& ={\hat{x}}_{n|n-1}+K_n(z_n-{\hat{z}}_{n|n-1})-x_n\\ & ={\hat{x}}_{n|n-1}+K_n(H{x}_n+v_n)-K_{n}H{\hat{x}}_{n|n-1}-x_n\\ & =(I-K_{n}H){\hat{x}}_{n|n-1}-(I-K_{n}H)x_n+K_n{v}_n\\ & =(I-K_{n}H)({\hat{x}}_{n|n-1}-x_n)+K_n{v}_n\\ & =(I-K_{n}H){\tilde{x}}_{n|n-1}+K_n{v}_n\end{array}$$
其中，$v_n$为$n$时刻的观测噪声，${\tilde{x}}_{n|n-1}$为预测误差。

由$P_{n|n}$与${\tilde{x}}_{n|n}$的关系可知
$$\begin{array}{rl}{P}_{n|n}& =E\left\{{\tilde{x}}_{n|n}{\tilde{x}}_{n|n}^T\right\}\\ & =E\left\{((I-K_{n}H){\tilde{x}}_{n|n-1}+K_n{v}_n)((I-K_{n}H){\tilde{x}}_{n|n-1}+K_{n}v{)}^T\right\}\\ & =(I-K_{n}H)P_{n|n-1}(I-K_{n}H{)}^T+K_n{R}_n{K}_n^T\\ & =P_{n|n-1}-P_{n|n-1}{H}^T{K}_n^T-K_{n}H{P}_{n|n-1}+K_{n}H{P}_{n|n-1}{H}^T{K}_n^T+K_n{R}_n{K}_n^T\end{array}$$

现在需要找到使${\tilde{x}}_{n|n}$最小情况下的卡尔曼增益$K_n$的表达式，最直接的方法就是求偏导。
$$\frac{\partial {\Vert {\tilde{x}}_{n|n}\Vert }_F^2}{\partial K_n^T}=\frac{\partial tr({\tilde{x}}_{n|n}{\tilde{x}}_{n|n}^T)}{\partial K_n^T}\approx \frac{\partial tr(E\left\{{\tilde{x}}_{n|n}{\tilde{x}}_{n|n}^T\right\})}{\partial K_n^T}=\frac{\partial tr(P_{n|n})}{\partial K_n^T}$$

取上式可得
$$\begin{array}{r}\frac{\partial tr\left\{P_{n|n}\right\}}{\partial K_n^T}=-P_{n|n-1}{H}^T-(H{P}_{n|n-1}{)}^T+2K_{n}H{P}_{n|n-1}{H}^T+2K_n{R}_n\end{array}$$
令导数等于0，可得
$$\begin{gathered} -P_{n|n-1}{H}^T-(H{P}_{n|n-1}{)}^T+2K_{n}H{P}_{n|n-1}{H}^T+2K_n{R}_n=0 \\ {K}_n=P_{n|n-1}{H}^T(H{P}_{n|n-1}{H}^T+R_n{)}^{-1} \end{gathered}$$
上式与式(100)等价。

此外，我们还得到了滤波协方差矩阵的另外两种表达式
$$P_{n|n}^{xx}=P_{n|n-1}^{xx}-P_{n|n-1}^{xz}(P_{n|n-1}^{zz}{)}^{-1}{P}_{n|n-1}^{zx}$$
$$P_{n|n}^{xx}=(I-K_{n}H)P_{n|n-1}(I-K_{n}H{)}^T+K_n{R}_n{K}_n^T$$

尾声

在第一章，我们从贝叶斯公式入手，先推导了在贝叶斯框架下，针对点估计问题，最优估计值的表达形式，得到了最优估计值与后验概率密度的关系式。接着，进一步考虑更实际的场景。由于实际问题中，我们是针对一个连续的动态过程做估计。因此，为了简化问题，假设该过程为一阶马尔可夫过程，并以此假设为前提，引出了过程方程和观测方程，同时，得到了相邻时刻的后验概率密度之间的递推公式。并发现，该递推过程，可大致分为预测与更新两个步骤。最后，再次利用最优估计值是后验概率密度的一阶矩的结论，分别推导了在预测与更新两个步骤中，后验概率密度的一阶矩与二阶矩的公式。

在第二章，我们做出了进一步的假设。假设概率密度函数都是高斯的，并在此假设下，修正了之前在预测与更新步骤中，所涉及的核心公式。并发现，在高斯假设下求得的更新公式与在最小均方准则和一阶马尔可夫过程假设下求得的更新公式是一样的。

在第三章，我们又一次做出了假设。假设动态方程与观测方程都是严格线性的，在多重假设下，我们进一步修正第二章中的预测步与更新步的核心公式，终于得到了大家最熟悉的卡尔曼滤波器。

不知有没有帮到各位读者？总之，完结撒花~