直交射影法と階段状分布近似による 　　　　　　　　離散形非線形フィルタ

直交射影法と階段状分布近似による

        離散形非線形フィルタ

（昭和50年5月23日 原舗 受理）

情報工学教室  高    田       等

ADiscrete Nonlinear Filter based on the Estimation of the    Orthogonal Projection］M［ethod and the Distribution      approximated by Staircase−formed Functions

bv Hitoshi TAKATA

  In this paper， a sequential ロonline且r I］1ter is presented ［or lhe discrete noisy llonlinear syste皿s． The featllres of the丘1ter are as follows． 1）The optimal estimation formula， which is derived by orthogonaUy projecting on the linear manifold generated by polynomials of inno−

vations， is modi丘ed such that thc probabiiity density o［it is normed andロonnegative．2）In order to have a sequential田tering algorithm， any probability density of state prediction is propeτ吋constructed as a staircase−formed funotion∫rom the available prediotion moments．

  Digital simulatiOnS Show that this創ter iS SUperiOr tO lhe Ordinary filters in acouracy for the systems wlth higl・n・nlinearities，

      できなくなった場合のみ，階段状分布関数の存在韻域と

t序 論        の船い柵靴の補正を行なう方法について撚し

一般に非線形フィルタ柵成には，ベィズ法の条件付確   た。さらにこうして補正された推定式を基にし1非正規 率密度閏数による立場から多く論ぜられている1川。し  分布の立場から1次と2次のみならず高次モーメントま かしながら，ペィズ定理にとらわれることなく，確率変   でをも考血したフィルタアルゴリズムについて考察し 数の直交性に基づいて非線形フィルタを排成する立場に   た。従って本フィルタにおいては，偶数次モーメントの たつ方法も皿要である。これに属するもののうち1次元   負値現象などは生じない。

連続系に対してはBalakrishnan3）の方法，および班者   最後に，本フィルタを状態値について2次の非線形系 ら4］の方法がある。後者の方法は，Balakri曲nanの方   に適用した場合のディジタルシミュレーション結果を示 法の内の多項式形推定に相当するものを，多次元離散値   す。これにより非線形性の強い系において 木フィルタ 系にまで拡張して，観抑］情報の多項式からなる直交基底   の有効性が認められた。

変数鱒入と誼交射影麟の適用1こより軸されたも 2問題の設定

のである。得られた推定基礎式は，観測情報のべき級数

により衷わされた逐次計算可能なものである。さて，こ   力学系および観測系がつぎの非線形差分方程式で3己述 の基礎式は無限項表示であるため，実際の物理系への適   される系について考察する。

用酬・舗限項で近似される必要がある・とくにぷ  力学系エ、，、−f（副やκ  （・）

㌶㌶㌔ζ苧ll㌫㌶：繊：竺  蹴・」系炉・（工π）÷・・  ②

メントの負値現象を避けるための確串性の配慮も必要で   ただしエκは〃次元状態ペクトル・〃κはm次元観測 ある5〕。      ベクトル，f（工κ）とρ（エκ）はエκの多項式形ベクトル  そこで本稿においては，まず任意の分布が十分に近似   箇関数．賑と苗κは互いに独立でそれぞれ平均値が できる階段状分布を用いて，状態予測分布を近似的に表   零，共分散がC・・曲（非負定値行列）とC「・・面（正定値 現することを考える。つぎに上述の有阻項近似の推定式   行列）の白色雑音である。以下・状態予測・齪測予測お に対し，これに関する偶数次モーメントの正値性が確保   よび状態推定の平均値と共分散をそれぞれ

主。嗣一五長／γ。．1），P．旧一C。。（x．佑一1） ては3および櫛の結果を凧・てフ・ルタアル訓ズム       （3）    を導出しよう。

   」加口・＝E（〃κノγκ＿D， Sκぽ＿ユ＝C四（〃κノγκ一t）    3．確率性を具備した推定式       （4）

       最適推定式（8）は無限項表示であるため，具体的な系

 ハ      ロ

 エκ1κ＝E（エκ／γκ） Pκぽ＝Cρ「（工五／γ口  （o）   への適用の際には有限項で近似される必要がある。この と記す。なお（1），（2）式において，もしf（喬），ρ（エD   場合有限項打ち切りにより偶数次モーメントなどが負値 が一般の微分可f屯な非線形関数の場合には，それぞれ   を示す現象が起りフィルタ発散の原囲につながる。そこ 三五1κ，三κ1κ一、のまわりでテーラー展開し，有限次まで   でこれを補憤するため（8）式において任意の次数lvまで 考慮した多項式関数で近似することにする。       ：号慮した推定式に対し・つぎのような確率性を保証する  文献め〜6）によれば，欺の非線形スカラー描閲数   重み係数λκ（0く㌔≦1）

占伝κ）に対する最適推定式は．観測情報       鬼…＝』刑κ一1十圭』C｛κ、，）ξ｛π．，，    （14）

       ノ ユ

  γκ全σ〔酬，〃ガ・㍉〃訂      （6）   を導入する。ただし（14）式，すなわちλκは具体的に が与えられたときに，推定誤差の2乗平均        は以下のように定められる。

      A       声     いまρ（エκ／γ五。1）が階段状分布として与えられてい

  Eぽ［（力（エκ）−A・1・班・］  の るとしよう．ここで観酬［ 。が得られたならば任意の

』舗が2乗平均収束の意味で有界なイ直をもち，さらに ξ，。．，，α、。．，コ（エ。）σ＝1β，3、…）カ｛定まる・・一・，。従・

観測情報のべき級数で表わされる場合には結果として    てイノベーション（槻測情i報）をN次まで考慮した（8）

  免。1戸￡。＿＋豆c、κ。，ξ，脳）  （8） 式は

         ヨロエ       ム       ば      ロ

粥られている．ただし条件付確率齪酬ρ（エ。／γ。）， 力・1…ぽ一・＋恩C一ξ…，  ㈹ ρ縮γ・−1）に対し       一∫占（工κ）・（エ・／γ・−1）輌

』…一∫力ω・（エ・／yDゴエ・ （9）  ＋曇［∫ゐ（エ。）α、 （エ・）・（エ・／γ・−1）

  孟一一∫力（工π）wγκ一1）ゴエ・ （1°）       x4エ・ユξ一

      高々ゴ乗のぬを含む」次元配列である・    ただし

輪≧霊；竺㌫一可測閏数の 飯三F（」：R）・㌣鋼  （・7）

c−一 逞ﾍω一・（エ。／γ。．1）輌（・・） F（エ・）＝1＋，苔α・エ・コωξ刷@      となる。（9）と（16）式の比較により（16）式㌫

として1己述される・        がρ（耐γ。）1・対応する．しかしながら一殿には（17）

さて緒においては・状態酬酬 ｱ工・ノγ・−11に 式口。は（、3）式の条件を齪しない．これが1聯

対しては 任蜘分布をト分に近似し得酬段状酬で

瓢ll蕊㌶慧㌶警警㌻麟め・ま撒状分布として与えらオ・て一（工鋼

するの端である・rでいう酬生とは・確轍の存在慧≦．＿≦免＋一≦」≦。〕

関数が本来馴せねば脳い｛生質

（i）∫ア（エ／肋一・   （12） ぎ1・β．が酬をとる場合のみF（エ・）を1β・櫛上

 （ii）すぺてのエに対し ρ（エ／γ）≧0  （13）   に平行移動させた関数

のことである．つぎの3および撒おいてはそれぞ F・ω一・塙α・指・コ（エ・）ξ一＋欄φ㈲

れ，（8）式に対し螂性を酬し雌定式，およ柵段         （ユ9）

状分布関数の構成法について考察し，さらに5節におい    ただし

⇒；：ll：1：遼：（2・）㌶㌶蕊當す㌘乱1≦∫⑳

を導入する・すなわちβ・が非負の場合1・はF・（エκ）＝    確串密度関数、ρ（エ）の性質（・2）式

鑑1才蕊；）（：㍑1：（認：袋1；：籔；三 ∫ρ（エ）ゴー・

なる．さらにまた（12）式の条件をも1菌足させるために， より

これを正規化したつぎ鵬

P琵＝Y蕊鷲㌃  であ＿原点のまわりの（＿．．＋み㌫．

  一（  ぐマ1十Σα［x．」］（エκ  」一］）ξ 、÷…β。・φ㈲）・（エ。／γ．．1）ントは麟

  一    1＋1β・1岬・）      巳・｛1・・（エ）ゴエーE、8・｛お

       ザ

  ＝（1十Σ｝λκα［x．」コ（x五）ξ（π．〆ぱ♪（エκノγκ一1）（21）    より      ノロエ

ただし

@      E騨一蕊…摯1H・．綿晶1

  拓＝1・0ノ（1÷1βκ1φ（βD）      一         （＿。。〈βば≦o，oくλ。≦ユ） （22）         x（°｛識・一α｛卿 （2°）

       である。

を導入する。ここで（工6）式のρκの代わりに（21）式

      4．2．階段状分布の構成

の確靴（12）・（13）式を齪する・鮭ft入すれぽ  いま・次から（M−、）次までの酬予測モー〃ト

確串

hピ轍     榊一G≦力＋五一五≦㍑

    一∫卉（      む喬X1十Σλκα￡担（エκ      」−1）ξ一）  が与えられているとしよう．撫は（26）式の欄の

      Xρ（エκ／γκ一1）ば叛     みから状態予謝分布ρ（欺．1／γκ）として近似的に階段

     ム        ふド

    ＝転1κ一1十Σ｝λκC｛κ、∫，ξ（エ．∫，     （23）    状分布を構成する。

         ゴ ユ

が得られる。（23）式はすなわち（14）式モのものであ   まず階段状分布の区間巾としては・任意の正の実数価 るo      「査L1に対し

      トrκ．11／E（埠＋1ω／γκ）÷五（x翻山／γκ），

4・階段状分布による端予測分布の近｛以     ，．、1〆E（報＿／γ，）÷恥＿／r。）］

 一一般に任意の分布閏数は階段状分布により十分近似で       （1≦1≦め  （27）

きる・そこで本節においては・与えられた有限次の状態   を仮定する。この間を上、（1≦∫≦の分割し，その値を

剛モー・ン陸雄し西みつき酬誤差の櫨小の  （・刷・，，。，…，・、岬）（・≦1≦ ・）

鰍から髄酬段状分布を構成Lこオ・をも・て樵 とする。

予測分布の近似とする・      つぎ嚥段状ク｝布の高さH，、臼，帽（・≦属い

       1≦1≦のについては（26）式の各種のモーメントと階        段状分布としてのモーメント（25）式との差の2乗和が       勘      最小になるように定める。この際，各差の2乗の値がそ       肱。  刊q圓      の絶対髄の大きさに関係なくほぼ同程度の評価をもつよ

      酬    うにつぎのπみつき評嚥埠；

  ㎏惜副即｝1仇。H…。網       〃一、11，蕊、．1κ・・ドJE〜躯［1川／γ・）

     図12次元階段状分撒    竜藁藻1丹・ド肺晶一

鐘変謬鷺燃㌘，…，細］．，醐了   聯1＿㈱い8）

を最小にするように選ぶ。ただし      とより，有限次（ルド次）の状態推定モーメントと雑音

脳一早＝

諱^甦拓㌘   玲懸忽）：：驚竺麗㌫

       〒幕1碧1 冨、H＝㌃、∫，＋1       （35）

       xl剛仁・紬｝2（29） 芸鷲㌫ご歴竺誓：：晋欝獲㍍2：

       ほ      ヘ      ココ

  段＝1・。μ（σ，，，、r・1，，，）      れbより42節に従し すべLの        」⇔1

      （一様分布に対応）  （30）      ｛ロ ∫ω，』ゴ〔 1 … ，・・い。ゴユ≦輌1≦上r・1≦1≦n｝

である。この場台（24）式の等式関係より（28）式の未  が求まり・結局階段状分布が定まることになる・ただ 知数∫」亡□の個数は（L1十L叶…÷L一1）個である。   し・任意に与えられた㌔・1・（Lβコ≦ ≦肩に対して また（田）式はH国に閏し2次式であるので，容易に   得られた階段状分布が・もし（13）式の確率性の条件を

  ・P・／∂恥一・     ㈹ 

^齪してし 鮒れば・耐一｛上β1≦「≦H｝の値

が解けて諜知数岳・についての解は（L1橘÷ ＋ を見咄すまで4．2胸操作を鋤返す．註ここで エ・−1）元避欺方程式とし碇まる・こうして勅 一般には真の澱確率密度関数〆エ。／γ。）は，時刻κ る階段状分布耀られる・      の増加ととも｛・・関雛近づいていくことに髄して，

5フィルタ構成         次時刻の観測1直〃阻が得られた後に（33）式のコ紛散

       に対し，そのノルムがもし  本節においては3，4節の結果を用いてフィルタを構

成しよう。      °＜llP−＋lll≦｝lp・1・｝1   （36）

 5．1， いま観測値〃κと，状態予測分布P（エx／γ・π一1）  を満たさなければ階段状分布ρ（エκ・】／γR）に対し・λ『・

が領域［ロ川，≦エκ【P≦屯r＋ユf」日1≦1≦η〕をム（1≦ …h仏；1≦ ≦・）の値を変えて4・2・節の操作を

∫≦め分割した階段状分布により与えられているとしよ   やり直すフィードバック機構を設ける。これにより・

う。       （36）式の条件をも満足させる。

 5，2． ここでイノペーシ目ンの次数1Vを定めれば，    5・4・先験情報として先験分布あるいは有限次の先験 すべてのα［κ」ユ（エκ），ξ｛凡，）（ゴ＝1，2，…，jv）が求まり   モーメントが与えられれば・5・1・〜5・3・節の繰り返しに 結局（18）式のF（エ。）力淀まったことになる。従。て よって逐次的に推定値力淀まる・このアルゴリズムを以 階段状分布の存在阻域に対しF伝κ）の最小倣βκ，す   下に示す。なお本推定ブロック図を図2に示す。

なわち〔22）武の娠が定まる。つぎに任意の正の整数

M卓に対し，1次からM牢次までの各種の推定モーメ         〈推定アルゴリズム〉

ントが（14）式から      ［0］与えられた先験情報に対する階段状分布♪（エ】ノ

カω一麺（・≦品…÷み≦鵬 田γ㌶驚布卿兵．Dおよ酬伽．が与

       （32）       えられる

とおくことにより求まる欄。この際，確率性の補正係   ｛2］ イノベーションの次数Nの選定

数λκの導入により．特に，推定共分倣       ［3］［1L［2］と観測方程式よりF（エπ）の導出

  Pκ1苦＝E【（     L工κ一工κ1κ）伝κ＿三泊κ）T／γκ］（33）    14］ βκすなわち㌔の決定

      ［51 Pκlxを導出し，0＜lIP川副≦lIP席一］1κ一111が満た

は正定齢列を示しさらに酬め偶数次モー・ント  され鮒れば塒刻前の17はで帰り［7］の値を

       ム

  E［（工古（∫）一工田κα））2ゲγκ1      選定し直す

         （f＝1，2，・‥；ユ≦∫≦η）（34）  ［6］三。、。．E（誌・。1ノγ。）（2≦計∫、＋…愉≦

      ユロエ は決して負倣を示すことはない。      Mつの決定

 5．3．つぎに逐次推定計算を可能にするため，次時刻   口1 兀・bM・｛」L日コ≦「≦π］の値の選定

の状解測分砲工。。ノγ。）を勅ておく。   18川61と蝉紡程式よりE〜掃・1…／r・）

 さて（1）式の充XK）がエκの多項式関数であるこ    （1≦五十∫2十…十み≦」1イー1）の導出

先験虹〔 @ Y、・9（x、｝・W、  観測値

階段状分布

P（X 汲倹痰g〕

z一l

製測予測モーメソト

Y綱・s

 最適推定式

倉一↑．鍵c．ξ

kIk   klk−1 j＝l  kj kj

P〔X^{求D11写、）}

H 1

H2

■●●●

al a2  a3白・・．  aL＋1

Y巴s

確率性の条件 No

β・㌔の決定

 硫串性を貝備した推定式         N

ハ     ハ

h ：＝h  ÷』 ∫ O ・ξ

klk     k止一l     k j＝1   kj    kj

階段状分布の高さ H

司一・1（1㌃㌫迎立1

ただし   M−1

PI一 bjj〔E（・㍊1亨、）

一⊥』H〔。jrL、」←｝〕2  」＋1il i  i41 1

階段状分布の区冊巾         ム_{・1−一「}

ﾎA．11、

  三

、 ＝r，一＋棄

 L＋1     k十llk  ｝［＋llk

Z     O〈P ≦：P

一1No

klk− k−11k−1

Yes

状態推定モーメント  A      i

x、1、・E（x、1宥k）

    〔 i＝ 21 3， 4， 1 ・ ．， M  ）

M，1、．γの選定

XH＝f（Xk）＋Vk

状態予測モーメソト

   lE｛x計11㌢、）（、＝1隅．，M−D

      図2  推定過程のプロック図（スカラー表示）

［9］｛α∫r｛」ハ；1≦ら≦」し」十1，1≦ ≦吋の選定        力学系と観測系がモれぞれ

［101 〔、甘〔h．．．．∫1．＿．，。コ：1≦ ≦L，，1≦1≦n｝の導出        エκq＝エR−0．1鳩十む      （37）

m］ 19L ［10］を基にして構成される階段状関数が，     」・圧＝β㌔十鳩十芒ψ       （38）

 確率触条件を満たさなければ［了］へ帰り・田の の・次元の例題に引、て撚する．初靱幽ま。F1．0

 値を選定し直す        で蠕はそれぞれ

1121畷状分布・』／mの構成    。〜N（，，。，、。一・），1。〜N（、，，、。，。．25）

11310蹴がκ＋・へ魅恥返し綻計算もエ11 な珀色正蹴音とした．ただしN（ロ：占，c）は螂変

 へ進む・         数。が正規分布をしその平均｛血胸で分麟、である

杜例 題         ことを示す。また先験値としては

      x1〜1V（泊10，2、0）

       ロ 

管⊥      R、，、

⊥｝       °

c6

      β＝O．5       。6

c

       口ロロユ

1＿ロ       。

コ⊂……    ．   。

：：1  昂 ㌦岬  1：1

：II −・…一   一        ・

．。！       一叫       ・110

5   埠   15   ・o       コ0

       5t叩

巾1 2nd Ord■r川1■r

図3表1の・L＝4でβニ0．5の        図4表1のL＝4でβ＝Lおよび

   ときの推定軌跡       L＝4でβ＝5のときの推定軌跡

P｛X・∫暫、．1）        O．5

APriori

t

｜ （K＝1）

1 1

一占  一3   −2

^一1

1      α5

1 l

I

→

（K＝2）

 仁       一4   −3   −2   −1    0    1    2    3    4

0

巴 器

q5

¹

〔K＝3）

1

＝4  −3 一2

一1

｜

t 2 3 4

．5

（K＝4）

X2

x，

一』L   −3    −2    −1    0     1     2     3    4

  （K＝5）

｛     1刀

  ．1。 1 2X・．1。

図5 図3において状態予測分布として用いられた階段状分布関数（破線は状態真値斯）

（K三6） （K＝10）

↓

（K＝15）

1．O 旧0 ¹⁰⁰     ｝oo

o．5 _5， 1 50 ^5．o

1 o

     X 1   2

o ^1X可 0 0

表1βをパうメjタにしたときの推趨差の 熾，本手融L−4で約・57m，（ミリ秒），L−5で      2乗P」＝君』一』）2砒較  系勺・96m・および2次近継フ・ルタ噺勺4m、であっ

＼〃1・。d．。，d。， S亡a 「cas鑑h・g・n・1 た゜

  β ・     田ter̲1 L−・1

10 0．0423 l I．357 _0．7097 5 0．1936 1 1．46フ _o．B493

1 芸 l  l．647 _1，544

0．5

x

o．1 x

_X

0．01

x

苦：Very large value       ・ なりの

L＝5     7・結   言

0．7097     非線形離散形系に対し．直交射影法によるイノペ＿シ 1．544     ヨンの多項式形推定を基にした逐次計算可能なフィルタ 且336     アルゴリズムを導出した。その際状態予測分布忙ついて 3．249     は・階段状関数による近似法を提案し確率性をも配慮し        た。これにより非線形性の強い系に対しても か        精度で推定が可能になった。本フィルタにおいては，イ

を選び，評価としては3・鋼融推趨呈差の2乗  ノベーシヨンの次数（め描段状1蝦の幽（L）・およ

〃一

を考えた・      最後iこ舶頃御欄を贈っている九州大学工学部辻節

 表1は観測線形項の係数βをパラメータにしたとき   三教授に心から謝意を表します。

の本フィルタ（Staircase−Orthogona1田ter）と2次近

似法フ・・レタ（2・d−0・d・・丘1…）・・1、よる（39）式P∫    参考文献

の比繊である・ただし本フ・ルタに担・ては囎状分 1）例えば埴穂髄システム制臓。。ナ社 紬醐巾は等分割され，標酬として凄内の（L＿  （1972）．

璽：：：議繧竃糠i藷1璽璽ll欝聯憲1憲

において上＝4醐合のそれぞれβ一〇．5およびβ＝1 3）A・T・B・1・k・i・h・…AG・・e・alTh…y。f

：鷲㌶㌫㌃（隠㌫麟：：藍  i欝：1「工慧ti蒜蠕霊；諸膿

 これらの結果よりβがかなり小さい場合において    め 辻・高田：直交射彫法による離散形非綜形フィル も・すなわち非線形性が極めて強い場合においても，本     タ1計測自動制鋤学会論文集・8−】・1／10（1972）

フ価タはかなり繍度幽態値を綻することができ 5）辻 語竺騨性鞭を伴う馴剛次モー・

た・さらに図5より離状分布（状態予訓扮布）剛 ；3慧㌶，；；；ヂ鵬動制…論文集・

刻の剛とともに状態真値を・1・心1・したδ関数へ近づ 6）且T・k・t・・Di・c・，t。 N。nli。。a， E、tim。ti。n いていくことが実際に確認された。      with Higher Order Innovations by using

F鑑㌫瓢は㌫㌻：耀籔㍊：ぎ：膿lfs搬