宣告第2

= rり

5 宣告第2

移流方程式計からのデータをもとにした降雨予測では

からレーダ雨従

ング手法を適用して予や移流拡散方程式にカルマンフィルターなどのフィルタリ

降水レーダを用いた水文現象の予測手法に関する研究

測が行われてきた。しかし、空間差分を用いる場合に降雨域が孤立波的形状で無降雨の部分が含まれるためフィルターの動作が不安定になりやすく、レーダの観測範囲のほとんどで降雨が観測されていないと計算が発散してしまう問題があった。このため、移流方程式の拡散項を無視したり、流速の空間分布を一次式で近似したりしてその物理性に問題があった。また、各観測点毎に雨量を計算する必要があった。

そこで、レーダのデータのように孤立波的形状をしている波形を計算するために、移流拡散方程式を波数空間に展開し、フィルタリング手法の適用を試みる。

まず最初は一次元のレーダデータに対して本手法を適用し、その精度を検証す

ると共に二次元のデータに対する予測l手法の理論的展開と問題点を述べることにする。

第2 . 5 . 1項カルマンフィルタ一理論

カルマンの予測理論におけるシステム方程式はそれぞれ次のように表すことが

できる2 1)。

x(k+ 1)=φ(k+ l/k)'x(k)+δ(k+ 1)

y(k+ 1 )=M(k+ 1) .x(k+ 1)+ε(k+ 1) ^(2.5.1)

ここで、 x(k)は k ステップでのシステムの真の状態量ベクトル、 φ(k+1jk)は時間ステップ k から k+1への遷移行列、 y (k+1)は状態量の変化を表す観測行列、お

よびS と ε は誤差ベクトjレ(白色性正規雑音)である。 (2.5.1)式の最初の式をシステム方程式、 2番目の式を観測方程式と呼ぶ。

ここで、カルマンフィルターの理論を展開する際の記号の説明を行う。

八は最適推定値及びーは推定誤差である。従って各記号の意味は下表のようになる。

降水レーダを用いた水文現象の予測手法に関する研究

x欣+1) k+1ステップでのシステムの真の状態量

;(k+ ljk) k ステップまでの情報を利用して求められたk+l ステ

ツプのxの最適推定値

Z欣+l jk+1) k+1 ステップまでの情報を利用して求められたその時

間ステップでのxの最適推定値

'i(k+ljk) =x(k+l) - ;(k+1jk) k スアップまでの情報によりk+1 ステップまでのxを

推定した場合の誤差

'X(k/k) =x(k) - ;(k/k) k ステップまでの情報によりその時刻でのxを推定し

た場合の誤差

y(k) k スアップでの観測量

，、、

k ステップまでの観測情報よりの、次のステップでの y (k+ ljk)

y (kjk)の最適推定値

お

か1) =y (k+ 1)・

y

(k+ 1jk) 観測値の推定誤差

ここで時刻 1，2，・・・・，kにおける観測量をy(1) ， y (2) ，・・・・，y (k) とする。われわれは Yiのみを観測しているのであり、時間ステップk+1における真の状態量の最

適推定値は回帰関係式

五(k+ 1/ k) =a( 1 ，k)y( 1 )+a(2，k)

(2)+・・+a(k，k)y(k)

(2.5.2) により求められるとするo k+ 1までの情報より予想される次のステップk + 1 の

観測量の推定値

デ

俳+ljk)は従って(2.5.1)式の観測方程式により

y(k+ l/k)=M(k+ 1 ).Jr7._k+ l/k)

である。

(2.5.3)

ところで、次の時間ステップk+1において、観測値 y(k+1)が得られたとする。

新しく得られた情報 y(k+ 1)はそれ以前までの観測量y(iXi=1，2，' "， k)と従属な成分

デ(k+1/k) (即ち Y (i)の線形結合で作り得る成分)と、それらと独立な成分子(k+1)とに分けられる。

このことから y(k+ 1) は

y(k+ 1 )=M(k+ 1 ).x"'{k + l)+y(k+ l/k)

と表される。即ち

y(k+ 1 )=y(k+ 1)手(k+ 1)

(2.5.4)

(2.5.5)

このデ(k+1)は最も新しく、信頼できる情報である。 k+1ステップにおいて新しい

降水レーダを用いた水文現象の予測手法に関する研究

観測情報が得られた場合の x の最適推定値�(k+1 /k+ 1)は、 k ステップまでの情

報からの x(k+l/k)の推定値;(k+l/k) にこの今までの情報と独立な成分 y(k+1) に重みを付けて加えると良いと考えられる。従って次の関係が仮定される。

五(k+ 1/ k+ 1) =x(k+ 1/ k)+Jf(_沢k+ 1)

_(2.5.6)

この荷重係数 K をKa1man Gain と呼ぶ。

次にk ステップまでの情報によりk+1ステップでの状態量 x(k+1)を推定するときの誤差を x(k+1/k) とする。

x(k+ l/k)=x(k+ 1)ー支(k+l/k)

あるいは

x(k+ l)=x(k+ l/k)+x(k+ l/k)

x (k+1) に変換 M を作用させると、 (2.5.1)式の観測方程式より

y(k+ 1) =M(k+ 1 ).x(k+ 1)+ε(k+ 1)

となる。これに(2.5.8)式を代入し次式を得る。

y(k+ l)=M(k+ l).x(k+ 1 )+M(k+ 1 ).x1_k+ 1)+ε(k+ 1)

これを(2.5.8)式に代入し、 (2.5.1)式の観測方程式を用いて表すと、

y(k+ l)=M(k+ 1 )-x1_k+ 1)+ε(k+ 1)

(2.5.7)

(2.5.8)

(2.5.9)

(2.5.10)

(2.5.11)

となる。即ちk ステップまでのデータで予測しきれなかったk+1 ステップでの予測誤差はk+1ステップでの観測量の中に情報として含まれている。

次に、 k+1ステップまでの情報を得たときの推定誤差を x(k+l/k+1)と表すとき、

(2.5.5)式の関係を用いて次のようになる。

x(k+ 1/ k+ 1) =x(k+ 1 )x(k+ 1/ k+ 1)

=x(k+ 1 )- [ えれl/k)+lK長k+ 1) ]

=x(k+ l/k)-lKy(k+ 1)

(2.5.12)

従って、王(k+1 /k+ 1) は k+lステップまでの観測情報では推定し得ない成分であるo よって、えれl/k+1)=五(k+l/k)-K.ヌk+l )はk+1 ステップでの新しい情報成分 y (k+ 1)に

降水レーダを用いた水文現象の予測l手法に関する研究

も独立である。

(双k+l/k)-!Kyてk+1))上知+ 1)

^(2.5.13)

よって確率変数の独立・直交の性質より、

E[x(k+l/k).yT(k+1)]-!K E [Y(k+1).yT(k+1)]=O

^(2.5.14)

あるいは(2.5.11)式の関係を用いて

抑制/k). (MT(k+l/k).xT(k+l)+εT(k+1))]

^(2.5.15)

={(M(k+ 1 }xTk+ 1) +ε(k+l)). (MT(k+l).xT(k+l)+εT(k+1))]

雑音 ε(k+1)と状態量とは無相関であるから

[玄(k+l/k)王T(k+ 1) .MT(k+川 (2.5凶)

= [M(k+l)沢k+ 1)玉川+l)).MT(k+l)+ε(k+ 1)εT(k+l)]

ここで E

I.x(k+ l/k).;?(k+ 1) I

^はx の推定誤差の相関行列(共分散行列) IP であ

り、観測誤差の分散行列を R とすれば

P(k+ l/k)三E[玄(k+1)玄T(k+ 1)]

R(炉開=E[e(k+ 1)εT(k+ 1)]

(2.5.17)

従ってか1ステップで新しい観測情報y(k+ 1)が得られたとき、最も信頼し得

る情報である「観測量の推定値 y(k+ 1)と実際の観測値 y(k+1)との差デ(k+1)Jによ

る補正係数 K ^は

!K(k+ 1 )=tP(k+ l/k)MT(k+ l)J. [M(k+ l)P(k+ 1/ k)MT (k+ l)+R(k+ 1 )]-1

^(2.5.18)

により与えられる。

次に、推定誤差'X(k+1/k)の共分散行列 lP(k+l/k)について詳しく説明する。

x (k+ l/k)の定義は

正'(k+ 1/ k)=x(k+ 1 )-x(k+ 1/ k)

(2.5.l9)

である。一方訂k/k)は

双k/k)=x(k)-x(k/ k)

(2.5.20)

」こで、

x(k+ 1)=φ(k+ l/k).x(k)+δ(k+ 1)

_(2.5.21)

降水レーダを用いた水文現象の予測手法に関する研究

五(k+ l/k)=�k+ 1/k}x1_k/ k)

^(2.5.22)

の関係を考慮し、 (2.5.20)式の両辺にφを作用させると

φ .xてk/k)=φ.x(k)ーφ.x{_k/k)

^(2.5.23)

=(x(k+

)-�(k+ 1))仰+ l/k)

:. i(k+ 1/ k)=φ.xT.k/k)+δ(k+ 1)

^(2.5.24)

の関係を得る。上式とその転置行列の積により、 lP(k+ l/k)は次のようになる。

P(k+附=E[双k+ l/k).xT.k+ 1 / k)]

^(2.5.25)

=E [x(k+ 1/ k) 双が同T]φl +E [ δ(か1)δ(k+ 1) T ]

:. P(k+ l/k)=φ.P(k/ k) .φ+Q T

^(2.5.26)

、... 、...

1-'- '- 官」

Q= & 伽1)δ(k+ 1) T ]

はシステム雑音の分散行列である。

最後に lP(k+1 /k+ 1)について考える。

(2.5.12)式より

五(k+1/k+1)=支(k+ l/k)-lK只k+1)

(2.5.27)

である。一方k+1ステップでの推定誤差X'(k+1 /k+1) とれlステップでの補正係数

//( y (k+ 1)とは独立である。

五(k+l/k+ l).llKj(k+ 1)

故に(2.5.27)式より

(双k+1/ k)-lK只k+1 )XK<只k+1 ))T =0

k+1ステップでの最適推定誤差の共分散行列 lP(k+1/ k+ 1)は

P(k+1/k+1)=E [双k+1/k+1)気k+ l/k+ 1)]

=再(双k+1/ k)-lK沢k+1)阪k+1/ k)-lK沢k+

)y]

=E l(双k+ l/k)-lK只k+ 1) )X(k+ 1/ k)TJ +E[(沢k+ l/k)lK沢k+1))XlK，ヌk+

)T]

=E[双k+ l/k)玖k+ l/k)TJ -KME[玖k+ l/k)双k+l/k)TJ

(2.5.28)

(2.5.29)

降水レーダを用いた水文現象の予測手法に関する研究

=P(k+l/k) -KMJP(k+ l/k)

=[J/-K MJ!P(k+l/k)

ここに fは単位行列である。

(2.5.18)式、 (2.5.26)式及び(2.5.30)式によりKa1man Gainが逐次決定される。

第2. 5.

2項一次フィルター及び二次フィルタ一理論

(2.5.30)

フィルタリングにおいて、状態方程式が非線形となる場合は、前項で記述したカルマンフィルターをそのまま使用することは出来ない。これに対応するフィルターとして、拡張カルマンフィルター・一次フィルター・二次フィルタ一等が開

発されている(片山(1983)22) ) 。本項では、一次フィルター及び二次フィルターの理論をレビューする。

( 1

_{)記号の説明}

前項で定義しなかった記号の説明を行う。

(

(t )

)は状態方程式のヤコビアン行列、 lFi(

î (t )

) は状態方程式のヘジアン行

列及びμ(t)はバイアスコレクションである。

( 2

)状態方程式・観測方程式

カルマンの予測理論における状態方程式はそれぞれ次のように表すこ

とができる。

主(t)= f (x (t) )+ &_t) ; X(t) = Xo y(t)= M(x(t)) +ε(t)

(2.5.31) (2.5.32) ここに、 &，ε は誤差ベクトル (白色性正規雑音) である。ずらし時間

を τ とすると、それらの自己相関関数はそれぞれ次式となるo

E[δ(t)δ(川

E [ ε(t )ε(t+'t') ] ^{= R(t} ^).L1 ^{(t )}

ここに、 L1(t)は D i ^{r a}c のデルタ関数である。

(2.5.33) (2.5.34)

降水レーダを用いた水文現象の予測手法に関する研究二次の項ま状態方程式をZのまわりでテイラー展開し、

χ(t) =ミ(t)+ヌ(t)より、

(2.5.35) で取る。

f (x(t ) )王f(五(t )}+A (.五(t )).e(t ) + t E ^{仰(t )}

^.lF

^(王山

ハU ・・・ハU

、，、�

^{I �}

'- '- V '-- 、

i番目

ハU ・・・ハU

中i=

である。

初一一仇祈一〈均二ι 一へd

弘一

仏

ヤコビアン行列

Ò!2 ÒXn

弘一仏側、，r h 『句B・・ザ一〈 x t 一円d

Ò!n ÒXn òfn ò!n

ÒX1 ÒX2 A

(i(t)

降水レーダを用いた水文現象の予測手法に関する研究

ヘジアン行列

Fi (.支(。

d2fl d2fl dXl dXl d _支

d2f2 d2f2 dX2d支1 己主2

a2ん d2fn dXnd支1 dXndX2

( 3

)時間更新アルゴリズム

状態、方程式が線形であれば次式が言える。

五(t) = ^f (X(t )) ;支(tk) =Xk

d2fl dXl d丸

d2f2 dX2dXn

d2ん a2f

(2.5.36)

しかし実際はf(*)が非線形であるとすると、状態方程式をテイラー展開した時に、二次の項の統計亘

にバイアスが生じるのでそれを補正する項(バイアスコレクション)μ(1)を導入する。

五(t)= ^f (X(t )) +μ(t) ;支(tk)= ^x (k\k) ; tk壬t < _tk+l e(tk) = e(k/k) =

(k/k)- X(tk) = e(k/k) - Xk _=吉(k/k)

(2.5.2.1)、(2.5.2.7)式より

ë(t)= i(t)- X(t ) = ^f (X(t )) +δ(t ) ^{- f} (れ) -μ(t )

L式に(2.5.35)式を代入する。

(2.5.37) (2.5.38)

(2.5.39)

è(t)=

(.仰(I)+t pieT(f)Fi(れ)峨t)ーん;μ t <恥1

(2.5.40) μ(1)を決定するために、独立変数の平均値を用いる。-;(k) がx (k)の不偏

推定値であるとすると次式が言える

E [e(tk) ] ₌

₀ ^(2.5.41)

;-(t)もx (t)の不偏推定値である必要があるので、次式が言える。

E[e(t)J=O E[ë(t)]=O ; tk

^�

t < _tk+l

(2.5.42)

降水レーダを用いた水文現象の予測手法に関する研究

(2.5.39)式より、

é (t)+μ(t ) = f(x (t )トf(五(t)) +δ(t )

両辺の期待値を取ると

山+附]= Eケ(x(t ) ⁾ ^)-1 (x (t川巾) ]1 (t) = E [t(x (t )+e (t ^{))-1 (i} (t)) ^]

(2.5.43)

(2.5.44)

(2.5.45)

(2.5.44)式を計算することが出来れば μ(t)は正確に算定できるが、一般的に

E[t(主(t)+e (t ))Jげ(E[x (t )+e (t ^)J)

であるので、 (2.5.40)，(2.5.42)，(2.5.45)式より近似する必要がある。

(2.5.40)式において、両辺の期待値を取ると

A(れ). I山叶ゆieT (t )Fi仏)+ 1E[δ(t )J|五(t)= 1山)]1 μ(t)= rエやiE [eT (t )Fi (Xk)e (t)J

、，、� . _

'- '- h ...

=;エゆi&tr{Fi(Xk)e(t ^)内) }J

=25仰{Fi(xk)&e (t )eT(t) J}

=宅仲{Fi ^(ik)pο

(t ) = E [e (t )e T(t )J

式を I について微分する。

P (t)= E ^[é (t ) e T(t ) + e (t ^)é T (t)J

(2.5.50)式に(2.5.40)式を代入する。

A(五(t)片(t)eT(t)

(t)= E + [{中ieT(t )附州t)い(t) + e (t )A市(t))eT(t) 五れん川(οω4小fけ1)]

+ e (t {{t �仰川

(2.5.46)

(2.5.47)

(2.5.48)

(2.5.49)

(2.5.50)

降水レーダを用いた水文現象の予測手法に関する研究

A収(t)片(t) eT(t) +δ(t )eT(t)

+ �IエØitr{lFï(五(t )){ e (t ) e T (t ) - P (t )日leT(t)

=E ^I ^ーL i=

¹ ^_j

+ e (t ) A T (五(t))eT(t) + e (t) δι(t)

+ e(t). t [ ^吉川i川(t)e T(t ) - P (t) )]1

(2.5.51)

。る〉えニ=ロカ式次シ」るす定 o JL ニ一シ」 1j

φ'&

、っ UK 従 1lp」

布 ο 分 P 規け正

UK

ドキュメント内降水レーダを用いた水文現象の予測手法に関する研究 (ページ 34-44)

= rり

5 宣告第2

，、、

お

y

五(k+ 1/ k) =a( 1 ，k)y( 1 )+a(2，k)

(2)+・ ・+a(k，k)y(k)

デ

y(k+ l/k)=M(k+ 1 ).Jr7._k+ l/k)

y(k+ 1 )=M(k+ 1 ).x"'{k + l)+y(k+ l/k)

y(k+ 1 )=y(k+ 1)手(k+ 1)

五(k+ 1/ k+ 1) =x(k+ 1/ k)+Jf(_沢k+ 1)

x(k+ l/k)=x(k+ 1)ー支(k+l/k)

x(k+ l)=x(k+ l/k)+x(k+ l/k)

y(k+ 1) =M(k+ 1 ).x(k+ 1)+ε(k+ 1)

y(k+ l)=M(k+ l).x(k+ 1 )+M(k+ 1 ).x1_k+ 1)+ε(k+ 1)

y(k+ l)=M(k+ 1 )-x1_k+ 1)+ε(k+ 1)

x(k+ 1/ k+ 1) =x(k+ 1 )x(k+ 1/ k+ 1)

=x(k+ 1 )- [ えれl/k)+lK長k+ 1) ]

=x(k+ l/k)-lKy(k+ 1)

(双k+l/k)-!Kyてk+1))上知+ 1)

E[x(k+l/k).yT(k+1)]-!K E [Y(k+1).yT(k+1)]=O

抑制/k). (MT(k+l/k).xT(k+l)+εT(k+1))]

={(M(k+ 1 }xTk+ 1) +ε(k+l)). (MT(k+l).xT(k+l)+εT(k+1))]

[玄(k+l/k)王T(k+ 1) .MT(k+川 (2.5凶)

= [M(k+l)沢k+ 1)玉川+l)).MT(k+l)+ε(k+ 1)εT(k+l)]

I.x(k+ l/k).;?(k+ 1) I

P(k+ l/k)三E[玄(k+1)玄T(k+ 1)]

R(炉開=E[e(k+ 1)εT(k+ 1)]

!K(k+ 1 )=tP(k+ l/k)MT(k+ l)J. [M(k+ l)P(k+ 1/ k)MT (k+ l)+R(k+ 1 )]-1

正'(k+ 1/ k)=x(k+ 1 )-x(k+ 1/ k)

双k/k)=x(k)-x(k/ k)

x(k+ 1)=φ(k+ l/k).x(k)+δ(k+ 1)

五(k+ l/k)=�k+ 1/k}x1_k/ k)

φ .xてk/k)=φ.x(k)ーφ.x{_k/k)

=(x(k+

)-�(k+ 1))仰+ l/k)

:. i(k+ 1/ k)=φ.xT.k/k)+δ(k+ 1)

P(k+附=E[双k+ l/k).xT.k+ 1 / k)]

=E [x(k+ 1/ k) 双が同T]φl +E [ δ(か1)δ(k+ 1) T ]

:. P(k+ l/k)=φ.P(k/ k) .φ+Q T

Q= & 伽1)δ(k+ 1) T ]

五(k+1/k+1)=支(k+ l/k)-lK只k+1)

五(k+l/k+ l).llKj(k+ 1)

(双k+1/ k)-lK只k+1 )XK<只k+1 ))T =0

P(k+1/k+1)=E [双k+1/k+1)気k+ l/k+ 1)]

=再(双k+1/ k)-lK沢k+1)阪k+1/ k)-lK沢k+

)y]

=E l(双k+ l/k)-lK只k+ 1) )X(k+ 1/ k)TJ +E[(沢k+ l/k)lK沢k+1))XlK，ヌk+

)T]

=E[双k+ l/k)玖k+ l/k)TJ -KME[玖k+ l/k)双k+l/k)TJ

=P(k+l/k) -KMJP(k+ l/k)

=[J/-K MJ!P(k+l/k)

第2. 5.

( 1

(

(t )

î (t )

( 2

主(t)= f (x (t) )+ &_t) ; X(t) = Xo y(t)= M(x(t)) +ε(t)

E[δ(t)δ(川

E [ ε(t )ε(t+'t') ] = R(t ).L1 (t )

降水レーダを用いた水文現象の予測手法に関する研究 二次の項ま 状態方程式をZのまわりでテイラー展開し、

χ(t) =ミ(t)+ヌ(t)より、

f (x(t ) )王f(五(t )}+A (.五(t )).e(t ) + t E 仰(t )

(王山

ハU ・・・ ハU

、， 、�

ハU ・・・ ハU

初一一仇 祈 一 〈 均 二ι 一 へd

弘一

仏

弘一仏側 、，r h 『 句B・・ ザ 一 〈 x t 一 円d

(i(t)

Fi (.支(。

d2fl d2fl dXl dXl d 支

d2f2 d2f2 dX2d支1 己主2

a2ん d2fn dXnd支1 dXndX2

( 3

(2)+・・+a(k，k)y(k)

E [ ε(t )ε(t+'t') ] ^{= R(t} ^).L1 ^{(t )}

降水レーダを用いた水文現象の予測手法に関する研究二次の項ま状態方程式をZのまわりでテイラー展開し、

f (x(t ) )王f(五(t )}+A (.五(t )).e(t ) + t E ^{仰(t )}

^(王山

ハU ・・・ハU

、，、�

ハU ・・・ハU

初一一仇祈一〈均二ι 一へd

弘一仏側、，r h 『句B・・ザ一〈 x t 一円d

d2fl d2fl dXl dXl d _支

五(t) = ^f (X(t )) ;支(tk) =Xk

五(t)= ^f (X(t )) +μ(t) ;支(tk)= ^x (k\k) ; tk壬t < _tk+l e(tk) = e(k/k) =

(k/k)- X(tk) = e(k/k) - Xk _=吉(k/k)

ë(t)= i(t)- X(t ) = ^f (X(t )) +δ(t ) ^{- f} (れ) -μ(t )

E [e(tk) ] ₌

t < _tk+l

山+附]= Eケ(x(t ) ⁾ ^)-1 (x (t川巾) ]1 (t) = E [t(x (t )+e (t ^{))-1 (i} (t)) ^]

E[t(主(t)+e (t ))Jげ(E[x (t )+e (t ^)J)

=;エゆi&tr{Fi(Xk)e(t ^)内) }J

=宅仲{Fi ^(ik)pο

P (t)= E ^[é (t ) e T(t ) + e (t ^)é T (t)J

=E ^I ^ーL i=

+ e(t). t [ ^吉川i川(t)e T(t ) - P (t) )]1

。る〉えニ=ロカ式次シ」るす定 o JL ニ一シ」 1j

布 ο 分 P 規け正