对 MMSE 估计器与 Kalman 滤波器的理解

对 MMSE 估计器与 Kalman 滤波器的理解

通常说 Kalman Filter (KF) 是线性动力学模型下的最优估计器，这一说法应当如何解释？或者更基本的，KF 的理论基础是什么？它和常用的低通滤波器 (LPF) 有什么本质上的差异？本文尝试对这些问题进行简单的说明。

KF 的特点

KF 与常用低通滤波器的差异

一般信号处理常用的滤波器 (如 FIR, IIR, Butterworth, Chebyshev) 基于拉普拉斯或傅立叶变换，它们都属于线性时不变 (LTI) 滤波器，滤波器的特性 (如频宽) 与输入信号是独立的，也就是说不随信号而变化。要计算输出信号，先要推出或定义滤波器的传递函数或脉冲响应，再与输入信号进行频域相乘或时域卷积得到。

与之对应，KF 是由随机信号的估计出发，输入信号与噪声的统计信息决定了 KF 的特性，另外，KF 基本上是时变的滤波器，这两点是 KF 与常用的 LPF 的主要区别。

KF 与常用的 LPF 的类比

KF 可以视为一个时变 Kalman 增益 ( K K) 的反馈控制器，因此可以看做一个时变的 IIR 滤波器，频宽随时间改变。

如果 KF 的稳态存在，那么在稳态时， K K 基本为常数，近似为一个时不变滤波器。但是， K K 由输入信号和噪声决定，和一般的 LPF 仍然有区别。

MMSE 估计器

随机变量 MMSE 估计器

先从最简单的单变量开始。随机变量  X X 的 MMSE 估计器是找到  x^ x^ 令  E(|X−x^|2) E(|X−x^|2) 最小。很容易证明： x^=E(X) x^=E(X)， MSE=E(|X−E(X)|2)=σ2 MSE=E(|X−E(X)|2)=σ2。

当有两个变量  X X,  Z Z 时，MMSE 估计器变为： z z 可观测，找到不可观测的  x x 的估计  x^ x^ 使得  E(|X−x^|2) E(|X−x^|2) 最小，其中： E(|X−x^|2)=Ex[Ez(|X−x^|2|z)] E(|X−x^|2)=Ex[Ez(|X−x^|2|z)]。

很明显如果  X X， Z Z 相互独立或互不相关， z z 是否可观测对  X X 的估计没有帮助：

x^=E(X)andMSE=σ2x x^=E(X)andMSE=σx2

如果  X X 与  Z Z 相关， z z 可以对  X X 的估计提供帮助，可以证明： x^=E(X|Z=z) x^=E(X|Z=z) 是最小 MSE 估计器，并且 MSE 比没有  z z 的情况小。直观上来讲可以解释为：由于  X X 与  Z Z 存在关联， z z 的信息对减小  X X 的估计误差可以提供帮助。注意表达式  x^=E(X|Z=z) x^=E(X|Z=z) 并没有假设任何  X X 和  Z Z 的统计特性（比如 Gaussian），这是 Bayesian 估计器的基础。

下面再假设  X X,  Z Z 均为 Gaussian 的情况来进一步说明。 X X 和  Z Z 的联合分布为

[x,z]∼N([E(X),E(Z)],[σ2x,σ(xz);σ(xz),σ2z]) [x,z]∼N([E(X),E(Z)],[σx2,σ(xz);σ(xz),σz2])

X X的估计和MSE分别为

x^=E(X|Z=z)=E(X)+σxz(z−E(Z))/σ2z x^=E(X|Z=z)=E(X)+σxz(z−E(Z))/σz2

MSE=σ2x–σ(xz)2/σ2z<σ2x MSE=σx2–σ(xz)2/σz2<σx2

由上式可以得到以下结论：

1. x^=E(X|Z=z)=a+bz x^=E(X|Z=z)=a+bz 是一个  z z 的线性函数；
2. z z 的引入使得 MSE 下降；
3. 如果  X X,  Z Z 互不相关，那么  σxz=0 σxz=0。这时  x^=E(X) x^=E(X)； MSE=σ2x MSE=σx2，这个前文的推论一致。

在  X X 与  Z Z 非 Gaussian 的情況下要得到  x^=E(X|Z=z) x^=E(X|Z=z) 就比较困难。上述 Gaussian 结论需要修正，这里不做展开，但第 2 和第 3 条依然满足。

线性观测的 MMSE

考虑下述特殊情况， Z Z 是  X X 的线性变换加 Gaussian 噪声：

z=Hx+v,x∼N(E(X),σ2x),v∼N(E(V),σ2v) z=Hx+v,x∼N(E(X),σx2),v∼N(E(V),σv2)

这时  X X 的估计为

x^=E(X)+K(z−E(Z)) x^=E(X)+K(z−E(Z))

其中  E(Z)=HE(X)+E(V) E(Z)=HE(X)+E(V)， K=σ2xHT(Hσ2xHT+σ2v)−1 K=σx2HT(Hσx2HT+σv2)−1， K K 称作 Kalman 增益。MSE 为

MSE=σ2x–KHσ2x=σ2x(1–KH) MSE=σx2–KHσx2=σx2(1–KH)

同样的，MSE 由于观测信息  z z 而减小。上式形成了 KF 的基础。

线性观测序列的 MMSE

接下来考虑一个数据流的序贯估计过程，可观测数据流为  z1,z2,z3,…zk z1,z2,z3,…zk， k k 表示离散时刻。测量过程为

zk=Hxk+vk zk=Hxk+vk

k k 时刻的估计由 k k时刻及其之前的观测数据决定：

x^k=E(X|z1,z2,z3,…,zk)=E(X|z1|z2|…|zk) x^k=E(X|z1,z2,z3,…,zk)=E(X|z1|z2|…|zk)

当 k=1 k=1 时，

x^1=E(X|z1)=E(X)+K1(z1–E(Z)) x^1=E(X|z1)=E(X)+K1(z1–E(Z))

MSE1=P1=σ2xHT(Hσ2xHT+σ2v)−1 MSE1=P1=σx2HT(Hσx2HT+σv2)−1

K1=σ2xHT(Hσ2xHT+σ2v)−1 K1=σx2HT(Hσx2HT+σv2)−1

当  k=2 k=2 时，

x^2=E(X|z1|z2)=E(X|z1)+K2(z2–E(Z2))=x^1+K2(z2–E(Z2)) x^2=E(X|z1|z2)=E(X|z1)+K2(z2–E(Z2))=x^1+K2(z2–E(Z2))

MSE2=P2=P21HT(HP21HT+σ2v)−1 MSE2=P2=P12HT(HP12HT+σv2)−1

K2=σ2xHT(Hσ2xHT+σ2v)−1 K2=σx2HT(Hσx2HT+σv2)−1

当  k=3 k=3 时，

x^3=x^2+K3(z3–E(Z3)) x^3=x^2+K3(z3–E(Z3))

…

其实上述公式就是 KF 的雏形。只要再加上动态模型

xk+1=Fxk+w xk+1=Fxk+w

就变成了 KF。反之 KF 的状态方程简化为

xk+1=xk xk+1=xk

就变成了序贯 MMSE 估计器。

序贯 MMSE 估计器和 KF 的一个很大的不同是没有状态过程噪声，因此测量值越多，MSE 会越小，最终接近于 0；KF 有状态过程噪声，在预测过程中，MSE 噪声会增大，在更新过程中，MSE 才会减小，因此 KF 的 MSE 最终会收敛到一个非零的稳态。

Kalman 滤波

Kalman 滤波的预测与更新过程

Kalman 滤波器用反馈控制的方法估计过程状态：滤波器估计过程某一时刻的状态，然后以含噪声的测量变量的方式获得反馈。Kalman 滤波可以分成两部分，更新与预测。

更新过程：
状态更新：

x^k=x¯k+Kk(zk−Hkx¯k) x^k=x¯k+Kk(zk−Hkx¯k)

Kalman 增益：

Kk=PkHTk(Rk+HkPkHTk)−1 Kk=PkHkT(Rk+HkPkHkT)−1

误差协方差更新：

Mk=Pk−KkHkPk Mk=Pk−KkHkPk

预测过程：
误差协方差预测：

Pk+1=FkMkFTk+Qk Pk+1=FkMkFkT+Qk

状态预测：

x¯k+1=Fkx^k x¯k+1=Fkx^k

上式中，更新方程负责反馈，先将先验估计和新的测量变量结合以构造改进的后验估计，预测方程负责推算当前状态变量和协方差估计的值，以便为下一个时间状态构造先验估计。更新方程也可视为校正方程，而预测方程可视为预估方程。

KF 是一种递归的估计，即只要获知上一时刻状态的估计值以及当前状态的观测值就可以计算出当前状态的估计值，因此不需要记录观测或者估计的历史信息。

下图说明了 Kalman 滤波器的递归过程：

Kalman 滤波的一个例子

下面用一个简单的例子说明 KF 的过程。假设待估计值为温度，初值为 23 度，随时间变化。KF 过程的 MATLAB 程序如下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

% Kalman滤波的一个简单的例子 clear all; close all; % 初始值 nT = 100; sz = [nT, 1]; Q = 0.1; % 状态噪声协方差 R = 1; % 测量噪声协方差 x(1) = 23; % 状态真值 for k=2:nT x(k,1) = x(k-1,1) + sqrt(Q)*randn(1); end z = x + sqrt(R)*randn(sz); % 观测值 (含噪声) % 分配存储空间 xhat = zeros(sz); % 后验估计 M = zeros(sz); % 后验估计误差 xhatminus = zeros(sz); % 先验估计 P = zeros(sz); % 先验估计误差 K = zeros(sz); % Kalman增益 % 初始值估计 xhat(1) = rand(1,1)*50; P(1) = 1.0; for k = 2:nT % 预测 xhatminus(k) = xhat(k-1); P(k) = M(k-1)+Q; % 更新 K(k) = P(k)/( P(k)+R ); xhat(k) = xhatminus(k)+K(k)*(z(k)-xhatminus(k)); M(k) = (1-K(k))*P(k); end figure; hold on; plot(x,'r-', 'LineWidth', 2); plot(z,'k+', 'LineWidth', 2); plot(xhat,'b-', 'LineWidth', 2); legend('真值', '测量值', '估计值'); xlabel('时间'); ylabel('温度'); hold off;

运行结果如下图：

参考文献：
[1] Kalman R E. A new approach to linear filtering and prediction problems [J]. Journal of Fluids Engineering, 1960, 82(1): 35–45.
[2] http://allenluadvance.blogspot.jp/2014/01/kalman-filter-mmse-estimator.html

链接： http://lucky6.coding.me/2016/08/19/kalman-filter/