3. 1.
緒⾔
本章では，動的環境推定モデルの理論的構築を⾏う．まず，本研究が⽬的とする動的環境モデルの概 念を図 3.1に⽰す．このモデルでは，軌道⼒学に基づいた数値シミュレーションによる推定環境の更新 と，観測データと推定環境を⽐較することによる修正を逐次的に繰り返すことで，常に最新かつ観測デ ータに対し最適な環境推定結果を得ることを⽬指す．
図 3.1 動的環境推定モデルの概念図
上 記の ように シミ ュレ ーショ ンと 観測デ ータ を融 合さ せる概 念は ⼀般 に「デ ータ 同化 (Data
Assimilation)」と呼ばれる．そのため，本研究においても3.2節で述べるデータ同化の⼿法に基づきモデ
ルの構築を⾏う．データ同化の実現には複数の⼿法があるが，それらのうち本研究では「逐次モンテカ ルロフィルタ（Sequential Monte Carlo Filter, SMC Filter）」を採⽤する．SMCフィルタを採⽤する理由お よびそのアルゴリズムについては，3.3節で述べる．そして3.4節において，第2章で得られた観測され 得る微⼩デブリの軌道の特性を活⽤し本研究で構築する動的環境推定モデルの具体的なモデリングを
⾏う．
3. 2.
データ同化
樋⼝ら [29]によれば，データ同化とは『シミュレーション科学のような演繹的な推論と，統計科学に 代表される帰納的な推論を融合するためのプラットフォーム』であるとされる．データ同化は計算機の 発達に伴うシミュレーション技術の向上と，⼈⼯衛星の利⽤や情報化の進展による⼤量データの蓄積を 背景として，1990年代より気象学・海洋学の分野を中⼼に発達してきた．前節でも述べたように，シミ ュレーションによる予測とデータによる修正を⾏う本研究の環境推定モデルにおいてもデータ同化の
システムモデル
まず，あるシステムのある時刻tにおける状態変数をまとめたものを𝒙…= (𝜉l, 𝜉u, … . 𝜉‰)Šとおき，これ を状態ベクトルと呼ぶ．状態ベクトル𝒙…の時間発展を計算機の上で数値計算するためには離散化が必要 である．ある時刻t -1からtへの離散化された更新操作は式（3.1）のように書ける．ここで， fは，𝒙…ol
から𝒙…への写像である．
𝒙…= 𝒇…(𝒙…ol) (3.1)
しかし，離散化による数値的問題や，全ての状態変数を知り得ないことから，現実世界の真のダイナミ クスを上記のような𝒇…によって計算機上で完璧に再現することは不可能である．すなわち，式（3.1）は 完全な等号ではなく近似等号であり，不確実性を含む式（3.2）となる．
𝒙…≃ 𝒇…(𝒙…ol) (3.2)
そこで，真のダイナミクスとシミュレーションモデルの乖離に対する緩衝材として，データ同化ではモ デルに確率変数を導⼊することでこの近似等号を表現する．この確率変数をシステムノイズと呼び，以 後𝒗…で表す．システムノイズを導⼊したシミュレーションは式（3.3）となり，この式はシステムモデル と呼ばれる．
𝒙
…= 𝒇
…(𝒙
…ol, 𝒗
…)
(3.3)システムモデルは確率変数を含むため，その計算結果も決定論的なものではなく複数の解を持つ確率 論的なものとなる．すなわち，時刻 t の状態ベクトル𝒙…は前時刻の状態ベクトル𝒙…olが所与のときの条 件付き確率として𝑝(𝒙…|𝒙…ol)と表される．
観測モデル
システムの状態ベクトル𝒙…に対して観測データ𝒚…が得られるとき，この関係を𝒚…= 𝒉…(𝒙…) と表す．𝒉… は観測演算⼦と呼ばれる．システムモデルの議論と同様に観測においても不確実性を想定すれば状態ベ クトルと観測データの関係は式（3.4）のように書ける．
𝒚… ≃ 𝒉…(𝒙…) (3.4)
そこでシステムモデルと同様に近似等号を表現する確率変数𝒘…を導⼊し式（3.5）のように表す．
𝒚… = 𝒉…(𝒙…, 𝒘…) (3.5)
この確率変数𝒘…を観測ノイズ，式（3.5）を観測モデルと呼ぶ．なお，観測ノイズは観測機器の特性に
ルの差をも含んだものであることに注意する必要がある．
状態空間モデル
データ同化では，システムモデルと観測モデルを連⽴させたモデルを考え，これを状態空間モデルと 呼ぶ．
𝒙…= 𝒇…(𝒙…ol, 𝒗…) [システムモデル] (3.3)
𝒚…= 𝒉…(𝒙…, 𝒘…) [観測モデル] (3.5)
現実世界における様々な応⽤を考えたとき𝒇…，𝒉…は⼀般に線形ではなく，また𝒗…，𝒘…もガウス分布で あるとは限らないため，式（3.3），（3.5）に⽰す状態空間モデルは⾮線形・⾮ガウス状態空間モデルと 呼ばれる．
上記の状態空間モデルを式（3.6），（3.7）のように条件付き分布を⽤いて表現することで，さらに⼀
般化された状態空間モデルを考えることができる．
𝒙… ~ 𝑝(𝒙…|𝒙…ol) [システムモデル] (3.6)
𝒚… ~ 𝑝(𝒚…|𝒙…) [観測モデル] (3.7)
システムモデルは，時刻t における状態ベクトル𝒙…の分布が前時刻の状態ベクトル𝒙…olが所与のとき の条件付き分布となることを意味する．また，観測モデルは時刻tにおける観測データ𝒚…の分布がその ときの状態ベクトル𝒙…が所与のときの条件付き分布となることを意味する．このように表現された状態 空間モデルを⼀般状態空間モデルと呼ぶ．
逐次ベイズフィルタ
⼀般状態空間モデルを考えることで，状態ベクトルの分布に関して以下のような漸化式[30]を⽤いる ことで推定することができる．なお，本推定の基礎となるベイズ推定の理論については付録Aに述べる．
ある時刻t – 1において，状態ベクトルの確率分布𝑝(𝒙…ol|𝒚l:…ol)が得られているとする．ここで𝒚l:…olは
時刻1からt – 1までの全ての観測データを意味し，𝑝(𝒙…ol|𝒚l:…ol)はそれまでの観測データに基づいて推
定された𝒙…olの分布を表している．
まず，𝑝(𝒙…ol|𝒚l:…ol)から漸化式（3.8）を⽤いて⼀期先の予測を⾏うことができる．
𝑝“𝒙…|
𝒚
1:…ol” = • 𝑝“𝒙…, 𝒙…−1|𝒚
1:…ol”𝑑𝒙…ol= •𝑝“𝒙𝑡|𝒙𝑡−1,
𝒚
1:𝑡−1”𝑝“𝒙𝑡−1|𝒚
1:𝑡−1”𝑑𝒙𝑡−1 = •𝑝(𝒙𝑡|𝒙𝑡−1)𝑝“𝒙𝑡−1|𝒚
1:𝑡−1”𝑑𝒙𝑡−1(3.8)
𝑝(𝒙…|𝒙…ol)は式（3.6）に⽰すシステムモデルであるから，⼀期前の確率分布とシステムモデルによって 現時刻の状態ベクトルの確率分布が推定できることがわかる．このようにして計算された前時刻までの データに基づく推定分布を，予測分布と呼ぶ．
次に，予測分布を計算した状態で最新のデータ𝒚…が得られると，ベイズの定理を使って漸化式（3.9）
により予測分布が修正される．この操作をフィルタリングと呼ぶ．
𝑝“𝒙…|
𝒚
1:𝑡” = 𝑝“𝒙…|𝒚l:…−1,𝒚
𝑡”=𝑝“𝒙…, 𝒚l:…−1|
𝒚
𝑡” 𝑝“𝒚
𝑡|𝒚l:…−1”=𝑝“
𝒚
𝑡|𝒙…, 𝒚l:…−1”𝑝(𝒙𝑡|𝒚l:…−1) 𝑝“𝒚
𝑡|𝒚l:…−1”=𝑝“
𝒚
𝑡|𝒙…”𝑝(𝒙𝑡|𝒚l:…−1) 𝑝“𝒚
𝑡|𝒚l:…−1”=𝑝“
𝒚
𝑡|𝒙…”𝑝(𝒙𝑡|𝒚l:…−1)∫ 𝑝“𝒚…, 𝒙…|
𝒚
1:…ol”𝑑𝒙…= 𝑝“
𝒚
𝑡|𝒙…”𝑝(𝒙𝑡|𝒚l:…−1)∫ 𝑝(𝒚…| 𝒙…)𝑝(𝒙𝑡|𝒚l:…−1)𝑑𝒙…
(3.9)
𝑝(𝒚…|𝒙…)は式（3.7）に⽰す観測モデルであるから，予測分布と観測モデルを⽤いることで最新の観測デ
ータに基づいて修正された状態ベクトルの分布を推定できることがわかる．この修正された状態ベクト ルを，フィルタ分布と呼ぶ．なお，𝑝(𝒚…|𝒙…)は状態ベクトル𝒙…が所与のときに観測データ𝒚…が得られる 確率であり，得られた観測データのもっともらしさを表しているため，尤度とも呼ばれる．
このようにして⼀期先予測とフィルタリングを繰り返すことにより，逐次的に最新の観測データに基 づくフィルタ分布を得ることが可能となる．以上が逐次ベイズフィルタに基づく状態ベクトル推定の原 理である．
3. 3.
逐次モンテカルロフィルタ
3.2 節の理論に基づくデータ同化の⼿法としては，カルマンフィルタおよびその応⽤形である拡張カ ルマンフィルタやアンサンブルカルマンフィルタが様々な分野で⽤いられている．しかし，それらの⼿
法はシステムの状態と観測の間に線形の関係が成り⽴つことを前提とする．本研究で対象とする問題で はシステムの状態であるデブリの分布と観測データである衝突位置の間に線形の関係を規定する事は できないため，カルマンフィルタを利⽤することはできない．
そのような⾮線形の時系列モデルにも適⽤可能なデータ同化⼿法としては，逐次モンテカルロフィル タ(Sequential Monte Carlo filter, SMC)が提案されている[31]．後述するようにSMCフィルタでは観測モデ ルとして尤度𝑝(𝒚…|𝒙…)の計算⽅法が与えられていればよく，デブリと衛星の衝突の期待値が衝突フラッ クスとして確率論的に記述されているデブリ環境のモデリングとの親和性は⾼いと考えられるため，本
研究の環境推定モデルではSMCフィルタを採⽤する．なお， SMCフィルタは粒⼦フィルタ（Particle filter, PF）と呼ばれることの⽅が⼀般的であるが，本論⽂ではデブリ粒⼦（Particle）との混同を避けるため，
SMCフィルタの名称を⽤いる．以降にSMCフィルタの基本的概念とアルゴリズムについて⽂献[32]に 基づき詳述する．
確率分布のモンテカルロ近似
データ同化を計算機上で実現するには，条件付き確率分布をどのように表現するかという問題が存在 する．いま，⼀次元の変数xに対する真の確率分布p(x)が図 3.2上段のような形で与えられているとす る．図 3.2中・下段は，真の分布を近似して計算機上で表現する様々な⼿法について概念的に⽰したも のである．
x p(x)
x p(x)
x p(x)
x p(x)
x
p(x)
計算上で確率分布を取り扱うもっとも単純な⽅法は，図 3.2中段左のように単⼀のガウス分布で近似 することである．カルマンフィルタのように分布に線形性・ガウス性を仮定する場合にはこの⼿法が取 られるが，より⼀般的な問題を考えたとき確率分布はあらゆる形を取りうるため，単⼀ガウス分布によ る近似では理論値から⼤きく外れてしまうことになる．
図 3.2中段右は，ヒストグラムによる近似を表している．これは，xの取り得る値の範囲を有限個に 分割し，それぞれの区間での⾼さを記録する⽅法である．縦棒の⾯積の総和が 1 となるようにすれば，
このヒストグラムは確率分布の近似となる．この⼿法は⼀⾒合理的であるが，状態ベクトルの次元が
⾼々4 次元程度までしか有効ではない．原因は，状態ベクトルの変数それぞれを分割して格⼦状にヒス トグラムの⾼さを記録しなければならず必要な計算機上のメモリが次元数のべき乗で増⼤することで ある．例えば各変数を10分割するならば，n次元の問題に対しては10n個のヒストグラムの⾼さを記録 しなくてはならない．さらに，ヒストグラムを多重積分する際のプログラミング上の困難も問題となる．
次に，図 3.2下段左に⽰すように，複雑な形状の分布を複数のガウス分布の和で近似する⼿法が考え られる．しかし，この⼿法では逐次計算を繰り返すうちにガウス分布の個数が指数的に爆発してしまう という短所があり，増殖したガウス分布のうち不要なものを取り除く作業も煩雑となる．このガウス和 フィルタは理論的には優れているが，そのような実⽤上の不便さによって普及していない．
最後に図 3.2下段右に⽰すのが，SMCフィルタで⽤いられるモンテカルロ近似である．モンテカルロ 近似では，真の分布に従うように配置された多数の実現値（サンプル）によって分布を表現する．この 実現値の１つ１つを，粒⼦，またはアンサンブルメンバと呼ぶ．実際にN個の実現値を⽤いて予測分布 およびフィルタ分布を表現することを考えたとき，１つ１つの実現値は次のように表される．
𝒙𝑡|𝑡−1(𝑖) [予測分布]
𝒙𝑡|𝑡(𝑖) [フィルタ分布]
ここで，(i)はN個ある実現値のうち何番⽬のものであるかを⽰す．また，下付き添字は，縦棒の左が状 態ベクトルの時刻，右側がどの時刻まで観測データを得られているかを表す．このN個の集団はアンサ ンブルと呼ばれ，式（3.10），（3.11）のように表記される．
𝑋
…|…ol≡ š
𝒙𝑡|𝑡−1(𝑖)›
Bœ𝟏
‰ [予測分布] (3.10)
𝑋
…|…≡ š
𝒙𝑡|𝑡(𝑖)›
Bœ𝟏
‰ [フィルタ分布] (3.11)
また，各粒⼦の重みは常に等しいとする．すなわち，
𝑃 }𝒙… = 𝒙…|…ol(B) ~ = 𝑃 }𝒙… = 𝒙…|…(B)~ = 1
𝑁 (3.12)