Kalman 滤波最优状态估计特性的几何解读


看到一些论文和资料, 将 Kalman 滤波置于 Hilbert Space 描述下进行几何解读 5-6. 这些几何解释让人有更直观和深刻的理解, 观点也可以更高屋建瓴. 我们借鉴这些资料, 试图更详细地对 Kalman Filter 进行几何意义的解读. 写下来也为了自己学习和理解.

我们先对要讨论的卡尔曼滤波相关内容进行简要回顾.

I. 系统方程

离散时间系统方程

x t = F t x t − 1 + B t u t + w t (I-1) \mathbf{x}_{t} = \mathbf{F}_{t} \mathbf{x}_{t-1} + \mathbf{B}_{t}\mathbf{u}_{t} + \mathbf{w}_{t} \tag{I-1} xt=Ftxt1+Btut+wt(I-1)

描述了系统的当前时刻状态 x t \mathbf{x}_t xt 由上一时刻状态 x t − 1 \mathbf{x}_{t-1} xt1 演变的动力学, 也称为预测方程. 其中

t t t 是离散时间步, t = 0 , 1 , 2 , … t=0,1,2,\dots t=0,1,2,

x t \mathbf{x}_t xt n n n 维状态向量;

F t \mathbf{F}_t Ft 为系统的状态转移矩阵;

u t \mathbf{u}_t ut 为当前控制输入;

B t \mathbf{B}_t Bt 为控制输入矩阵;

w t \mathbf{w}_t wt n n n 维零均值过程白噪声, 即
E [ w t w l T ] = Q t δ t l = { Q t ( t = l ) 0 ( t ≠ l ) (I-2) \mathbf{E}[\mathbf{w}_t \mathbf{w}_l^{\rm \scriptsize T}] =\mathbf{Q}_t \delta_{tl} =\left \{ \begin{aligned} &\mathbf{Q}_t & (t=l) \\ &0 & (t \neq l) \end{aligned} \right. \tag{I-2} E[wtwlT]=Qtδtl={Qt0(t=l)(t=l)(I-2)
其中 Q t \mathbf{Q}_t Qt 为正定的协方差矩阵.


测量方程

y t = H t x t + v t (I-3) \mathbf{y}_t = \mathbf{H}_t \mathbf{x}_t +\mathbf{v}_t \tag{I-3} yt=Htxt+vt(I-3)

其中

y t \mathbf{y}_t yt m m m 维观测向量;

H t \mathbf{H}_t Ht 观测矩阵, 描述状态向量到观测向量的变换;

v t \mathbf{v}_t vt n n n 维零均值的观测白噪声, 即
E [ v t v t T ] = R t δ t l = { R t ( t = l ) 0 ( t ≠ l ) (I-4) \mathbf{E}[\mathbf{v}_t \mathbf{v}_t^{\rm \scriptsize T}] = \mathbf{R}_t \delta_{tl}= \left \{ \begin{aligned} &\mathbf{R}_t & (t=l) \\ &\mathbf{0} & (t \neq l) \end{aligned} \right. \tag{I-4} E[vtvtT]=Rtδtl={Rt0(t=l)(t=l)(I-4)
其中 R t \mathbf{R}_t Rt 为正定的协方差矩阵.

需要说明:

- 我们本文的分析都建立在线性系统的基础上;

- 一步时间内的控制输入提供一个确定的偏置 (后面分析时再讨论);

- 随机变量 x t \mathbf{x}_t xt w t \mathbf{w}_t wt v t \mathbf{v}_t vt 之间独立不相关.

在系统模型的基础上, 继续回顾对系统状态进行估计的卡尔曼滤波方法.

II. Kalman Filter

卡尔曼滤波算法分为两个阶段, 分别是预测阶段和更新阶段1.


预测阶段

x ^ t ∣ t − 1 = F t x ^ t − 1 ∣ t − 1 + B t u t (II-1) \hat{\mathbf{x}}_{t \vert t-1} = \mathbf{F}_t \hat{\mathbf{x}}_{t-1 \vert t-1} + \mathbf{B}_t \mathbf{u}_t \tag{II-1} x^tt1=Ftx^t1∣t1+Btut(II-1)

P t ∣ t − 1 = F t P t − 1 ∣ t − 1 F t T + Q t (II-2) \mathbf{P}_{t \vert t-1} = \mathbf{F}_{t} \mathbf{P}_{t-1\vert t-1} \mathbf{F}_t^{\rm \scriptsize T} + \mathbf{Q}_t \tag{II-2} Ptt1=FtPt1∣t1FtT+Qt(II-2)

其中

x ^ t ∣ t − 1 \hat{\mathbf{x}}_{t \vert t-1} x^tt1 是基于预测方程 (II-1) 对未知的真实状态向量 x t \mathbf{x}_t xt 的预测, 是对真实状态向量的先验估计 (Priori Estimation).

P t ∣ t − 1 \mathbf{P}_{t \vert t-1} Ptt1 是预测状态的误差协方差矩阵.


更新阶段

x ^ t ∣ t = x ^ t ∣ t − 1 + K t ( y t − H t x ^ t ∣ t − 1 ) (II-3) \hat{\mathbf{x}}_{t \vert t} = \hat{\mathbf{x}}_{t \vert t-1} + \mathbf{K}_{t} \left(\mathbf{y}_t - \mathbf{H}_t \hat{\mathbf{x}}_{t\vert t-1}\right) \tag{II-3} x^tt=x^tt1+Kt(ytHtx^tt1)(II-3)

P t ∣ t = P t ∣ t − 1 − K t H t P t ∣ t − 1 (II-4) \mathbf{P}_{t\vert t} = \mathbf{P}_{t\vert t-1} - \mathbf{K}_t \mathbf{H}_t \mathbf{P}_{t\vert t-1} \tag{II-4} Ptt=Ptt1KtHtPtt1(II-4)

K t = P t ∣ t − 1 H t T ( H t P t ∣ t − 1 H t T + R t ) − 1 (II-5) \mathbf{K}_t = \mathbf{P}_{t\vert t-1} \mathbf{H}_t^{\rm \scriptsize T}\left(\mathbf{H}_t \mathbf{P}_{t\vert t-1} \mathbf{H}_t^{\rm\scriptsize T} + \mathbf{R}_t\right)^{-1} \tag{II-5} Kt=Ptt1HtT(HtPtt1HtT+Rt)1(II-5)

其中

x ^ t ∣ t \hat{\mathbf{x}}_{t \vert t} x^tt 称为修正的状态估计, 是对状态预测值 x ^ t ∣ t − 1 \hat{\mathbf{x}}_{t \vert t-1} x^tt1 的更新, 是对真实状态向量的后验估计 (Posteriori Estimation).

y ~ t = y t − H t x ^ t ∣ t − 1 \tilde{\mathbf{y}}_t=\mathbf{y}_t - \mathbf{H}_t \hat{\mathbf{x}}_{t\vert t-1} y~t=ytHtx^tt1 为新息 (Innovation Variable), 也叫测量残差, 是实际传感器测量值与传感器预测值之间的差.

K t \mathbf{K}_t Kt 为卡尔曼滤波增益 (Kalman Gain), 是该算法的核心.

H t P t ∣ t − 1 H t T + R t \mathbf{H}_t \mathbf{P}_{t\vert t-1} \mathbf{H}_t^{\rm\scriptsize T} + \mathbf{R}_t HtPtt1HtT+Rt 称为新息协方差矩阵.
P t ∣ t \mathbf{P}_{t\vert t} Ptt 为更新状态的误差协方差矩阵.

卡尔曼滤波的预测阶段的重点是动力学建模, 需要对系统的动力学特性进行准确刻画, 得到的状态预测值 x ^ t ∣ t − 1 \hat{\mathbf{x}}_{t \vert t-1} x^tt1 作为更新阶段的先验信息.

预测阶段相对明确和直观, 没有融合与优化. 预测阶段只作为更新阶段的先验输入对待.

卡尔曼滤波多传感器信息的融合、最优状态估计等特性都体现在更新阶段, 卡尔曼滤波器的特性是在更新阶段实现的.

从最优化的角度, 卡尔曼滤波是最优估计, 是一种最小均方误差估计器 (Minimum Mean-Square Error Estimator).

下面我们从几何正交投影角度解释卡尔曼滤波的优化特性、新息以及更新过程的几何意义等.

III. Hilbert Space

我们已经熟悉了欧式空间, 维度是确定和有限的, 内积计算也简单. 我们将要解释处理的卡尔曼滤波器的状态向量或者观测向量的维度并不一定限于某一值, 不同的系统维度变化不同; 另外, 状态变量和观测数据都是被噪声/扰动污染的随机向量. 甚至, 更深入的研究利用 Kalman Filter 进行时间-空间过程的最优估计 (无穷维). 为了得到一致和有力的数学工具, 引入可以描述无限维内积空间的 Hilbert Space 的定义. 当然, 此处我们在处理卡尔曼滤波数据时, 还是基于有限维内积空间中的情形进行考虑.

Hilbert Space2

A Hilbert space is a vector space H \mathcal{H} H with an inner product ⟨ f , g ⟩ \langle f, g \rangle f,g such that the norm defined by
∥ f ∥ = ⟨ f , f ⟩ (III-1) \Vert f\Vert =\sqrt{\langle f,f \rangle} \tag{III-1} f=f,f (III-1)
turns H \mathcal{H} H into a complete metric space. If the metric defined by the norm is not complete, then H \mathcal{H} H is instead known as an inner product space.

已知 n n n 维实的随机向量
x = [ x 1 x 2 ⋮ x n ] ∈ R n (III-2) \mathbf{x} = \begin{bmatrix}x_{1}\\ x_{2}\\ \vdots \\ x_{n} \end{bmatrix} \in \mathbb{R}^n \tag{III-2} x= x1x2xn Rn(III-2)
其中每个元素 x i x_{i} xi 为随机变量 (随机的标量).

假设 { y i } i ∈ N \{ \mathbf{y}_i \}_{i\in\mathbb{N}} {yi}iN 序列组成了 (III-2) 形式的 n n n 维随机向量的集合. 该集合内任意元素的所有线性组合构成了一个随机向量空间 H \mathcal{H} H. 我们在该空间上构造内积.

假设 x \mathbf{x} x y \mathbf{y} y 是随机向量空间 H \mathcal{H} H 上的任意元素, x \mathbf{x} x y \mathbf{y} y 都为 n n n 维随机向量, 则两者的内积定义为
⟨ x , y ⟩ = E [ x T y ] = E [ ∑ i = 1 n x i y i ] = T r ( E [ x y T ] ) (III-3) \langle \mathbf{x},\mathbf{y}\rangle = \mathbf{E}[\mathbf{x}^{\rm\scriptsize T} \mathbf{y}] = \mathbf{E}\left[\sum_{i=1}^{n} {x_i} {y_i}\right] = {\rm Tr}(\mathbf{E}[\mathbf{x}\mathbf{y}^{\rm\scriptsize T}]) \tag{III-3} x,y=E[xTy]=E[i=1nxiyi]=Tr(E[xyT])(III-3)
其中 Tr 即为矩阵的迹 (Trace). 由内积导出的随机向量 ∥ x ∥ \Vert \mathbf{x} \Vert x 的范数的定义
∥ x ∥ = E [ x T x ] = E [ x 1 2 + x 2 2 + ⋯ + x n 2 ] = E [ x 1 2 ] + E [ x 2 2 ] + ⋯ + E [ x n 2 ] = T r a c e ( E [ x x T ] ) (III-4) \begin{aligned} \Vert x \Vert &= \sqrt{\mathbf{E}[\mathbf{x}^{\rm \scriptsize T} \mathbf{x}]} \\ &= \sqrt{\mathbf{E}[x_1^2 + x_2^2+\dots+x_n^2]}\\ &= \sqrt{\mathbf{E}[x_1^2] + \mathbf{E}[x_2^2]+\dots + \mathbf{E}[x_n^2]}\\ &=\sqrt{\rm{Trace}(\mathbf{E}[\mathbf{x}\mathbf{x}^{\rm\scriptsize T}])} \end{aligned} \tag{III-4} x=E[xTx] =E[x12+x22++xn2] =E[x12]+E[x22]++E[xn2] =Trace(E[xxT]) (III-4)
关于随机向量空间 H \mathcal{H} H 的完备性, 此处只能给个大概的参考描述. 随机向量 x \mathbf{x} x 中独立随机变量 x i x_{i} xi 处于完备概率空间 (完备化得到), 由参考文献可知3, 随机向量 x \mathbf{x} x 所处内积空间是完备的. 故装备了以上 (III-3) 内积定义的随机向量空间 H \mathcal{H} H 就是本文需要构建的那个 Hilbert Space.

下文讨论的距离就是式 (III-4) 定义的范数.

IV. Orthogonal Porjection

Orthogonal Decomposition4

If Y \mathcal{Y} Y is a closed subspace of Hilbert space H \mathcal{H} H, then we have:
H = Y ⊕ Y ⊥ (IV-1) \mathcal{H} = \mathcal{Y} \oplus \mathcal{Y}^{\bot} \tag{IV-1} H=YY(IV-1)
where Y ⊥ = { z ∈ H : ⟨ z , y ⟩ = 0 , ∀ y ∈ Y } \mathcal{Y}^{\bot} = \{\mathbf{z} \in \mathcal{H}: \langle \mathbf{z}, \mathbf{y} \rangle = 0, \forall \mathbf{y} \in \mathcal{Y}\} Y={zH:z,y=0,yY}.

In other words, for any x ∈ H \mathbf{x} \in \mathcal{H} xH, x \mathbf{x} x can be expressed uniquely as
x = y + z \mathbf{x} = \mathbf{y} + \mathbf{z} x=y+z
where y ∈ Y \mathbf{y} \in \mathcal{Y} yY and z ∈ Y ⊥ \mathbf{z} \in \mathcal{Y}^{\bot} zY.
Moreover, y \mathbf{y} y and z \mathbf{z} z are the unique elements of Y \mathcal{Y} Y and Y ⊥ \mathcal{Y}^{\bot} Y whose distance to x \mathbf{x} x is smallest.


Fig.1 正交分解示意图

为了说明正交分解得到的 y \mathbf{y} y z \mathbf{z} z 是各自空间内到达 x \mathbf{x} x 距离最近的点, 我们只需要借用勾股定理 (毕达哥拉斯定理), 如 Fig.1 所示.

简单证明:

M \mathcal{M} M 上任一不同于正交分解所得向量 y \mathbf{y} y 的另一向量 y ′ \mathbf{y}^{'} y, 则由勾股定理
∥ x − y ′ ∥ = ∥ x − y ∥ + ∥ y ′ − y ∥ > ∥ x − y ∥ \Vert \mathbf{x} - \mathbf{y}^{'}\Vert = \Vert \mathbf{x} - \mathbf{y}\Vert + \Vert \mathbf{y}^{'} - \mathbf{y}\Vert > \Vert \mathbf{x} - \mathbf{y}\Vert xy=xy+yy>xy
M ⊥ \mathcal{M}^{\bot} M 上的 z \mathbf{z} z 同理. 得证. 有这个简单证明, 可以知道正交分解是唯一的.

为下面说明方便, 定义记号如下:

- 假设向量 x ∈ H \mathbf{x} \in \mathcal{H} xH 和一个闭的线性子空间 Y ⊂ H \mathcal{Y} \subset \mathcal{H} YH.

- 沿着 Y ⊥ \mathcal{Y}^{\bot} Y 投影于 Y \mathcal{Y} Y 上的正交投影子 (Orthogonal Projector) 记为 P Y ( ⋅ ) \mathcal{P}_{\mathcal{Y}}(\cdot) PY().

- 向量 x \mathbf{x} x 向子空间 Y \mathcal{Y} Y 进行正交投影变换得到投影向量记为 x Y \mathbf{x}_{\mathcal{Y}} xY, 即
x Y = P Y ( x ) (IV-2) \mathbf{x}_{\mathcal{Y}} = \mathcal{P}_{\mathcal{Y}} (\mathbf{x}) \tag{IV-2} xY=PY(x)(IV-2)
- 另记距离平方 σ Y ( x ) = ∥ x − x Y ∥ 2 \sigma_{\mathcal{Y}}(\mathbf{x}) = \Vert \mathbf{x} - \mathbf{x}_{\mathcal{Y}} \Vert^{2} σY(x)=xxY2.

V. Minimum Variance Estimate

由测量方程
y = H x + v (V-1) \mathbf{y} = \mathbf{H} \mathbf{x} + \mathbf{v} \tag{V-1} y=Hx+v(V-1)
估计系统状态 x ^ = K y \hat{\mathbf{x}} = \mathbf{K} \mathbf{y} x^=Ky 称为基于测量 y \mathbf{y} y 的线性估计.

Minimum Variance Linear Estimate5

Let y \mathbf{y} y and x \mathbf{x} x be random vectors. Assume that [ E [ y y T ] ] − 1 \left[\mathbf{E} [\mathbf{y} \mathbf{y}^{\rm \scriptsize T} ]\right]^{−1} [E[yyT]]1 exists. The linear estimate x ^ \hat{\mathbf{x}} x^ of x {\mathbf{x}} x based on y \mathbf{y} y that minimizes ∥ x − x ^ ∥ 2 \Vert {\mathbf{x}}-\hat{\mathbf{x}} \Vert^2 xx^2 is
x ^ = E [ x y T ]   [ E [ y y T ] ] − 1   y (V-2) \hat{\mathbf{x}} = \mathbf{E}[\mathbf{x} \mathbf{y}^{\rm \scriptsize T} ] \, \left[\mathbf{E} [\mathbf{y} \mathbf{y}^{\rm \scriptsize T} ]\right]^{−1} \,\mathbf{y} \tag{V-2} x^=E[xyT][E[yyT]]1y(V-2)

式 (V-2) 得到的 x ^ \hat{\mathbf{x}} x^ 使得下式成立
x ^ = arg ⁡ min ⁡ m ∈ Y ∥ x − m ∥ 2 = arg ⁡ min ⁡ m ∈ Y   E [ ( x − m ) T   ( x − m ) ] (V-3) \begin{aligned} \hat{\mathbf{x}} &= \underset{\mathbf{m}\in \mathcal{Y}}{\arg\min} { \Vert {\mathbf{x}}-{\mathbf{m}} \Vert^2} \\ &= \underset{\mathbf{m}\in \mathcal{Y}}{\arg\min}\,\sqrt{\mathbf{E}\left[(\mathbf{x} - {\mathbf{m}})^{\rm\scriptsize T} \, (\mathbf{x} - {\mathbf{m}}) \right]} \end{aligned} \tag{V-3} x^=mYargminxm2=mYargminE[(xm)T(xm)] (V-3)
根据前面 “Orthogonal Decomposition” 定理可知, 极值条件 (V-3) 等价于估计误差 x − x ^ \mathbf{x}-\hat{\mathbf{x}} xx^ 与测量向量 y \mathbf{y} y 垂直 (内积为零), 即
⟨ x − x ^ , y ⟩ = E [ ( x − x ^ ) T   y ] = 0 & x ^ = K y (V-4) \langle \mathbf{x}-\hat{\mathbf{x}}, \mathbf{y} \rangle = \mathbf{E}\left[(\mathbf{x}-\hat{\mathbf{x}})^{\rm\scriptsize T}\, \mathbf{y}\right] = 0 \quad \& \quad \hat{\mathbf{x}}= \mathbf{K}\mathbf{y}\tag{V-4} xx^,y=E[(xx^)Ty]=0&x^=Ky(V-4)
此处 y \mathbf{y} y 张成子空间 Y \mathcal{Y} Y. 最小方差估计 x ^ \hat{\mathbf{x}} x^ 就是让状态 x \mathbf{x} x 正交投影到测量 y \mathbf{y} y 张成的子空间上得到 x Y \mathbf{x}_{\mathcal{Y}} xY, 即
x Y = x ^ (V-5) \mathbf{x}_{\mathcal{Y}} = \hat{\mathbf{x}} \tag{V-5} xY=x^(V-5)
∥ x − x ^ ∥ \Vert \mathbf{x} - \hat{\mathbf{x}}\Vert xx^ 视作估计误差, 记为 ϵ \boldsymbol{\epsilon} ϵ.

需要说明,

- 无偏估计情况下, 最小均方误差估计式 (V-3) 和最小方差估计等价;

- Kalman Filter 初值的设定以及迭代过程, 使其是无偏估计.

Minimum Variance Linear Estimate of a Linear Function5

The minimum variance linear estimate of a linear function of x \mathbf{x} x is equal to the linear function of the minimum variance estimate of x \mathbf{x} x. In other words given a matrix T \mathbf{T} T (a linear function), the minimum variance estimate of T x \mathbf{T} \mathbf{x} Tx is
T x ^ = T   E [ x y T ]   [ E [ y y T ] ] − 1 y (V-6) {\mathbf{T}}\hat{\mathbf{x}} = \mathbf{T} \, \mathbf{E}[\mathbf{x}\mathbf{y}^{\rm\scriptsize T}]\,\left[\mathbf{E}[\mathbf{y} \mathbf{y}^{\rm \scriptsize T} ]\right]^{−1} \mathbf{y} \tag{V-6} Tx^=TE[xyT][E[yyT]]1y(V-6)

简单证明:

T x \mathbf{T}\mathbf{x} Tx 看作整体替换 “Minimum Variance Linear Estimate” 中式 (V-2) 中的 x \mathbf{x} x,
( T x ) ^ = E [ ( T x ) y T ]   [ E [ y y T ] ] − 1   y = T   E [ x y T ]   [ E [ y y T ] ] − 1   y = T x ^ (V-7) \begin{aligned} \hat{(\mathbf{T}\mathbf{x})} &= \mathbf{E}[(\mathbf{T}\mathbf{x}) \mathbf{y}^{\rm \scriptsize T} ] \, \left[\mathbf{E} [\mathbf{y} \mathbf{y}^{\rm \scriptsize T} ]\right]^{−1} \,\mathbf{y}\\ &= \mathbf{T} \, \mathbf{E}[\mathbf{x} \mathbf{y}^{\rm \scriptsize T} ] \, \left[\mathbf{E} [\mathbf{y} \mathbf{y}^{\rm \scriptsize T} ]\right]^{−1} \,\mathbf{y}\\ &= \mathbf{T} \hat{\mathbf{x}} \end{aligned} \tag{V-7} (Tx)^=E[(Tx)yT][E[yyT]]1y=TE[xyT][E[yyT]]1y=Tx^(V-7)
线性转换后的最小方差估计是最小方差估计的线性变换, 即
P Y ( T x ) = T P Y ( x ) (V-8) \mathcal{P}_{\mathcal{Y}} (\mathbf{T}\mathbf{x}) = \mathbf{T} \mathcal{P}_{\mathcal{Y}} (\mathbf{x}) \tag{V-8} PY(Tx)=TPY(x)(V-8)

推论 V-1:

[ E [ y y T ] ] − 1 [\mathbf{E}[\mathbf{y}\mathbf{y}^{\rm\scriptsize T}]]^{−1} [E[yyT]]1 存在情况下, 任意向量 x \mathbf{x} x 和线性变换 T \mathbf{T} T 满足
- T P Y ( x ) ∈ Y \mathbf{T} \mathcal{P}_{\mathcal{Y}} (\mathbf{x})\in \mathcal{Y} TPY(x)Y;

- T x − T P Y ( x ) ∈ Y ⊥ \mathbf{T}\mathbf{x} - \mathbf{T} \mathcal{P}_{\mathcal{Y}} (\mathbf{x})\in \mathcal{Y}^{\bot} TxTPY(x)Y.

证明: 因为 P Y ( T x ) ∈ Y \mathcal{P}_{\mathcal{Y}} (\mathbf{T}\mathbf{x})\in \mathcal{Y} PY(Tx)Y 及式 (V-8);
\qquad 因为 T x − P Y ( T x ) ∈ Y ⊥ \mathbf{T}\mathbf{x} - \mathcal{P}_{\mathcal{Y}} (\mathbf{T}\mathbf{x})\in \mathcal{Y}^{\bot} TxPY(Tx)Y 及式 (V-8).

在这一章结束前, 我们补充定义沿着 Y \mathcal{Y} Y 投影于 Y ⊥ \mathcal{Y}^{\bot} Y 上的正交投影子为 P Y ⊥ ( ⋅ ) \mathcal{P}_{\mathcal{Y}^{\bot}}(\cdot) PY().
由 Orthogonal Decomposition 定义, ∀ x ∈ H \forall \mathbf{x} \in \mathcal{H} xH
x = P Y ( x ) + P Y ⊥ ( x ) (V-9) \mathbf{x} = \mathcal{P}_{\mathcal{Y}}(\mathbf{x}) + \mathcal{P}_{\mathcal{Y}^{\bot}}(\mathbf{x}) \tag{V-9} x=PY(x)+PY(x)(V-9)
考虑线性变换 T \mathbf{T} T P Y ⊥ ( ⋅ ) \mathcal{P}_{\mathcal{Y}^{\bot}}(\cdot) PY() 的影响
T x = P Y ( T x ) + P Y ⊥ ( T x ) ⇒ ( V − 8 ) P Y ⊥ ( T x ) = T x − T P Y ( x ) ⇒ P Y ⊥ ( T x ) = T [ x − P Y ( x ) ] ⇒ ( V − 9 ) P Y ⊥ ( T x ) = T P Y ⊥ ( x ) (V-10) \begin{aligned} \mathbf{T} \mathbf{x} & = \mathcal{P}_{\mathcal{Y}}(\mathbf{T}\mathbf{x}) + \mathcal{P}_{\mathcal{Y}^{\bot}}(\mathbf{T}\mathbf{x})\\ \overset{\rm{(V-8)}}{\Rightarrow} \quad \mathcal{P}_{\mathcal{Y}^{\bot}}(\mathbf{T}\mathbf{x}) &= \mathbf{T} \mathbf{x} - \mathbf{T} \mathcal{P}_{\mathcal{Y}}(\mathbf{x}) \\ {\Rightarrow} \quad \mathcal{P}_{\mathcal{Y}^{\bot}}(\mathbf{T}\mathbf{x}) &= \mathbf{T}\left[ \mathbf{x} - \mathcal{P}_{\mathcal{Y}}(\mathbf{x}) \right]\\ \overset{\rm{(V-9)}}{\Rightarrow} \quad \mathcal{P}_{\mathcal{Y}^{\bot}}(\mathbf{T}\mathbf{x}) &= \mathbf{T} \mathcal{P}_{\mathcal{Y}^{\bot}}(\mathbf{x}) \end{aligned}\tag{V-10} Tx(V8)PY(Tx)PY(Tx)(V9)PY(Tx)=PY(Tx)+PY(Tx)=TxTPY(x)=T[xPY(x)]=TPY(x)(V-10)

推论 V-2:
[ E [ y y T ] ] − 1 [\mathbf{E}[\mathbf{y}\mathbf{y}^{\rm\scriptsize T}]]^{−1} [E[yyT]]1 存在情况下, 任意向量 x \mathbf{x} x 和线性变换 T \mathbf{T} T, 满足 T P Y ⊥ ( x ) ∈ Y ⊥ \mathbf{T} \mathcal{P}_{\mathcal{Y}^{\bot}} (\mathbf{x})\in \mathcal{Y}^{\bot} TPY(x)Y.

证明: 因为 P Y ⊥ ( T x ) ∈ Y ⊥ \mathcal{P}_{\mathcal{Y}^{\bot}} (\mathbf{T}\mathbf{x})\in \mathcal{Y}^{\bot} PY(Tx)Y 及式 (V-10).

VI. Kalman Filter 的正交投影

我们使用归纳法证明 Kalman Filter 的几何意义 —— 基于几何正交投影实现最优状态估计.

1. 时刻 t = 0 t=0 t=0 的初值

x ^ 0 ∣ 0 = E [ x ] (VI-1-1) \hat{\mathbf{x}}_{0 \vert 0} = \mathbf{E}[\mathbf{x}] \tag{VI-1-1} x^0∣0=E[x](VI-1-1)

P 0 ∣ 0 = E [ ( x − x ^ 0 ∣ 0 ) ( x − x ^ 0 ∣ 0 ) T ] (VI-1-2) \mathbf{P}_{0\vert 0} = \mathbf{E}[(x-\hat{\mathbf{x}}_{0 \vert 0})(x-\hat{\mathbf{x}}_{0 \vert 0})^{\rm\scriptsize T}] \tag{VI-1-2} P0∣0=E[(xx^0∣0)(xx^0∣0)T](VI-1-2)

以上两式的成立使得 Kalman Filter 是无偏估计. P 0 ∣ 0 \mathbf{P}_{0\vert 0} P0∣0 是正定对称矩阵.

2. 时刻 t = n − 1 t=n-1 t=n1 的假设

假设 n − 1 n-1 n1 时刻, 给定测量值序列 y 1 , y 2 , … , y n − 1 \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1} y1,y2,,yn1关于状态向量 x n − 1 \mathbf{x}_{n-1} xn1 的最小均方误差估计是 x ^ n − 1 ∣ n − 1 \hat{\mathbf{x}}_{n-1 \vert n-1} x^n1∣n1, 记作 x A \mathbf{x}_{\mathcal{A}} xA. 即
x A = x ^ n − 1 ∣ n − 1 = arg ⁡ min ⁡ x ^ n − 1 ∈ A   E [ ( x n − 1 − x ^ n − 1 ) T   ( x n − 1 − x ^ n − 1 ) ] = arg ⁡ min ⁡ x ^ n − 1 ∈ A   ∥ x n − 1 − x ^ n − 1 ∥ (VI-2-1) \begin{aligned} \mathbf{x}_{\mathcal{A}} &=\hat{\mathbf{x}}_{n-1 \vert n-1} \\ &= {\underset{\hat{\mathbf{x}}_{n-1}\in \mathcal{A}} {\arg \min} } \,\sqrt{\mathbf{E}\left[(\mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1})^{\rm\scriptsize T} \, (\mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1}) \right]}\\ &={\underset{\hat{\mathbf{x}}_{n-1}\in \mathcal{A}} {\arg \min} } \, \Vert \mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1} \Vert \end{aligned} \tag{VI-2-1} xA=x^n1∣n1=x^n1AargminE[(xn1x^n1)T(xn1x^n1)] =x^n1Aargminxn1x^n1(VI-2-1)
估计误差向量为
ϵ n − 1 = x n − 1 − x ^ n − 1 ∣ n − 1 (VI-2-2) {\boldsymbol{\epsilon}}_{n-1} = \mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1\vert n-1} \tag{VI-2-2} ϵn1=xn1x^n1∣n1(VI-2-2)
误差协方差矩阵为
P n − 1 ∣ n − 1 = E [ ( x n − 1 − x ^ n − 1 ∣ n − 1 ) ( x n − 1 − x ^ n − 1 ∣ n − 1 ) T ] = E [ ϵ n − 1 ϵ n − 1 T ] (VI-2-3) \begin{aligned} \mathbf{P}_{n-1 \vert n-1} &= \mathbf{E}[(\mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1\vert n-1}) (\mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1\vert n-1})^{\rm\scriptsize T}] \\ &= \mathbf{E}[{\boldsymbol{\epsilon}}_{n-1} {\boldsymbol{\epsilon}}_{n-1}^{\rm \scriptsize T}] \end{aligned} \tag{VI-2-3} Pn1∣n1=E[(xn1x^n1∣n1)(xn1x^n1∣n1)T]=E[ϵn1ϵn1T](VI-2-3)


Fig.2 上一时刻的正交投影示意图

以 Fig.2 来进行解释, A ⊂ H \mathcal{A} \subset \mathcal{H} AH 是由 y 1 , y 2 , … , y n − 1 \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1} y1,y2,,yn1 张成的子空间. x A \mathbf{x}_{\mathcal{A}} xA x n − 1 \mathbf{x}_{n-1} xn1 沿着 A ⊥ \mathcal{A}^{\bot} A A \mathcal{A} A 正交投影得到的向量, 即
x A = P A ( x n − 1 ) (VI-2-4) \mathbf{x}_{\mathcal{A}} = \mathcal{P}_{\mathcal{A}} (\mathbf{x}_{n-1}) \tag{VI-2-4} xA=PA(xn1)(VI-2-4)
最小状态估计的误差

ϵ n − 1 ∈ A ⊥ (VI-2-5) {\boldsymbol{\epsilon}}_{n-1} \in \mathcal{A}^{\bot} \tag{VI-2-5} ϵn1A(VI-2-5)
由内积的定义式 (III-3) 及式 (V-2-3) 得
x A = arg ⁡ min ⁡ x ^ n − 1 ∈ A   T r ( P n − 1 ∣ n − 1 ) (VI-2-6) \mathbf{x}_{\mathcal{A}} = {\underset{\hat{\mathbf{x}}_{n-1}\in \mathcal{A}} {\arg \min} } \, {\rm{Tr}}( \mathbf{P}_{n-1 \vert n-1}) \tag{VI-2-6} xA=x^n1AargminTr(Pn1∣n1)(VI-2-6)
t = n − 1 t=n-1 t=n1 时刻的代价函数可以写成
J n − 1 = T r ( P n − 1 ∣ n − 1 ) (VI-2-7) \mathbf{J}_{n-1} ={\rm{Tr}}(\mathbf{P}_{n-1 \vert n-1}) \tag{VI-2-7} Jn1=Tr(Pn1∣n1)(VI-2-7)
这里我们假设了 t = n − 1 t=n-1 t=n1 时刻, 关于估计误差的最小二乘法意义下的最优解, 就是状态向量 x n − 1 \mathbf{x}_{n-1} xn1 向由测量向量序列 y 1 , y 2 , … , y n − 1 \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1} y1,y2,,yn1 张成的子空间 A \mathcal{A} A 投影, 得到的投影向量 x A \mathbf{x}_{\mathcal{A}} xA. 换句话说, Kalman Filter 这个最优估计的滤波过程就是计算正交投影子 P A ( ⋅ ) \mathcal{P}_{\mathcal{A}}(\cdot) PA() 的过程.

3. 时刻 t = n t=n t=n 的证明

原有向量 y 1 , y 2 , … , y n − 1 \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1} y1,y2,,yn1 组成子空间 A \mathcal{A} A. 在此基础上新增测量向量 y n \mathbf{y}_n yn, 由 y 1 , y 2 , … , y n − 1 , y n \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1}, \mathbf{y}_{n} y1,y2,,yn1,yn 张成扩展后的子空间 B ⊂ H \mathcal{B} \subset \mathcal{H} BH. 即因为新测量 y n \mathbf{y}_n yn 的加入, B \mathcal{B} B A \mathcal{A} A 的扩展, 写成 A ⊂ B ⊂ H \mathcal{A} \subset \mathcal{B} \subset \mathcal{H} ABH.

Fig.3 所示是 t = n t=n t=n 时刻真实状态向量 x n \mathbf{x}_n xn 与两个子空间 A \mathcal{A} A B \mathcal{B} B 之间的正交投影关系的描述, 梳理这里面的几何关系来解释 Kalman Filter 过程是我们的目标.

3.1. 预测推导

x A \mathbf{x}_{\mathcal{A}} xA 是上一时刻状态向量在子空间 A \mathcal{A} A 中的投影向量, 是上一时刻的最优估计值 (我们这个归纳法的假设条件).

由预测阶段的系统方程 (II-1) 可知
x ^ n ∣ n − 1 = F n x A + B n u n (VI-3-1-1) \hat{\mathbf{x}}_{n \vert n-1} = \mathbf{F}_n {\mathbf{x}}_{\mathcal{A}} + \mathbf{B}_n \mathbf{u}_n \tag{VI-3-1-1} x^nn1=FnxA+Bnun(VI-3-1-1)
t = n t=n t=n 时刻, B n u n \mathbf{B}_n \mathbf{u}_n Bnun 是控制输入中明确的已知的部分, 可能存在的输入扰动通过白噪声 w t \mathbf{w}_t wt 建模 (参考式 (I-1)). 所以式 (VI-3-1-1) 中的 B n u n \mathbf{B}_n \mathbf{u}_n Bnun 在当前时刻看作是常向量, 提供给整个线性系统一个整体的偏移. 这个整体偏移不会对我们的优化和几何描述产生影响, 故为了简化描述, 我们省略这个外部输入的向量项 B n u n \mathbf{B}_n \mathbf{u}_n Bnun.

原始的预测方程 (I-1) 也可以简化为

x t = F t x t − 1 + w t (VI-3-1-2) \mathbf{x}_{t} = \mathbf{F}_{t} \mathbf{x}_{t-1} + \mathbf{w}_{t} \tag{VI-3-1-2} xt=Ftxt1+wt(VI-3-1-2)
式 (VI-3-1-1) 就可以简写为
x A ′ ≜ x ^ n ∣ n − 1 = F n x A (VI-3-1-3) {\mathbf{x}}_{\mathcal{A}}^{'} \triangleq {\hat{\mathbf{x}}}_{n \vert n-1} = \mathbf{F}_n {\mathbf{x}}_{\mathcal{A}} \tag{VI-3-1-3} xAx^nn1=FnxA(VI-3-1-3)

x A ′ {\mathbf{x}}_{\mathcal{A}}^{'} xA 为根据上一时刻的最优状态估计, 利用系统动力学预测的当前时刻的先验估计. 由前面的定理 “Minimum Variance Linear Estimate of a Linear Function” 可知, x A ′ {\mathbf{x}}_{\mathcal{A}}^{'} xA 是系统状态 x n \mathbf{x}_n xn 基于 1 到 t = n − 1 t=n-1 t=n1 时刻测量 y 1 , y 2 , … , y n − 1 \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1} y1,y2,,yn1 得到的最优估计 (预测).

由推论 V-1, x n \mathbf{x}_n xn 正交投影于 A \mathcal{A} A 得到 x A ′ ∈ A {\mathbf{x}}_{\mathcal{A}}^{'} \in \mathcal{A} xAA, 且由定理 “Orthogonal Decomposition” 可知
x n − x A ′ ∈ A ⊥ (VI-3-1-4) \mathbf{x}_n - {\mathbf{x}}_{\mathcal{A}}^{'} \in \mathcal{A}^{\bot} \tag{VI-3-1-4} xnxAA(VI-3-1-4)
对应于 x ^ n ∣ n − 1 {\hat{\mathbf{x}}}_{n \vert n-1} x^nn1 的误差协方差矩阵
P n ∣ n − 1 = E [ ( x n − x ^ n ∣ n − 1 ) ( x n − x ^ n ∣ n − 1 ) T ] = E [ ( F n x n − 1 + w t − F n x ^ n − 1 ∣ n − 1 ) ( F n x n − 1 + w t − F n x ^ n − 1 ∣ n − 1 ) T ] = E [ ( F n ( x n − 1 − x ^ n − 1 ∣ n − 1 ) + w n ) ( F n ( x n − 1 − x ^ n − 1 ∣ n − 1 ) + w n ) T ] = F n P n − 1 ∣ n − 1 F n T + Q n (VI-3-1-5) \begin{aligned} &\mathbf{P}_{n \vert n-1} \\ =& \mathbf{E}[(\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n-1}) (\mathbf{x}_{n} - \hat{\mathbf{x}}_{n \vert n-1})^{\rm\scriptsize T}]\\ =& \mathbf{E}[(\mathbf{F}_{n}\mathbf{x}_{n-1} + \mathbf{w}_t- \mathbf{F}_n \hat{\mathbf{x}}_{n-1\vert n-1} ) (\mathbf{F}_{n}\mathbf{x}_{n-1} + \mathbf{w}_t - \mathbf{F}_n \hat{\mathbf{x}}_{n-1\vert n-1} )^{\rm\scriptsize T}]\\ =& \mathbf{E}[(\mathbf{F}_{n}(\mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1\vert n-1}) + \mathbf{w}_n) (\mathbf{F}_{n}(\mathbf{x}_{n-1} - \hat{\mathbf{x}}_{n-1\vert n-1}) + \mathbf{w}_n)^{\rm\scriptsize T}]\\ =& \mathbf{F}_{n} \mathbf{P}_{n-1 \vert n-1} \mathbf{F}_{n}^{\rm\scriptsize T} + \mathbf{Q}_n \end{aligned} \tag{VI-3-1-5} ====Pnn1E[(xnx^nn1)(xnx^nn1)T]E[(Fnxn1+wtFnx^n1∣n1)(Fnxn1+wtFnx^n1∣n1)T]E[(Fn(xn1x^n1∣n1)+wn)(Fn(xn1x^n1∣n1)+wn)T]FnPn1∣n1FnT+Qn(VI-3-1-5)

P n ∣ n − 1 \mathbf{P}_{n \vert n-1} Pnn1 是正定对称矩阵.


Fig.3 当前时刻的正交投影示意图

有了状态向量的预测 x A ′ {\mathbf{x}}_{\mathcal{A}}^{'} xA 之后, 可以进一步根据测量方程 (I-3) 预测该状态向量条件下的测量向量, 即 H n x A ′ \mathbf{H}_{n} {\mathbf{x}}_{\mathcal{A}}^{'} HnxA.

根据 “Minimum Variance Linear Estimate of a Linear Function”, H n x A ′ \mathbf{H}_{n} {\mathbf{x}}_{\mathcal{A}}^{'} HnxA 是测量数据 y n \mathbf{y}_n yn 基于 1 时刻到 t = n − 1 t=n-1 t=n1 时刻测量 y 1 , y 2 , … , y n − 1 \mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_{n-1} y1,y2,,yn1 得到的最优估计 (预测), 即由 y n \mathbf{y}_n yn 正交投影于 A \mathcal{A} A 得到. 定义符号
y A ≜ H n x A ′ (VI-3-1-6) \mathbf{y}_{\mathcal{A}} \triangleq \mathbf{H}_{n} {\mathbf{x}}_{\mathcal{A}}^{'} \tag{VI-3-1-6} yAHnxA(VI-3-1-6)

将已测量到的 y n \mathbf{y}_n yn 与预测得到 y A \mathbf{y}_\mathcal{A} yA 之间的误差定义为新息 (Innovation variable)
y ~ n ≜ y n − y A (VI-3-1-7) \tilde{\mathbf{y}}_n \triangleq {\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}} \tag{VI-3-1-7} y~nynyA(VI-3-1-7)

对应于 y A \mathbf{y}_\mathcal{A} yA 的误差协方差矩阵
E [ ( y n − y A ) ( y n − y A ) T ] = E [ ( H n x n + v n − H n x ^ n ∣ n − 1 ) ( H n x n + v n − H n x ^ n ∣ n − 1 ) T ] = E [ ( H n ( x n − x ^ n ∣ n − 1 ) + v n ) ( H n ( x n − x ^ n ∣ n − 1 ) + v n ) T ] = H n P n ∣ n − 1 H n T + R n (VI-3-1-8) \begin{aligned} &\mathbf{E}[(\mathbf{y}_{n} - \mathbf{y}_{\mathcal{A}}) (\mathbf{y}_{n} - \mathbf{y}_{\mathcal{A}})^{\rm\scriptsize T}]\\ =& \mathbf{E}[( \mathbf{H}_n \mathbf{x}_n +\mathbf{v}_n -\mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1}) (\mathbf{H}_n \mathbf{x}_n +\mathbf{v}_n -\mathbf{H}_{n}{\hat{\mathbf{x}}}_{n \vert n-1})^{\rm\scriptsize T}] \\ =& \mathbf{E}[(\mathbf{H}_{n}(\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n-1}) + \mathbf{v}_n) (\mathbf{H}_{n}(\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n-1}) + \mathbf{v}_n)^{\rm\scriptsize T}]\\ =& \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_n \end{aligned} \tag{VI-3-1-8} ===E[(ynyA)(ynyA)T]E[(Hnxn+vnHnx^nn1)(Hnxn+vnHnx^nn1)T]E[(Hn(xnx^nn1)+vn)(Hn(xnx^nn1)+vn)T]HnPnn1HnT+Rn(VI-3-1-8)
这里定义符号
P y ≜ H n P n ∣ n − 1 H n T + R n (VI-3-1-9) \mathbf{P}_y \triangleq \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_n \tag{VI-3-1-9} PyHnPnn1HnT+Rn(VI-3-1-9)
P y \mathbf{P}_y Py 就是新息协方差矩阵, 也是正定对称矩阵 (保证逆矩阵存在).

3.2. 更新推导

我们分析一下新息, 由定理 “Orthogonal Decomposition” 可知
y ~ n = y n − y A ∈ A ⊥ (VI-3-2-1) \tilde{\mathbf{y}}_n = {\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}} \in \mathcal{A}^{\bot} \tag{VI-3-2-1} y~n=ynyAA(VI-3-2-1)
实时上 y A \mathbf{y}_\mathcal{A} yA 是有 A \mathcal{A} A 上原有信息构成的, 并不提供新的信息, 而 y ~ n \tilde{\mathbf{y}}_n y~n 是无法由 A \mathcal{A} A 上原有信息构成的新的信息, 故而命名. 类似于 Gram-Schmidt正交化, 新的子空间 B \mathcal{B} B 是由原有子空间 A \mathcal{A} A 扩展加上新的向量基 y ~ n \tilde{\mathbf{y}}_n y~n 得到的. 换一种说法是 B = A ⊕ A ⊥ \mathcal{B} = \mathcal{A} \oplus A^{\bot} B=AA.

进一步, 对 t = n t=n t=n 时刻的状态向量 x n \mathbf{x}_n xn B \mathcal{B} B B ⊥ \mathcal{B}^{\bot} B 两个正交子空间内进行直和分解, 得到在子空间 B \mathcal{B} B 的投影 (最优估计) 为 x B \mathbf{x}_\mathcal{B} xB. 参考 [6] 构造新的变量
z n = x n − x A ′ + K ‾ n ( y n − y A ) (VI-3-2-2) \mathbf{z}_n = \mathbf{x}_{\mathcal{n}} - \mathbf{x}_{\mathcal{A}}^{'} + \overline{\mathbf{K}}_n({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) \tag{VI-3-2-2} zn=xnxA+Kn(ynyA)(VI-3-2-2)
下面试图选取合适的 K ‾ n \overline{\mathbf{K}}_n Kn 使得 z n ⊥ B \mathbf{z}_n \bot \mathcal{B} znB.

由式 (VI-3-1-4) 、式 (VI-3-2-1) 以及推论 V-2 可知, 任意的 K ‾ n \overline{\mathbf{K}}_n Kn 都使得
z n   ⊥   A (VI-3-2-3) \mathbf{z}_n\, \bot\, \mathcal{A} \tag{VI-3-2-3} znA(VI-3-2-3)
故要使 z n ⊥ B \mathbf{z}_n \bot \mathcal{B} znB, 只要再让 z n \mathbf{z}_n zn 垂直于新扩展的向量基 y ~ n \tilde{\mathbf{y}}_n y~n, 即
⟨ z n , y n − y A ⟩ = 0 (VI-3-2-4) \langle \mathbf{z}_n, \mathbf{y}_n - \mathbf{y}_{\mathcal{A}}\rangle = 0 \tag{VI-3-2-4} zn,ynyA=0(VI-3-2-4)
上式左边计算展开, 得到
⟨ z n , y n − y A ⟩ = ⟨ x n − x A ′ + K ‾ n ( y n − y A ) , y n − y A ⟩ = ⟨ x n − x A ′ , y n − y A ⟩ + ⟨ K ‾ n ( y n − y A ) , y n − y A ⟩ = T r { E [ ( x n − x ^ n ∣ n − 1 ) ( H n ( x n − x ^ n ∣ n − 1 ) + v t ) T ] } + T r { K ‾ n E [ ( y n − y A ) ( y n − y A ) T ] } = T r [ P n ∣ n − 1 H n T + K ‾ n ( H n P n ∣ n − 1 H n T + R n ) ] (VI-3-2-5) \begin{aligned} &\langle \mathbf{z}_n, \mathbf{y}_n - \mathbf{y}_{\mathcal{A}}\rangle \\ =& \langle \mathbf{x}_{\mathcal{n}} - \mathbf{x}_{\mathcal{A}}^{'} + \overline{\mathbf{K}}_n ({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) , \mathbf{y}_n - \mathbf{y}_{\mathcal{A}}\rangle\\ =& \langle \mathbf{x}_{\mathcal{n}} - \mathbf{x}_{\mathcal{A}}^{'} , \mathbf{y}_n - \mathbf{y}_{\mathcal{A}}\rangle + \langle \overline{\mathbf{K}}_n ({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) , \mathbf{y}_n - \mathbf{y}_{\mathcal{A}}\rangle \\ =&{\rm Tr} \{\mathbf{E}[( \mathbf{x}_{\mathcal{n}} -{\hat{\mathbf{x}}}_{n \vert n-1}) (\mathbf{H}_{n}(\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n-1}) + \mathbf{v}_t)^{\rm\scriptsize T} ] \}+{\rm Tr}\{ \overline{\mathbf{K}}_{n} \mathbf{E}[(\mathbf{y}_{n} - \mathbf{y}_{\mathcal{A}}) (\mathbf{y}_{n} - \mathbf{y}_{\mathcal{A}})^{\rm\scriptsize T}]\} \\ =& {\rm Tr}[ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} +\overline{\mathbf{K}}_{n}( \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_n )] \end{aligned} \tag{VI-3-2-5} ====zn,ynyAxnxA+Kn(ynyA),ynyAxnxA,ynyA+Kn(ynyA),ynyATr{E[(xnx^nn1)(Hn(xnx^nn1)+vt)T]}+Tr{KnE[(ynyA)(ynyA)T]}Tr[Pnn1HnT+Kn(HnPnn1HnT+Rn)](VI-3-2-5)
一个矩阵的迹为零的一个充分条件是该矩阵为零矩阵, 故式 (VI-3-2-4) 成立的充分条件
K ‾ n = − P n ∣ n − 1 H n T   ( H n P n ∣ n − 1 H n T + R n ) − 1 (VI-3-2-6) \overline{\mathbf{K}}_{n} = - \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T}\, ( \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_n ) ^{-1} \tag{VI-3-2-6} Kn=Pnn1HnT(HnPnn1HnT+Rn)1(VI-3-2-6)
与标准 Kalman Filter 中式 (II-5) 比较, 此处 K ‾ n \overline{\mathbf{K}}_{n} Kn 与式 (II-5) 中的 K t \mathbf{K}_{t} Kt 之间相差负号, 即
K n = − K ‾ n = P n ∣ n − 1 H n T   ( H n P n ∣ n − 1 H n T + R n ) − 1 = P n ∣ n − 1 H n T P y − 1 (VI-3-2-7) \begin{aligned} {\mathbf{K}}_{n} &= -\overline{\mathbf{K}}_{n} \\ &= \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T}\, ( \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_n ) ^{-1} \\ & = \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} \end{aligned} \tag{VI-3-2-7} Kn=Kn=Pnn1HnT(HnPnn1HnT+Rn)1=Pnn1HnTPy1(VI-3-2-7)
已知 K ‾ n \overline{\mathbf{K}}_{n} Kn 取式 (VI-3-2-6) 时, z n ⊥ B \mathbf{z}_n \bot \mathcal{B} znB. 变换等式 (VI-3-2-2) 得到
x n = x A ′ + P n ∣ n − 1 H n T   ( H n P n ∣ n − 1 H n T + R t ) − 1 ( y n − y A ) + z n (VI-3-2-8) \begin{aligned} \mathbf{x}_{\mathcal{n}} = \mathbf{x}_{\mathcal{A}}^{'} + \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T}\, ( \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_t ) ^{-1} ({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) +\mathbf{z}_n \end{aligned}\tag{VI-3-2-8} xn=xA+Pnn1HnT(HnPnn1HnT+Rt)1(ynyA)+zn(VI-3-2-8)
对式 (VI-3-2-8) 右手边每一项分析他们的所处几何状态:

- 第三项 z n ∈ B ⊥ \mathbf{z}_n \in \mathcal{B}^{\bot} znB

- 第一项 x A ′ ∈ A ⊂ B {\mathbf{x}}_{\mathcal{A}}^{'} \in \mathcal{A} \subset \mathcal{B} xAAB

- 已知 y n − y A ∈ A ⊥ ⊂ B {\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}} \in \mathcal{A}^{\bot} \subset \mathcal{B} ynyAAB, 由推论 V-2 可知, 第二项 K n ( y n − y A ) ∈ A ⊥ ⊂ B {\mathbf{K}}_{n} ({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) \in \mathcal{A}^{\bot} \subset \mathcal{B} Kn(ynyA)AB.

由式 (VI-3-2-8) 各项几何状态, 我们定义 x n \mathbf{x}_n xn 在子空间 B \mathcal{B} B 上的正交投影所得向量 x B \mathbf{x}_{\mathcal{B}} xB. 本身属于 B \mathcal{B} B 的向量在正交投影 P B ( ⋅ ) \mathcal{P}_{\mathcal{B}}(\cdot) PB() 作用下不变, 而属于 B ⊥ \mathcal{B}^{\bot} B 的向量在正交投影 P B ( ⋅ ) \mathcal{P}_{\mathcal{B}}(\cdot) PB() 作用下消失, 即
x B = P B ( x n ) = x A ′ + K n ( y n − y A ) = x A ′ + P n ∣ n − 1 H n T   ( H n P n ∣ n − 1 H n T + R t ) − 1 ( y n − y A ) (VI-3-2-9) \begin{aligned} \mathbf{x}_{\mathcal{B}} &= \mathcal{P}_{\mathcal{B}}(\mathbf{x}_n)\\ &= \mathbf{x}_{\mathcal{A}}^{'} + {\mathbf{K}}_{n} ({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) \\ &= \mathbf{x}_{\mathcal{A}}^{'} + \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T}\, ( \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} + \mathbf{R}_t ) ^{-1} ({\mathbf{y}}_n - \mathbf{y}_{\mathcal{A}}) \end{aligned}\tag{VI-3-2-9} xB=PB(xn)=xA+Kn(ynyA)=xA+Pnn1HnT(HnPnn1HnT+Rt)1(ynyA)(VI-3-2-9)
x B \mathbf{x}_{\mathcal{B}} xB 即为基于上一时刻的最优估计 x A \mathbf{x}_{\mathcal{A}} xA 及当前时刻的测量 y n \mathbf{y}_n yn 得到的当前时刻上的最优估计 (当前时刻 t = n t=n t=n). x B \mathbf{x}_{\mathcal{B}} xB 也记作 x ^ n ∣ n \hat{\mathbf{x}}_{n \vert n} x^nn, 式 (VI-3-2-9) 可以写成
x ^ n ∣ n = x ^ n ∣ n − 1 + K n ( y n − H n x ^ n ∣ n − 1 ) (VI-3-2-10) \hat{\mathbf{x}}_{n \vert n} = {\hat{\mathbf{x}}}_{n \vert n-1} + \mathbf{K}_{n} ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) \tag{VI-3-2-10} x^nn=x^nn1+Kn(ynHnx^nn1)(VI-3-2-10)

几何推导得到的式 (VI-3-2-10) 与标准 Kalman Filter 更新方程 (II-5) 完全一致.

做到这里基本上已经实现基于正交投影的几何推导与标准 Kalman Filter 之间殊途同归了, 下面再验算一下更新阶段的误差协方差矩阵.

3.3. 误差协方差

n = t n=t n=t 时刻, 估计值与真实值之间的误差为 ϵ n = x n − x ^ n ∣ n {\boldsymbol{\epsilon}}_{n} =\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n} ϵn=xnx^nn. 误差协方差矩阵
P n ∣ n = E [ ( x n − x ^ n ∣ n ) ( x n − x ^ n ∣ n ) T ] = E [ { ( x n − x ^ n ∣ n − 1 ) − K n ( y n − H n x ^ n ∣ n − 1 ) }   { ( x n − x ^ n ∣ n − 1 ) − K n ( y n − H n x ^ n ∣ n − 1 ) } T ] = E [ ( x n − x ^ n ∣ n − 1 ) ( x n − x ^ n ∣ n − 1 ) T ] − E [ ( x n − x ^ n ∣ n − 1 ) ( y n − H n x ^ n ∣ n − 1 ) T ] K n T − K n E [ ( y n − H n x ^ n ∣ n − 1 ) ( x n − x ^ n ∣ n − 1 ) T ] + K n E [ ( y n − H n x ^ n ∣ n − 1 ) ( y n − H n x ^ n ∣ n − 1 ) T ] K n T (VI-3-3-1) \begin{aligned} \mathbf{P}_{n \vert n} &= \mathbf{E}[(\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n}) (\mathbf{x}_{n} - \hat{\mathbf{x}}_{n\vert n})^{\rm\scriptsize T}] \\ &=\mathbf{E}\left[ \{(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1}) - \mathbf{K}_{n} ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} )\}\, \{(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1}) - \mathbf{K}_{n} ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) \}^{\rm\scriptsize T}\right] \\ &=\mathbf{E}[(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1})(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1})^{\rm\scriptsize T}] \\ & \quad - \mathbf{E}[(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1}) ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) ^{\rm\scriptsize T}]\mathbf{K}_{n}^{\rm\scriptsize T} \\ & \quad - \mathbf{K}_{n} \mathbf{E}[ ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) (\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1}) ^{\rm\scriptsize T}] \\ & \quad + \mathbf{K}_{n} \mathbf{E}[ ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} )^{\rm\scriptsize T}] \mathbf{K}_{n}^{\rm\scriptsize T}\\ \end{aligned} \tag{VI-3-3-1} Pnn=E[(xnx^nn)(xnx^nn)T]=E[{(xnx^nn1)Kn(ynHnx^nn1)}{(xnx^nn1)Kn(ynHnx^nn1)}T]=E[(xnx^nn1)(xnx^nn1)T]E[(xnx^nn1)(ynHnx^nn1)T]KnTKnE[(ynHnx^nn1)(xnx^nn1)T]+KnE[(ynHnx^nn1)(ynHnx^nn1)T]KnT(VI-3-3-1)
我们分开计算,

- 根据式 (VI-3-1-5)可知, 式 (VI-3-2-10) 中的第一项为 P n ∣ n − 1 \mathbf{P}_{n \vert n-1} Pnn1.

- 式 (VI-3-2-10) 中的第二项为
− E [ ( x n − x ^ n ∣ n − 1 ) ( y n − H n x ^ n ∣ n − 1 ) T ] K n T = − E [ ( x n − x ^ n ∣ n − 1 ) ( H n x n + v n − H n x ^ n ∣ n − 1 ) T ] K n T = − P n ∣ n − 1 H n T K n T \begin{aligned} &-\mathbf{E}[(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1}) ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) ^{\rm\scriptsize T}]\mathbf{K}_{n}^{\rm\scriptsize T} \\ =& -\mathbf{E}[(\mathbf{x}_{n} - {\hat{\mathbf{x}}}_{n \vert n-1}) (\mathbf{H}_n \mathbf{x}_n +\mathbf{v}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) ^{\rm\scriptsize T}]\mathbf{K}_{n}^{\rm\scriptsize T} \\ =& -\mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \mathbf{K}_{n}^{\rm\scriptsize T} \end{aligned} ==E[(xnx^nn1)(ynHnx^nn1)T]KnTE[(xnx^nn1)(Hnxn+vnHnx^nn1)T]KnTPnn1HnTKnT
- 式 (VI-3-2-10) 中的第三项为 − K n H n P n ∣ n − 1 - \mathbf{K}_{n} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} KnHnPnn1.

- 式 (VI-3-2-10) 中的第四项为
K n   E [ ( y n − H n x ^ n ∣ n − 1 ) ( y n − H n x ^ n ∣ n − 1 ) T ]   K n T = K n   E [ ( H n x n + v n − H n x ^ n ∣ n − 1 ) ( H n x n + v n − H n x ^ n ∣ n − 1 ) T ]   K n T = K n H n P n ∣ n − 1 H n T K n T + K n R n K n T \begin{aligned} &\mathbf{K}_{n}\, \mathbf{E}[ ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) ({\mathbf{y}}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} )^{\rm\scriptsize T}]\, \mathbf{K}_{n}^{\rm\scriptsize T} \\ =& \mathbf{K}_{n} \, \mathbf{E}[ (\mathbf{H}_n \mathbf{x}_n +\mathbf{v}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} ) (\mathbf{H}_n \mathbf{x}_n +\mathbf{v}_n - \mathbf{H}_{n} {\hat{\mathbf{x}}}_{n \vert n-1} )^{\rm\scriptsize T}] \,\mathbf{K}_{n}^{\rm\scriptsize T}\\ =& \mathbf{K}_{n} \mathbf{H}_n \mathbf{P}_{n \vert n-1} \mathbf{H}_n^{\rm\scriptsize T} \mathbf{K}_{n}^{\rm\scriptsize T} + \mathbf{K}_{n} \mathbf{R}_{n} \mathbf{K}_{n}^{\rm\scriptsize T} \end{aligned} ==KnE[(ynHnx^nn1)(ynHnx^nn1)T]KnTKnE[(Hnxn+vnHnx^nn1)(Hnxn+vnHnx^nn1)T]KnTKnHnPnn1HnTKnT+KnRnKnT
这四项分别代入式 (VI-3-2-10) 得
P n ∣ n = P n ∣ n − 1 − P n ∣ n − 1 H n T K n T − K n H n P n ∣ n − 1 + K n H n P n ∣ n − 1 H n T K n T + K n R n K n T = ( I − K n H n ) P n ∣ n − 1 ( I − K n H n ) T + K n R n K n T (VI-3-3-2) \begin{aligned} \mathbf{P}_{n \vert n} =& \mathbf{P}_{n \vert n-1} - \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \mathbf{K}_{n}^{\rm\scriptsize T} - \mathbf{K}_{n} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} + \mathbf{K}_{n} \mathbf{H}_n \mathbf{P}_{n \vert n-1} \mathbf{H}_n^{\rm\scriptsize T} \mathbf{K}_{n}^{\rm\scriptsize T} + \mathbf{K}_{n} \mathbf{R}_{n} \mathbf{K}_{n}^{\rm\scriptsize T} \\ =& (\mathbf{I} - \mathbf{K}_n \mathbf{H}_n) \mathbf{P}_{n \vert n-1}(\mathbf{I} - \mathbf{K}_n \mathbf{H}_n)^{\rm\scriptsize T}+ \mathbf{K}_{n} \mathbf{R}_{n} \mathbf{K}_{n}^{\rm\scriptsize T} \end{aligned} \tag{VI-3-3-2} Pnn==Pnn1Pnn1HnTKnTKnHnPnn1+KnHnPnn1HnTKnT+KnRnKnT(IKnHn)Pnn1(IKnHn)T+KnRnKnT(VI-3-3-2)
这是更新阶段误差协方差矩阵的一种形式.

由式 (VI-3-2-7) K n = P n ∣ n − 1 H n T P y − 1 \mathbf{K}_{n} = \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} Kn=Pnn1HnTPy1 . 又有相关协方差矩阵的对称性, 可知 K n T = P y − 1 H n P n ∣ n − 1 \mathbf{K}_{n}^{\rm\scriptsize T} = \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} KnT=Py1HnPnn1.

代入式 (VI-3-3-2) 得到
P n ∣ n = P n ∣ n − 1 − P n ∣ n − 1 H n T P y − 1 H n P n ∣ n − 1 ‾ − P n ∣ n − 1 H n T P y − 1 ‾ H n P n ∣ n − 1 + P n ∣ n − 1 H n T P y − 1 ‾ H n P n ∣ n − 1 H n T P y − 1 H n P n ∣ n − 1 ‾ + P n ∣ n − 1 H n T P y − 1 ‾ R n P y − 1 H n P n ∣ n − 1 ‾ = P n ∣ n − 1 − 2 P n ∣ n − 1 H n T P y − 1 H n P n ∣ n − 1 + P n ∣ n − 1 H n T P y − 1 ( H n P n ∣ n − 1 H n T + R n ) ‾ P y − 1 H n P n ∣ n − 1 = P n ∣ n − 1 + 2 P n ∣ n − 1 H n T P y − 1 H n P n ∣ n − 1 + P n ∣ n − 1 H n T P y − 1 P y ‾ P y − 1 H n P n ∣ n − 1 = P n ∣ n − 1 − 2 P n ∣ n − 1 H n T P y − 1 H n P n ∣ n − 1 + P n ∣ n − 1 H n T P y − 1 H n P n ∣ n − 1 = P n ∣ n − 1 − P n ∣ n − 1 H n T P y − 1 ‾ H n P n ∣ n − 1 = P n ∣ n − 1 − K n H n P n ∣ n − 1 (VI-3-3-3) \begin{aligned} \mathbf{P}_{n \vert n} =& \mathbf{P}_{n \vert n-1} \\ &- \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \underline{\mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} }\\ &- \underline{ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} }\mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \\ &+ \underline{ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} } \mathbf{H}_n \mathbf{P}_{n \vert n-1} \mathbf{H}_n^{\rm\scriptsize T} \underline{ \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} }\\ &+ \underline{ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} } \mathbf{R}_{n} \underline{ \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1}} \\ =& \mathbf{P}_{n \vert n-1} \\ &- 2 \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \\ &+ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} \underline{\left(\mathbf{H}_n \mathbf{P}_{n \vert n-1} \mathbf{H}_n^{\rm\scriptsize T} + \mathbf{R}_{n}\right)} \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1}\\ =& \mathbf{P}_{n \vert n-1} \\ &+ 2 \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \\ &+ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} \underline{\mathbf{P}_y} \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1}\\ =& \mathbf{P}_{n \vert n-1} \\ &- 2 \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \\ &+ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm \scriptsize T} \mathbf{P}_y^{-1} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1}\\ =& \mathbf{P}_{n \vert n-1} - \underline{ \mathbf{P}_{n \vert n-1} \mathbf{H}_{n}^{\rm\scriptsize T} \mathbf{P}_y^{-1}} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \\ = & \mathbf{P}_{n \vert n-1} - {\mathbf{K}_n} \mathbf{H}_{n} \mathbf{P}_{n \vert n-1} \end{aligned} \tag{VI-3-3-3} Pnn======Pnn1Pnn1HnTPy1HnPnn1Pnn1HnTPy1HnPnn1+Pnn1HnTPy1HnPnn1HnTPy1HnPnn1+Pnn1HnTPy1RnPy1HnPnn1Pnn12Pnn1HnTPy1HnPnn1+Pnn1HnTPy1(HnPnn1HnT+Rn)Py1HnPnn1Pnn1+2Pnn1HnTPy1HnPnn1+Pnn1HnTPy1PyPy1HnPnn1Pnn12Pnn1HnTPy1HnPnn1+Pnn1HnTPy1HnPnn1Pnn1Pnn1HnTPy1HnPnn1Pnn1KnHnPnn1(VI-3-3-3)
与标准 Kalman Filter 中更新阶段误差协方差矩阵 (II-4) 一致.

3.4. 卡尔曼增益

t = n t=n t=n 时刻的代价函数为
J n = T r ( P n ∣ n ) = T r [ ( I − K n H n ) P n ∣ n − 1 ( I − K n H n ) T + K n R n K n T ] (VI-3-4-1) \begin{aligned} \mathbf{J}_{n} &={\rm{Tr}}(\mathbf{P}_{n \vert n})\\ &={\rm{Tr}}\left[ (\mathbf{I} - \mathbf{K}_n \mathbf{H}_n) \mathbf{P}_{n \vert n-1}(\mathbf{I} - \mathbf{K}_n \mathbf{H}_n)^{\rm\scriptsize T}+ \mathbf{K}_{n} \mathbf{R}_{n} \mathbf{K}_{n}^{\rm\scriptsize T}\right] \end{aligned} \tag{VI-3-4-1} Jn=Tr(Pnn)=Tr[(IKnHn)Pnn1(IKnHn)T+KnRnKnT](VI-3-4-1)
直接从代价函数也可以求出卡尔曼增益 K n \mathbf{K}_n Kn,
∂ J n ∂ K n = − 2 ( I − K n H n ) P n ∣ n − 1 H n T + 2 K n R n = 0 (VI-3-4-2) \frac{\partial{\mathbf{J}_n}}{\partial{\mathbf{K}_n}} = -2 (\mathbf{I} - \mathbf{K}_n \mathbf{H}_n) \mathbf{P}_{n \vert n-1} \mathbf{H}_n^{\rm\scriptsize T} + 2 \mathbf{K}_{n} \mathbf{R}_{n} = \mathbf{0} \tag{VI-3-4-2} KnJn=2(IKnHn)Pnn1HnT+2KnRn=0(VI-3-4-2)
一样可以计算得到式 (VI-3-2-7) 中的 K n \mathbf{K}_n Kn.

这样我们就可以对卡尔曼滤波的关键项 —— 卡尔曼增益 Kalman gain 进行如下解读:

- 以信息融合的观点, 卡尔曼滤波增益 K n \mathbf{K}_n Kn 是对预测与测量融合过程中不确定性 (置信程度) 的反映.

- 以最优化的角度, 卡尔曼滤波增益 K n \mathbf{K}_n Kn 是最小二乘法意义下的最优解, 可以获得最小均方误差估计.

- 几何直观的解释, 卡尔曼滤波增益 K n \mathbf{K}_n Kn 是为了获得状态向量 x n \mathbf{x}_n xn 在扩展子空间 B \mathcal{B} B 上的正交投影 x B \mathbf{x}_{\mathcal{B}} xB.

到此完成归纳法第三步, 当 t = n t=n t=n 时刻时, 利用几何正交投影获得的最优状态估计和标准 Kalman Filter 完全一致.
这就完成了对 Kalman Filter 的几何直观解读.

VII 总结

以上我们利用了 Hilbert 空间中正交投影的几何概念, 借助归纳法证明 (解读) 了卡尔曼滤波器最优估计的几何意义:

- 预测和更新均为正交投影意义下的最优估计;

- 新的测量是对原有子空间的扩展;

- 每一步状态估计都是在新扩展空间上的正交投影;

- 新息 (Innovation Variable) 在历史测量张成子空间的正交空间上, 可看做新的扩展维度上的向量基;

- 卡尔曼滤波增益 (Kalman Gain) 是为了获得状态向量在扩展子空间上的正交投影.

本文结束, 欢迎指正交流.

参考文献

[1] Ramsey Faragher, “Understanding the Basis of the Kalman Filter Via a Simple and Intuitive Derivation [Lecture Notes]”, IEEE Signal Processing Magazine, 2012

[2] Wolfram MathWorld, “Hilbert Space”, https://mathworld.wolfram.com/HilbertSpace.html

[3] Onno van Gaans, “Probability measures on metric spaces”, “https://www.math.leidenuniv.nl/~vangaans/jancol1.pdf”

[4] Yang Chen, Shaoshu Li, “A Brief Introduction to Hilbert Space”, https://sites.math.washington.edu/~farbod/teaching/cornell/math6210pdf/math6210Hilbert.pdf, 2016

[5] David G. Luenberger, “Optimization by Vector Space Methods”, Wiley-Interscience, 1969

[6] https://dornsife.usc.edu/assets/sites/1185/docs/507b/Linear.pdf

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
Kalman滤波是一种常用的状态估计方法,在嵌入式系统、机器人、自动控制等领域得到广泛应用。下面是一个基于Matlab实现的Kalman滤波状态估计代码示例。 ``` % 系统参数 Ts = 0.1; % 采样时间 A = [1 Ts; 0 1]; % 系统矩阵 B = [Ts^2/2; Ts]; % 输入矩阵 C = [1 0]; % 观测矩阵 Q = [0.1 0; 0 Ts^2/3]; % 系统噪声协方差 R = 1; % 观测噪声协方差 % 状态初始化 x0 = [0; 0]; % 初始状态 P0 = [1 0; 0 1]; % 初始协方差 % 生成随机信号 t = 0:Ts:10; u = 0.1*randn(1,length(t)); % 输入信号 z = C*x0 + sqrt(R)*randn(1); % 初始观测信号 x_true = zeros(2,length(t)); % 真实状态 x_est = zeros(2,length(t)); % 估计状态 x_true(:,1) = x0; x_est(:,1) = x0; % Kalman滤波 for k = 2:length(t) % 真实状态 x_true(:,k) = A*x_true(:,k-1) + B*u(k-1) + sqrt(Q)*randn(2,1); % 观测信号 z(k) = C*x_true(:,k) + sqrt(R)*randn(1); % 状态估计 x_pred = A*x_est(:,k-1) + B*u(k-1); P_pred = A*P0*A' + Q; K = P_pred*C'/(C*P_pred*C' + R); x_est(:,k) = x_pred + K*(z(k) - C*x_pred); P0 = (eye(2) - K*C)*P_pred; end % 绘图 figure; subplot(211); plot(t,x_true(1,:),'r',t,x_est(1,:),'b'); legend('真实状态','估计状态'); subplot(212); plot(t,x_true(2,:),'r',t,x_est(2,:),'b'); legend('真实状态','估计状态'); ``` 其中,`Ts`为采样时间,`A`、`B`、`C`分别为系统矩阵、输入矩阵和观测矩阵,`Q`为系统噪声协方差,`R`为观测噪声协方差,`x0`为初始状态,`P0`为初始协方差,`u`为输入信号,`z`为观测信号,`x_true`为真实状态,`x_est`为估计状态。在循环中,先计算真实状态和观测信号,然后进行状态估计,更新状态和协方差。最后绘制真实状态估计状态的图像。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值