[科普中國]-卡爾曼方程

卡爾曼濾波器是一種由卡爾曼（Kalman）提出的用于時變線性系統(tǒng)的遞歸濾波器。這個系統(tǒng)可用包含正交狀態(tài)變量的微分方程模型來描述，這種濾波器是將過去的測量估計(jì)誤差合并到新的測量誤差中來估計(jì)將來的誤差。1

簡介Dr Kalman 的卡爾曼濾波器。下面的描述，會涉及一些基本的概念知識，包括概率（Probability），隨機(jī)變量（Random Variable），高斯或正態(tài)分配（Gaussian Distribution）還有State-space Model等等。但對于卡爾曼濾波器的詳細(xì)證明，這里不能一一描述。
首先，要引入一個離散控制過程的系統(tǒng)。該系統(tǒng)可用一個線性隨機(jī)微分方程（Linear Stochastic Difference equation）來描述：
X(k)=A X(k-1)+B U(k)+W(k)
再加上系統(tǒng)的測量值：
Z(k)=H X(k)+V(k)

上兩式子中，X(k)是k時刻的系統(tǒng)狀態(tài)，U(k)是k時刻對系統(tǒng)的控制量。A和B是系統(tǒng)參數(shù)，對于多模型系統(tǒng)，他們?yōu)榫仃?。Z(k)是k時刻的測量值，H是測量系統(tǒng)的參數(shù)，對于多測量系統(tǒng)，H為矩陣。W(k)和V(k)分別表示過程和測量的噪聲。他們被假設(shè)成高斯白噪聲(White Gaussian Noise)，COVariance 分別是Q，R（這里假設(shè)他們不隨系統(tǒng)狀態(tài)變化而變化）。

估算對于滿足上面的條件(線性隨機(jī)微分系統(tǒng)，過程和測量都是高斯白噪聲)，卡爾曼濾波器是最優(yōu)的信息處理器。下面結(jié)合covariances 來估算系統(tǒng)的最優(yōu)化輸出。

首先我們要利用系統(tǒng)的過程模型，來預(yù)測下一狀態(tài)的系統(tǒng)。假設(shè)某刻的系統(tǒng)狀態(tài)是k，根據(jù)系統(tǒng)的模型，可以基于系統(tǒng)的上一狀態(tài)而預(yù)測出某刻狀態(tài)：

X(k|k-1)=A X(k-1|k-1)+B U(k)………(1)
式(1)中，X(k|k-1)是利用上一狀態(tài)預(yù)測的結(jié)果，X(k-1|k-1)是上一狀態(tài)最優(yōu)的結(jié)果，U(k)為某刻狀態(tài)的控制量，如果沒有控制量，它可以為0。
到某刻為止，系統(tǒng)結(jié)果已經(jīng)更新了，可是，對應(yīng)于X(k|k-1)的covariance還沒更新。用P表示covariance：

P(k|k-1)=A P(k-1|k-1) A’+Q………(2)

式(2)中，P(k|k-1)是X(k|k-1)對應(yīng)的covariance，P(k-1|k-1)是X(k-1|k-1)對應(yīng)的covariance，A’表示A的轉(zhuǎn)置矩陣，Q是系統(tǒng)過程的covariance。式子1，2就是卡爾曼濾波器5個公式當(dāng)中的前兩個，也就是對系統(tǒng)的預(yù)測。

某刻有了某刻狀態(tài)的預(yù)測結(jié)果，然后再收集某刻狀態(tài)的測量值。結(jié)合預(yù)測值和測量值，可以得到某刻狀態(tài)(k)的最優(yōu)化估算值X(k|k)：

X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H X(k|k-1)) ……… (3)
其中Kg為卡爾曼增益(Kalman Gain)：
Kg(k)= P(k|k-1) H’ / (H P(k|k-1) H’ + R) ……… (4)
到某刻為止，已經(jīng)得到了k狀態(tài)下最優(yōu)的估算值X(k|k)。但是為了要另卡爾曼濾波器不斷的運(yùn)行下去直到系統(tǒng)過程結(jié)束，還要更新k狀態(tài)下X(k|k)的covariance：
P(k|k)=（I-Kg(k) H）P(k|k-1) ……… (5)
其中I 為1的矩陣，對于單模型單測量，I=1。當(dāng)系統(tǒng)進(jìn)入k+1狀態(tài)時，P(k|k)就是式子(2)的P(k-1|k-1)。這樣，算法就可以自回歸的運(yùn)算下去。
卡爾曼濾波器的原理基本描述了，式子1，2，3，4和5就是他的5 個基本公式。根據(jù)這5個公式，可以很容易的實(shí)現(xiàn)計(jì)算機(jī)的程序。2

舉例下面，用程序舉一個實(shí)際運(yùn)行的例子。

舉一個非常簡單的例子來說明卡爾曼濾波器的工作過程。所舉的例子是進(jìn)一步描述第二節(jié)的例子，而且還會配以程序模擬結(jié)果。

把房間看成一個系統(tǒng)，然后對這個系統(tǒng)建模。當(dāng)然，見的模型不需要非常地精確。所知道的這個房間的溫度是跟前一時刻的溫度相同的，所以A=1。沒有控制量，所以U(k)=0。因此得出：

X(k|k-1)=X(k-1|k-1) ……….. (6)
式子（2）可以改成：
P(k|k-1)=P(k-1|k-1) +Q ……… (7)
因?yàn)闇y量的值是溫度計(jì)的，跟溫度直接對應(yīng)，所以H=1。式子3，4，5可以改成以下：
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-X(k|k-1)) ……… (8)
Kg(k)= P(k|k-1) / (P(k|k-1) + R) ……… (9)
P(k|k)=（1-Kg(k)）P(k|k-1) ……… (10)
模擬一組測量值作為輸入。假設(shè)房間的真實(shí)溫度為25度，模擬了200個測量值，這些測量值的平均值為25度，但是加入了標(biāo)準(zhǔn)偏差為幾度的高斯白噪聲（在圖中為藍(lán)線）。
為了令卡爾曼濾波器開始工作，需要告訴卡爾曼兩個零時刻的初始值，是X(0|0)和P(0|0)。他們的值不用太在意，隨便給一個就可以了，因?yàn)殡S著卡爾曼的工作，X會逐漸的收斂。但是對于P，一般不要取0，因?yàn)檫@樣可能會令卡爾曼完全相信你給定的X(0|0)是系統(tǒng)最優(yōu)的，從而使算法不能收斂。選了X(0|0)=1度，P(0|0)=10。
該系統(tǒng)的真實(shí)溫度為25度，圖中用黑線表示。圖中紅線是卡爾曼濾波器輸出的最優(yōu)化結(jié)果（該結(jié)果在算法中設(shè)置了Q=1e-6，R=1e-1）。

方程matlab下面的kalman濾波程序：

clear

N=200;

w(1)=0;

w=randn(1,N)

x(1)=0;

a=1;

for k=2:N;

x(k)=a*x(k-1)+w(k-1);

end

V=randn(1,N);

q1=std(V);

RVV=q1.^2;

q2=std(x);

Rxx=q2.^2;

q3=std(w);

Rww=q3.^2;

c=0.2;

Y=c*x+V;

p(1)=0;

s(1)=0;

for t=2:N;

p1(t)=a.^2*p(t-1)+Rww;

b(t)=c*p1(t)/(c.^2*p1(t)+Rvv);

s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1));

p(t)=p1(t)-c*b(t)*p1(t);

end

t=1:N;

plot(t,s,'r',t,Y,'g',t,x,'b');

function [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V,varargin)

% Kalman filter.

% [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, ...)

%

% INPUTS:

% y(:,t) - the observation at time t

% A - the system matrix

% C - the observation matrix

% Q - the system covariance

% R - the observation covariance

% init_x - the initial state (column) vector

% init_V - the initial state covariance

%

% OPTIONAL INPUTS (string/value pairs [default in brackets])

% 'model' - model(t)=m means use params from model m at time t [ones(1,T) ]

% In this case, all the abovematricestake an additional final dimension,

%i.e., A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m).

% However, init_x and init_V are independent of model(1).

% 'u' - u(:,t) the control signal at time t [ [] ]

% 'B' - B(:,:,m) the input regression matrix for model m

%

% OUTPUTS (where X is the hidden state being estimated)

% x(:,t) = E[X(:,t) | y(:,1:t)]

% V(:,:,t) = Cov[X(:,t) | y(:,1:t)]

% VV(:,:,t) = Cov[X(:,t), X(:,t-1) | y(:,1:t)] t >= 2

% loglik = sum{t=1}^T log P(y(:,t))

%

% If an input signal is specified, we also condition on it:

% e.g., x(:,t) = E[X(:,t) | y(:,1:t), u(:, 1:t)]

% If a model sequence is specified, we also condition on it:

% e.g., x(:,t) = E[X(:,t) | y(:,1:t), u(:, 1:t), m(1:t)]

[os T] = size(y);

ss = size(A,1); % size of state space

% set default params

model = ones(1,T);

u = [];

B = [];

ndx = [];

args = varargin;

nargs = length(args);

for i=1:2:nargs

switch args

case 'model', model = args{i+1};

case 'u', u = args{i+1};

case 'B', B = args{i+1};

case 'ndx', ndx = args{i+1};

otherwise, error(['unrecognized argument ' args])

end

end

x =zeros(ss, T);

V = zeros(ss, ss, T);

VV = zeros(ss, ss, T);

loglik = 0;

for t=1:T

m = model(t);

if t==1

%Prevx= init_x(:,m);

%prevV = init_V(:,:,m);

prevx = init_x;

prevV = init_V;

initial = 1;

else

prevx = x(:,t-1);

prevV = V(:,:,t-1);

initial = 0;

end

if isempty(u)

[x(:,t), V(:,:,t), LL, VV(:,:,t)] = ...

kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, 'initial', initial);

else

if isempty(ndx)

[x(:,t), V(:,:,t), LL, VV(:,:,t)] = ...

kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, ...

'initial', initial, 'u', u(:,t), 'B', B(:,:,m));

else

i = ndx;

% copy over all elements; only some will get updated

x(:,t) = prevx;

prevP = inv(prevV);

prevPsmall = prevP(i,i);

prevVsmall = inv(prevPsmall);

[x(i,t), smallV, LL, VV(i,i,t)] = ...

kalman_update(A(i,i,m), C(:,i,m), Q(i,i,m), R(:,:,m), y(:,t), prevx(i), prevVsmall, ...

'initial', initial, 'u', u(:,t), 'B', B(i,:,m));

smallP = inv(smallV);

prevP(i,i) = smallP;

V(:,:,t) = inv(prevP);

end

end

loglik = loglik + LL;

end

本詞條內(nèi)容貢獻(xiàn)者為:

王寧 - 副教授 - 西南大學(xué)