カテゴリー: 大学数学数学・統計学

カルマンフィルターと直交射影

カルマンフィルターは、観測できない確率変数(潜在変数)があるモデル(状態空間モデル)において、観測可能な確率変数に基づいて観測不可能な変数の状態を推定する手法です。カルマンフィルターは行列計算を用いて初等的に導出することも可能ですが、計算が結構複雑になる上に、逆行列の存在を仮定しなければならない部分もあったりと難点もあります。この記事では、ヒルベルト空間における直交射影の理論を用いてカルマンフィルターの漸化式を導出します。こちらのほうが理論がすっきりしており、数学の知識がある方には理解しやすいと思います。

問題設定

時刻k=0,1,2,に対し、xk観測できない状態量yk観測可能な観測値uk既知の入力ベクトルとし、それらが以下の条件を満たしているとします。 {xk+1=Akxk+Bkuk+ξkyk=Ckxk+Dkuk+ηk  (k=1,2,)
ここで、

  • 状態変数xkおよびノイズξkRn値確率変数、 観測ベクトルykおよび ノイズηkRp値確率変数 、入力ベクトルukRm
  • AkMn(R),BkMm,n(R),CkMn,p(R),DkMm,p(R)
  • x0,ξ1,ξ2,,η1,η2,は独立な確率変数で、  E[ξk]=0, V[ξk]=Qk, E[ηk]=0, V[ηk]=Rk
  • さらに、 Rkは正定値であるとする。

この設定の下で、既知のパラメータAk,Bk,Ck,Dk,Qk,Rkおよび観測値y1,y2,,から、xkたちの最適な推定量を構成する方法を考えます。さらに、ここでは初期値x0も既知であるとします。

まずは、入力uk=0の場合を考えます。この場合、xk,ykの満たす方程式は、{xk+1=Akxk+ξkyk=Ckxk+ηk  (k=1,2,)となります。

直交射影による推定量の定義

ここでは、2つの確率変数ベクトルの「内積」もどきを定義し、直交射影を構成します。 以下では、行列Aに対しその転地行列をAで表します。また、ほとんどいたるところ等しい確率変数は同一視することにします。

定義
Rd1値確率変数xRd2値確率変数yに対し、(x,y)およびx(x,y)=E[xy], x=(x,x)12 と定める。

ここで、(x,y)はスカラーではなく行列であることに注意します。 これに対しても、次に示すような内積に似た性質が成り立ちます。

命題1
Rd1値確率変数x,x1,x2Rd2値確率変数y,y1,y2AMn,d1(R),BMn,d2(R)に対し、次が成り立つ。
(1) (y,x)=(x,y)
(2) (x1+x2,y)=(x1,y)+(x2,y), (x,y1+y2)=(x,y1)+(x,y2)
(3) (Ax,y)=A(x,y), (x,By)=(x,y)B
(4) trx=0  x=0

(証明)  (1) (y,x)=E[yx]=E[xy]=(x,y)
(2) (x1+x2,y)=E[(x1x2)y]=E[x1y]+E[x2y]=(x1,y)+(x2,y) であり、第2式も同様である。
(3) (Ax,y)=E[(Axy]=AE[xy]=A(x,y) であり、また(1)より、 (x,By)=(By,x)={B(y,x)}=(x,y)B である。
(4) trx=trE[xx]=trE[xx]=E[|x|2]だから示すべきことを得る。 □

次に、Lk={A0x0+kj=1Ajyk;A0Mn,n(R),AjMn,p(R)}, L1={0} と置きます。
上の命題より、Rn値確率変数全体はtr(,)を内積とするHilbert空間で、Lkはその閉部分空間となるので、直交射影の存在と一意性から次の命題が得られます。

命題2
Rn値確率変数xに対し、yLktryx2=E[|yx|2] を最小にするもの(xLkへの直交射影)が一意的に存在する。
以下では、これをPkxと書く。特にP1x=0である。

これを用いて、「観測値y1,,yjに基づくxkの推定量」xk|jを次のように定義します。

定義
k,j=0,1,2,に対し、ˆxk|jˆxk|j=Pjxk によって定める。

これは、x0およびy1,,yjの線形結合で表されるxkの推定量のうち、平均2乗誤差を最小にするものということです。

入力がない場合のカルマンフィルタ

次に、上で構成した推定量xk|jを具体的に計算するための式を導きます。以下では、時刻が進むたびにy1,y2,たちが順次観測されてゆき、ykが観測された時点で、その時の状態量xkを推定したい状況を考えます。つまり、xk|kを知りたいということです。
この時、時刻kにおいて新たに観測されたykと、 時刻k1で求めたxk1|k1の値を使って、新たな推定量xk|kを計算するような漸化式が得られれば、計算のコストが省けてうれしいわけです。そのような漸化式を導くために、いくつかの補題を証明します。いずれも補題も、直交射影の性質から類推されるものです。

以下記述を簡単にするため、y0=x0と定めます。

補題1
ˆyLkおよびRn値確率変数xに対し、次は同値である。
(1) ˆy=Pkx  
(2) j=0,1,,kに対し、(ˆyx,yj)=O

(証明) (1)(2):背理法で示す。(1)を満たすˆyに対し、ある0jkが存在して(ˆyx,yj)Oであるとする。この時、yj0であるから、tryj2=E[|yj|2]>0である。そこで、 y=ˆy(ˆyx,yj)tryj2yj と置くと、 tryx2=trˆyx(ˆyx,yj)tryj2yj2=trˆyx2tr((ˆyx,yj)tryj2yj,ˆyx)tr(ˆyx,(ˆyx,yj)tryj2yj)+tr(ˆyx,yj)tryj2yj2=trˆyx2tr(ˆyx,yj)(yj,ˆyx)tryj2tr(ˆyx,yj)(ˆyx,yj)tryj2+tr(ˆyx,yj)(ˆyx,yj)tryj2=trˆyx2tr(ˆyx,yj)(ˆyx,yj)tryj2 であるが、(ˆyx,yj)Oから上と同様にtr(ˆyx,yj)(ˆyx,yj)>0となり、(1)の仮定に反する。よって、(2)が成り立つ。
(2)(1):yLkとすると、仮定より(ˆyx,yˆy)=Oであるから、 tryx2=trˆyx+yˆy2=trˆyx2+tryˆy2trˆyx2 となり(1)が成り立つ。 □

上の補題により、次のような結果が得られます。

補題2
(1) 任意のAMn,p(R)およびk=0,1,2,に対し、  PkAx=APkx  が成り立つ。
(2) 任意の0jkに対し、Pjξk=0が成り立つ。

(証明) (1) 任意のj=0,1,,kに対し、  (APkx,yj)=A(Pkx,yj)=O  だから上の補題により示すべきことを得る。
(2) (2) 任意のi=0,1,,jに対し、yiξi1までの式で表されるから、ξkとは独立である。よって、(ξk,yi)=0となるから、上の補題により結論が従う。 □

次に、{yk}の「直交化」であるイノベーション過程を定義します。

定義
zk=ykCkPk1xk=ykCkˆxk|k1  で定まる{zk}イノベーション過程と呼ぶ。特に、P1=Oにより、z0=y0=x0である。

次の補題により、これが{yk}の「直交化」になっていることがわかります。

補題3
イノベーション過程{zk}に対し、  (zi,zj)={CixiPi1xi2Ci+Ri(i=j1)O(ij)Hk={kj=1Ajzk;A0Mn,n(R),AjMn,p(R)}  である。

(証明) zk=Ckxk+ηkCkPk1xk=Ck(xkPk1xk)+ηk  であるから、i=jの時xiPi1xiηiの独立性により、  (zi,zi)=Ci(xiPi1xi)2+ηi2=CixiPi1xi2Ci+Ri  となる。一方i<jとすると、ηjηiおよびxiPi1xiと独立であるから、  (zi,zj)=(Ci(xiPi1xi)+ηi,Cj(xjPj1xj))=(zi,Cj(xjPj1xj)=(yiCiPi1xi,Cj(xjPj1xj))=O  となる。最後の等号は、補題1より(yiCiPi1xi,xjPj1xj)=Oとなることを用いた。よって前半の主張が示された。後半の主張は帰納的に明らかである。 □

上の補題により、k1の時zk2は半正定値+正定値の形なので、正定値となることがわかります。そこで、ek=zk1zkと定めると、{ek}k1{yk}k1の「正規直交化」となり、次の補題が成り立ちます。

補題4
ek=zk1zk (k1)とすると、Rn値確率変数xに対し、  Pjx=P0x+ji=1(x,ei)ei  が成り立つ。

(証明) y=P0x+ji=1(xi,ei)eiと置くと、前の補題よりHjであr任意のei (1ij)に対し、  (xy,ei)=(x,ej)(y,ej)=(x,ei){(P0x,ej)+(x,ei)(ei,ei)}=0(xy,x0)=(P0xy,x0)=0  となる。これと前の補題により、任意のyi (0ij)に対し、(xy,yi)=0となるから、補題1によりy=Pjxである。 □

以上の準備の下で、xk|kを計算するための漸化式が証明できます。

定理5
k1に対し、Gk=(xk,ek)zk1と置くと、  ˆxk|k1=Ak1ˆxk1|k1ˆxk|k=ˆxk|k1+Gk(ykCkˆxk|k1)  が成り立つ。

(証明) 補題2により、  ˆxk|k1=Pk1xk=Pk1(Ak1xk1+ξk1)=Ak1Pk1xk1=Ak1ˆxk1|k1  となるから最初の式が従う。
また、補題4により  ˆxk|k=Pkxk=Pk1xk+(xk,ek)ek=xk|k1+(xk,ek)zk1zk=xk|k1+(xk,ek)zk1(ykCkˆxk|k1)  となり2式目も従う。

上で登場したGkカルマンゲインと呼ばれます。これを計算するための漸化式を求めます。

定理6
Pk,j=xkˆxk|j2とすると、  Pk,k1=Ak1Pk1,k1Ak1+Qk1Pk,k=(IGkCk)Pk,k1Gk=Pk,k1Ck(CkPk,k1Ck+Rk)1  が成り立つ。

(証明) 前の定理を用いると、(xk1ˆxk1|k1)ξk1の独立性により、  Pk,k1=xkˆxk|k12=A(xk1ˆxk1|k1)+ξk12=A(xk1ˆxk1|k1)2+ξk12=Ak1Pk1,k1Ak1+Qk1  となり最初の式が得られる。
また、ηkxkの独立性および補題7から(xk,ηk)=0,(xk|k1,xkxk|k1)=0となることに注意すると、  Gk=(xk,ek)zk1=(xk,zk)zk2=(xk,Ck(xkxk|k1)+ηk)zk2=(xk,Ck(xkxk|k1))zk2=(xkxk|k1,Ck(xkxk|k1))zk2=xkxk|k12Ck(zk,zk)1=Pk,k1Ck(CkPk,k1Ck+Rk)1  となり、3番目の式を得る。
  さらにこれより、  Pk,k=xkˆxk|k2=xkˆxk|k1Gk(ykCkˆxk|k1)2=xkˆxk|k1Gk(CkxkCkˆxk|k1+ηk)2=(IGkCk)(xkˆxk|k1)2+Gkηk2=(IGkCk)Pk,k1(IGkCk)+GkRkGk=(IGkCk)Pk,k1Pk,k1CkGk+GkCkPk,k1CkGk+GkRkGk=(IGkCk)Pk,k1Gk(CkPk,k1Ck+Rk)Gk+GkCkPk,k1CkGk+GkRkGk=(IGkCk)Pk,k1  となって、2番目の式を得る。 □

以上により、xk|kを計算するための漸化式が得られました。定義により初期推定量x0|0x0|0=P0x0|0=x0となることも合わせると、次のような結果になります。

カルマンフィルターによる状態推定(入力uk=0の場合)
Pk,k1=Ak1Pk1,k1Ak1+Qk1Pk,k=(IGkCk)Pk,k1Gk=Pk,k1Ck(CkPk,k1Ck+Rk)1ˆx0|0=x0ˆxk|k=ˆxk|k1+Gk(ykCkˆxk|k1)ˆxk|k1=Ak1ˆxk1|k1

入力がある場合

次にukが一般の場合を考えます。この時は{x(1)0=x0,x(1)k+1=Akx(1)k+ξkx(1)1=0,x(2)k+1=Akx(2)k+Ckuk, {y(1)k=Bkx(1)k+ηky(2)k=Bkx(2)k+Dkukと定めると、xk=x(1)k+x(2)k, yk=y(1)k+y(2)kとなります。今x(2)kについては直接計算することができますので、x(1)kの値を推定できれば十分です。そこで、前節でxkx(1)k,yky(1)kとして得られる推定量をˆx(1)k|jとし、ˆxk|j=ˆx(1)k|j+x(2)kと定めると、定理5のGkに対し、ˆxk|k=ˆx(1)k|k+x(2)k=ˆx(1)k|k1+Gk(y(1)kCkˆx(1)k|k1)+x(2)k=ˆxk|k1+Gk{(yky(2)k)Ck(ˆxk|k1ˆx(2)k|k1)}=ˆxk|k1+Gk(ykCkˆxk|k1Dkuk)ˆxk|k1=ˆx(1)k1|k+x(2)k=Ak1ˆx(1)k1|k1+x(2)k=Ak1(ˆxk1|k1x(2)k1)+Ak1x(2)k1+Ck1uk1=Ak1ˆxk1|k1+Ck1uk1となります。以上をまとめると、次の結果になります。

カルマンフィルターによる状態推定(入力\bm u_kがある場合)
\begin{align*}  &P_{k,k-1}=A_{k-1}P_{k-1,k-1}A_{k-1}’+Q_{k-1}\\  &P_{k,k}=(I-G_kC_k)P_{k,k-1}\\  &G_k=P_{k,k-1}C_k'(C_kP_{k,k-1}C_k’+R_k)^{-1}\\  &\hat{\bm x}_{0|0}=\bm x_0\\  &\hat{\bm{x}}_{k|k}=\hat{\bm{x}}_{k|k-1}+G_k(\bm{y}_k-C_k\hat{\bm{x}}_{k|k-1}-D_k\bm u_k)\\  &\hat{\bm{x}}_{k|k-1}=A_{k-1}\hat{\bm{x}}_{k-1|k-1}+C_{k-1}\bm u_{k-1}\end{align*}

P_{k,j}G_kは入力があってもなくても同じ値になります。さらに、これらは観測データ\bm y_kがなくても事前に計算することができます。そのため、G_kをあらかじめ計算しておき、\bm y_kが観測された時点で\bm x_kを計算することにより、計算リソースを節約することができます。

+4

作成者:

数学科の学生で、確率論、統計学を専攻しています。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です