卡爾曼方程

卡爾曼濾波器是一種由卡爾曼（Kalman）提出的用於時變線性系統的遞歸濾波器。這個系統可用包含正交狀態變數的微分方程模型來描述，這種濾波器是將過去的測量估計誤差合併到新的測量誤差中來估計將來的誤差。

基本介紹

• 中文名：卡爾曼方程
• 外文名：Kalman equation
• 人物：Dr Kalman
• 儀器：卡爾曼濾波器
• 公式：　X(k)=A X(k-1)+B U(k)+W(k)

簡介

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

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

估算

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

首先我們要利用系統的過程模型，來預測下一狀態的系統。假設某刻的系統狀態是k，根據系統的模型，可以基於系統的上一狀態而預測出某刻狀態：

X(k|k-1)=A X(k-1|k-1)+B U(k)………(1)
式(1)中，X(k|k-1)是利用上一狀態預測的結果，X(k-1|k-1)是上一狀態最優的結果，U(k)為某刻狀態的控制量，如果沒有控制量，它可以為0。
到某刻為止，系統結果已經更新了，可是，對應於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)對應的covariance，P(k-1|k-1)是X(k-1|k-1)對應的covariance，A’表示A的轉置矩陣，Q是系統過程的covariance。式子1，2就是卡爾曼濾波器5個公式當中的前兩個，也就是對系統的預測。

某刻有了某刻狀態的預測結果，然後再收集某刻狀態的測量值。結合預測值和測量值，可以得到某刻狀態(k)的最最佳化估算值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)
到某刻為止，已經得到了k狀態下最優的估算值X(k|k)。但是為了要另卡爾曼濾波器不斷的運行下去直到系統過程結束，還要更新k狀態下X(k|k)的covariance：
P(k|k)=（I-Kg(k) H）P(k|k-1) ……… (5)
其中I 為1的矩陣，對於單模型單測量，I=1。當系統進入k+1狀態時，P(k|k)就是式子(2)的P(k-1|k-1)。這樣，算法就可以自回歸的運算下去。
卡爾曼濾波器的原理基本描述了，式子1，2，3，4和5就是他的5 個基本公式。根據這5個公式，可以很容易的實現計算機的程式。

舉例

下面，用程式舉一個實際運行的例子。

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

把房間看成一個系統，然後對這個系統建模。當然，見的模型不需要非常地精確。所知道的這個房間的溫度是跟前一時刻的溫度相同的，所以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)
因為測量的值是溫度計的，跟溫度直接對應，所以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)
模擬一組測量值作為輸入。假設房間的真實溫度為25度，模擬了200個測量值，這些測量值的平均值為25度，但是加入了標準偏差為幾度的高斯白噪聲（在圖中為藍線）。
為了令卡爾曼濾波器開始工作，需要告訴卡爾曼兩個零時刻的初始值，是X(0|0)和P(0|0)。他們的值不用太在意，隨便給一個就可以了，因為隨著卡爾曼的工作，X會逐漸的收斂。但是對於P，一般不要取0，因為這樣可能會令卡爾曼完全相信你給定的X(0|0)是系統最優的，從而使算法不能收斂。選了X(0|0)=1度，P(0|0)=10。
該系統的真實溫度為25度，圖中用黑線表示。圖中紅線是卡爾曼濾波器輸出的最最佳化結果（該結果在算法中設定了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