6/16（月）

信号解析第９回講義録
日時：２０１４年６月１６日
講義内容：状態空間モデルとカルマンフィルタ
担当者：情報知能工学科小島史男
1
はじめに
先週は AR モデルの推定について学習しました。データの AR モデルによるあてはめがレビンソンのア
ルゴリズムを用いて効率よく実現できることを示しました。教科書では第８章でさらに詳細な議論があ
りますが、講義では時間の関係でとばします。今週は状態空間モデルによる推定問題について考察しま
す。カルマンフィルタはシステムの状態量の推定値を逐次求めていく計算アルゴリズムです。ここでは、
データがモデルのなかに時時刻々組み込まれていきます。教科書第９章の内容です。
2
状態空間モデル
今回の講義では混乱を避けるため、取得される時系列信号は
YN = {y[1] = y˜1 , y[2] = y˜2 , · · · , y[N ] = y˜N }
のようにデータ長 N と冠記号をつけて、明示的に区別します。またこれまで、時系列信号はスカラー量
としてきましたが、ここでは観測データは各時刻で l 次元のベクトル表現
(
y[n] = y 1 [n], y 2 [n], · · · y l [n]
)T
∈ Rl
をとることとします。この時系列 y[n] のモデルを記述します。まずシステムの状態量 x[n] は k 次元の状
態ベクトル
(
)
x[n] = x1 [n], x2 [n], · · · xk [n]
T
∈ Rk
で与え、システムの不規則性は、m 次元の互いに独立な正規型白色雑音ベクトル
(
)T
v[n] = v 1 [n], v 2 [n], · · · v m [n]
∈ Rm
で与えます。ただし、平均値ベクトルは 0 また分散共分散行列は Q[n] ∈ Rn×m で与えられるとします。
このときシステムの状態量 x[n] は次の逐次更新式で与えられるものとします。
x[n] = F [n]x[n − 1] + G[n]v[n]
ただし、 F [n] および G[n] は時刻 n におけるそれぞれ k × k および k × m の推移行列、すなわち

F [n] =






G[n] =






F 12 [n] · · ·
F 22 [n] · · ·
..
..
.
.
F k1 [n] F k2 [n] · · ·
F 1k [n]
F 2k [n]
..
.
G12 [n] · · ·
G22 [n] · · ·
..
..
.
.
Gk1 [n] Gk2 [n] · · ·
G1m [n]
G2m [n]
..
.
F 11 [n]
F 21 [n]
..
.
G11 [n]
G21 [n]
..
.





F kk [n]
∈ Rk×k






∈ Rk×m
Gkm [n]
とします。これを状態空間モデルと呼びます。一方観測信号 YN = {˜
y [n]}N
n=1 に対応する観測モデルは、
この状態ベクトル x[n] と観測に混入する雑音を l 次元の互いに独立な正規型白色雑音ベクトル
(
)T
w[n] = w1 [n], w2 [n], · · · , wl [n]
との和によって
y[n] = H[n]x[n] + w[n]
で記述することにします。ただし w[n] は平均値ベクトルは 0、かつ分散共分散行列が R[n] ∈ Rk×l で与
えられるとし、H[n] は時刻 n における l × k の観測行列、すなわち



H[n] = 


H 12 [n] · · ·
H 22 [n] · · ·
..
..
.
.
H k1 [n] H k2 [n] · · ·
H 11 [n]
H 21 [n]
..
.
H 1l [n]
H 2l [n]
..
.






∈ Rk×l
H kl [n]
で記述します。すなわちシステム全体の状態量に関する情報の次元 k に対して l 次元の観測情報に縮約
されたことを意味します。結局状態空間モデルは
x[n] = F [n]x[n − 1] + G[n]v[n]
y[n] = H[n]x[n] + w[n]
で記述される構造を持っています。これが線形システムの標準系であり、線形システム解析の基本モデ
ルとなっています。このシステムの入力はシステムノイズ v[n] およびこれとは独立な観測ノイズ w[n] で
あり、モデルの出力が y[n] に相当します。モデルの基本構成はシステムの次元 (k, m, l) およびダイナミ
クスと観測機構を定める行列 (F [n], G[n], H[n]) となります。それではこのモデルと時系列信号との関係
づけを以降おこなっていきます。
3
カルマンフィルタの構成
ここでは求めるのは、観測データ YN = {˜
y1 , y˜2 , · · · , y˜N } にもとづいて状態量 x[n] のフィルタリングを
行います。カルマンフィルタ (Kalman filter) のアルゴリズムの導出には各種の方法がありますが、詳し
くは講義科目「確率過程論」で説明されるとおもいますので、ここでは省略します。フィルタの特性を
求めるために以下の記法を用いることにします。
x
ˆn|n := E(x[n]|Yn )
Vn|n := E{(x[n] − x
ˆn|n )(x[n] − x
ˆn|n )T |Yn )
ここで、Yn = {˜
yk }nk=1 は時刻 n までの観測量の集積とします。カルマンフィルタは以下の予測問題と
フィルタリングの問題を各時間ステップで繰り返し計算していくアルゴリズムです。
予測問題（１ステップ）
：
x
ˆn|n−1 = F [n]ˆ
xn−1|n−1
Vn|n−1 = F [n]Vn−1|n−1 F [n]T + G[n]Q[n]G[n]T
濾波問題（フィルタリング）
：
(
)−1
K[n] = Vn|n−1 H[n]T H[n]Vn|n−1 H[n]T + R[n]
(
x
ˆn|n = x
ˆn|n−1 + K[n] y˜n − H[n]ˆ
xn|n−1
)
Vn|n = (I − K[n]H[n]) Vn|n−1
(
)
観測データの情報は予測誤差の計算 y˜[n] − H[n]ˆ
xn|n−1 に組み込まれています。この予測誤差のこと
をイノベーション過程と呼びます。この修正はカルマンゲイン K[n] によって制御されています。更新式
を少し変形すると
x
ˆn|n = K[n]˜
yn + (I − K[n]H[n])ˆ
xn|n−1
とも書けます。つまり、フィルタリング x
ˆn|n は新しい観測データ y˜[n] と予測ベクトル x
ˆn|n−1 の線形和
となっています。また分散共分散行列についてみると
Vn|n = Vn|n−1 − K[n]H[n]Vn|n−1
となって、誤差が減少する方向、つまり繰り返し演算を通じて精度が改善されいく様子がみてとれます。
4
カルマンフィルタを動かすには
このフィルタを実行するには下記の情報を事前に与える必要があることがわかります。
(1) 信号のデータ長 N 、状態空間表示の次元 (k, m, l)
(2) 状態空間を構成する各時刻 n における行列要素
F [n] ∈ Rk×k ,
G[n] ∈ Rk×m ,
H[n] ∈ Rl×k
(n = 1, · · · , N )
(3) 各時刻 n における不規則雑音の分散共分散行列
Q[n] ∈ Rm×m ,
R[n] ∈ Rl×l
(n = 1, · · · , N )
(4) 初期状態の平均値ベクトルと２乗の平均値行列
x
ˆ0|0 = E(x[0])
V0|0 = E{(x[0] − x
ˆ0|0 )(x[0] − x
ˆ0|0 )T }
= E(x[0]x[0]T ) − x
ˆ0|0 x
ˆT0|0
(5) 観測信号 YN = {˜
yi } N
i=1
ところで、前回、前々回に学習した ARMA モデルは本日示した状態空間モデルの構造を持っている
のです。そのとき行列の構造 F [n]、G[n]、H[n] には係数パラメータや白色雑音の分散 σ 2 が含まれます。
どのように含まれるかについては来週説明することにします。
宿題：先週は AR モデルの推定について学習しました。来週のカルマンフィルタの導入の理解を深める
ために、講義予定表にある時系列信号をダウンロードして、レビンソンのアルゴリズムを用いて AR モ
デルを構成してみましょう。先週の講義録にある R のプログラムを参考にしてください。
(1) データの自己相関関数、ペリオドグラムをみて AR モデルに適合するか検討してください。
(2) AIC 基準から、適切な次数およびそのときの係数を求めてください
(3) 近似された AR モデルのシミュレーション結果
（これはレポートごとに異なるはずです。念のため）ともとの時系列信号とを比較してみてください。
補遺：本日紹介いたしました、Kalman の論文の出典は
R.E. Kalman, A New Approach to Linear Filtering and Prediction Problems, ASME Transactions, Volume 82, Part D (Journal of Basic Engineering), 1960, pp. 35-45
です。図書館には当該論文が保管されていますので、閲覧可能です。

Download Report