知財求人 - 知財ポータルサイト「IP Force」
特開2025-21308人工天体の軌道推定装置、軌道推定方法、及び軌道推定プログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2025021308
(43)【公開日】2025-02-13
(54)【発明の名称】人工天体の軌道推定装置、軌道推定方法、及び軌道推定プログラム
(51)【国際特許分類】
B64G 3/00 20060101AFI20250205BHJP
【ＦＩ】
B64G3/00
【審査請求】未請求
【請求項の数】13
【出願形態】ＯＬ
(21)【出願番号】P 2023125125
(22)【出願日】2023-07-31
(71)【出願人】
【識別番号】523155489
【氏名又は名称】デロイトトーマツリスクアドバイザリー合同会社
(74)【代理人】
【識別番号】110002516
【氏名又は名称】弁理士法人白坂
(72)【発明者】
【氏名】満田 和真
(57)【要約】
【課題】掃天観測による短い円弧の観測からでも、精度よく軌道要素を推定可能な人工天体の軌道推定装置、軌道推定方法、及び軌道推定プログラムを提供する。
【解決手段】掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定方法であって、掃天観測によって得られた、複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出ステップと、掃天観測による観測のうち、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定ステップと、同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出ステップとを備える。
【選択図】図５
【解決手段】掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定方法であって、掃天観測によって得られた、複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出ステップと、掃天観測による観測のうち、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定ステップと、同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出ステップとを備える。
【選択図】図５
【特許請求の範囲】
【請求項1】
掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定装置であって、
前記掃天観測によって得られた、前記複数の人工天体のそれぞれに関するデータを記憶する記憶部と、
演算部であって、
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出部と、
前記掃天観測のうち、異なる観測について前記事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定部と、
同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出部と、
を含む演算部と
を備えることを特徴とする軌道推定装置。
前記掃天観測によって得られた、前記複数の人工天体のそれぞれに関するデータを記憶する記憶部と、
演算部であって、
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出部と、
前記掃天観測のうち、異なる観測について前記事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定部と、
同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出部と、
を含む演算部と
を備えることを特徴とする軌道推定装置。
【請求項2】
前記低精度事後分布算出部は、
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの時刻と位置とから、前記軌道要素を算出し、
算出された前記軌道要素をサンプリングする
ことを特徴とする請求項１に記載の軌道推定装置。
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの時刻と位置とから、前記軌道要素を算出し、
算出された前記軌道要素をサンプリングする
ことを特徴とする請求項１に記載の軌道推定装置。
【請求項3】
前記低精度事後分布算出部は、前記軌道要素を算出する際、
前記複数の人工天体のそれぞれの時刻と位置とから、接線要素として平均時刻における平均位置と平均速度とを算出し、
前記算出された平均位置と平均速度とから前記軌道要素を算出する
ことを特徴とする請求項２に記載の軌道推定装置。
前記複数の人工天体のそれぞれの時刻と位置とから、接線要素として平均時刻における平均位置と平均速度とを算出し、
前記算出された平均位置と平均速度とから前記軌道要素を算出する
ことを特徴とする請求項２に記載の軌道推定装置。
【請求項4】
前記低精度事後分布算出部は、
前記複数の人工天体のそれぞれの明るさから前記天体サイズの低精度事後分布を得る
ことを特徴とする請求項１に記載の軌道推定装置。
前記複数の人工天体のそれぞれの明るさから前記天体サイズの低精度事後分布を得る
ことを特徴とする請求項１に記載の軌道推定装置。
【請求項5】
前記人工天体特定部は、
任意の３視野の観測の組み合わせを選択し、
選択された前記組み合わせについて、前記低精度事後分布を比較し、前記３視野の観測の組み合わせの一致度を算出し、
前記算出された一致度が所定の閾値を超えた場合に、選択された前記組み合わせを同一の人工天体のものであるとみなす
ことを特徴とする請求項１に記載の軌道推定装置。
任意の３視野の観測の組み合わせを選択し、
選択された前記組み合わせについて、前記低精度事後分布を比較し、前記３視野の観測の組み合わせの一致度を算出し、
前記算出された一致度が所定の閾値を超えた場合に、選択された前記組み合わせを同一の人工天体のものであるとみなす
ことを特徴とする請求項１に記載の軌道推定装置。
【請求項6】
前記高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングすることを特徴とする請求項５に記載に軌道推定装置。
【請求項7】
前記高精度事後分布算出部は、前記ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする際、ＭｕｌｔｉＮｅｓｔを用いて実行することを特徴とする請求項６に記載の軌道推定装置。
【請求項8】
前記高精度事後分布算出部は、前記ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのサンプリングを、ＭｕｌｔｉＮｅｓｔを用いて実行する際、Ｐｒｉｏｒを、前記３視野の観測の組み合わせのそれぞれの観測点を得ることによって得ることを特徴とする請求項６に記載の軌道推定装置。
【請求項9】
前記高精度事後分布算出部は、Ｐｒｉｏｒを、前記３視野の観測の組み合わせのそれぞれの観測点を得たのちに、
得られたそれぞれの前記観測点の接線要素をサンプリングし、
前記ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布を足し合わせる
ことを特徴とする請求項８に記載の軌道推定装置。
得られたそれぞれの前記観測点の接線要素をサンプリングし、
前記ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布を足し合わせる
ことを特徴とする請求項８に記載の軌道推定装置。
【請求項10】
前記高精度事後分布算出部は、前記ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのサンプリングを、ＭｕｌｔｉＮｅｓｔを用いて実行する際、Ｐｒｉｏｒを荒い事後分布を用いた複数視野のマッチングにおいて得られた、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの共通部分の面積、又は体積から得られる各パラメータの分布から得ることを特徴とする請求項６に記載の軌道推定装置。
【請求項11】
前記高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする際、前記３視野の観測の組み合わせのすべてについて前記サンプリングを繰り返すことを特徴とする請求項５に記載の軌道推定装置。
【請求項12】
掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定方法であって、
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出ステップと、
前記掃天観測のうち、異なる観測について前記事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定ステップと、
同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出ステップと
を備えることを特徴とする方法。
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出ステップと、
前記掃天観測のうち、異なる観測について前記事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定ステップと、
同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出ステップと
を備えることを特徴とする方法。
【請求項13】
掃天観測によって得られた、複数の人工天体のそれぞれに関するデータを記憶する記憶部と、演算部とを備え、前記掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定装置において、前記軌道推定装置に、
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出機能と、
前記掃天観測のうち、異なる観測について前記事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定機能と、
同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出機能と
を実現させるための軌道推定プログラム。
前記掃天観測によって得られた、前記複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出機能と、
前記掃天観測のうち、異なる観測について前記事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定機能と、
同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出機能と
を実現させるための軌道推定プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、人工天体の軌道推定装置、軌道推定方法、及び軌道推定プログラムに関する。
【背景技術】
【0002】
現在、地球の周りには数多くの宇宙ゴミが分布している。図１に示すように、人工衛星１を運用するにあたり、これらの宇宙ゴミ２を回避する際、宇宙ゴミ２の軌道を推定し、予報円６を避けるように、人工衛星１を本来の軌道３から、一時的に、例えば軌道４、５へ変更する必要がある。一時的な軌道変更により、余分な燃料消費による寿命短縮、軌道変更期間中のビジネスの停止等による経済的な損失が生じる。図１（ａ）、および図１（ｂ）に示すように、予報円６が小さい程一時的な軌道変更の距離は小さく、変更期間が短くなり、したがって経済的損失は小さくなる。精度良く人工天体の軌道を推定することにより、一時的な軌道変更の距離を小さくすることができる。
【0003】
従来、人工衛星や宇宙ゴミといった人工天体の軌道を推定する際、光学観測や電波レーダーによる観測からＫｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓ（ケプラーの法則に基づく軌道要素）を推定していた。光学観測では、電波レーダーと比べて位置決定精度が高い、静止軌道以遠などの高い軌道の人工天体でも観測しやすい等のメリットがある。一方で、距離の情報は得られず、赤経（Ｒｉｇｈｔ Ａｓｃｅｎｓｉｏｎ；ＲＡ）、赤緯（Ｄｅｃｌｉｎａｔｉｏｎ；Ｄｅｃ）の２つの角度のみが得られる。人工天体の運動の角度のみが得られる観測に基づく軌道決定においては、一般にＧａｕｓｓの方式やＧｏｄｄｉｎｇｓの方法といった天体力学に基づく手法が用いられる。近年では、これらの方法では軌道決定が難しい場合でも適用できる方法も開発されている（例えば、特許文献１参照）。特許文献１は、角度のみの観測による初期起動決定の方法を開示している。
【0004】
一方、人工天体を発見、観測する手段として、広視野動画サーベイ観測といった天文学の掃天観測データの利用が注目されている。広視野動画サーベイ観測は、従来、超新星爆発等の突発現象、重力波の光学対応天体の検出、惑星、彗星、地球近接天体等の移動天体が従来の主な観測対象である。広視野動画サーベイ観測は、地上の所定の観測位置から、夜間、上空の広い範囲を短時間の動画撮影により観測し、多数の人工天体を検出する。この観測では、毎晩、広い範囲を短時間の動画で観測することで、一晩に数千の移動天体を検出し、これらの検出された移動天体の大部分は人工天体と考えられる。このような観測データを利用して人工天体を発見し、軌道決定できれば人工天体の安全な運用などに役立てることができる。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特許第６１０６３４１号明細書
【発明の概要】
【発明が解決しようとする課題】
【0006】
しかしながら、従来の光学観測や電波レーダーによる観測からＫｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓを推定する方法では、短い円弧からの軌道決定や観測点同士が離れている場合にロバストな軌道決定を行うことが難しく、また摂動が反映された軌道要素を求めることが出来ない。広視野動画サーベイといった掃天観測では、それぞれの動画観測の時間が短く、したがって人工天体が描く円弧も１度よりも遥かに小さいため、従来の手法を適用して軌道決定を行うことが困難である。
【0007】
上記問題点を鑑み、本発明は、掃天観測による短い円弧の観測からでも、精度よく軌道要素を推定可能な人工天体の軌道推定装置、軌道推定方法、及び軌道推定プログラムを提供することを目的とする。
【課題を解決するための手段】
【0008】
本発明の第１の態様は、掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定装置であって、掃天観測によって得られた、複数の人工天体のそれぞれに関するデータを記憶する記憶部と、演算部であって、掃天観測によって得られた、複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出部と、掃天観測のうち、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定部と、同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出部と、を含む演算部とを備ええることを要旨とする。
【0009】
本発明の第１の態様において、低精度事後分布算出部は、掃天観測によって得られた、複数の人工天体のそれぞれの時刻と位置とから、軌道要素を算出し、算出された軌道要素をサンプリングしてもよい。
【0010】
本発明の第１の態様において、低精度事後分布算出部は、軌道要素を算出する際、複数の人工天体のそれぞれの時刻と位置とから、接線要素として平均時刻における平均位置と平均速度とを算出し、算出された平均位置と平均速度とから軌道要素を算出してもよい。
【0011】
本発明の第１の態様において、低精度事後分布算出部は、複数の人工天体のそれぞれの明るさから天体サイズの低精度事後分布を得てもよい。
【0012】
本発明の第１の態様において、人工天体特定部は、任意の３視野の観測の組み合わせを選択し、選択された組み合わせについて、低精度事後分布を比較し、３視野の観測の組み合わせの共通部分の面積を算出し、算出された面積が所定の閾値を超えた場合に、選択された前記組み合わせを同一の人工天体のものであるとみなしてもよい。
【0013】
本発明の第１の態様において、高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングしてもよい。
【0014】
本発明の第１の態様において、高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする際、ＭｕｌｔｉＮｅｓｔを用いて実行してもよい。
【0015】
本発明の第１の態様において、高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのサンプリングを、ＭｕｌｔｉＮｅｓｔを用いて実行する際、Ｐｒｉｏｒを、３視野の観測の組み合わせのそれぞれの観測点を得ることによって得てもよい。
【0016】
本発明の第１の態様において、高精度事後分布算出部は、Ｐｒｉｏｒを、３視野の観測の組み合わせのそれぞれの観測点を得たのちに、得られたそれぞれの観測点の接線要素をサンプリングし、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布を足し合わせてもよい。
【0017】
本発明の第１の態様において、高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのサンプリングを、ＭｕｌｔｉＮｅｓｔを用いて実行する際、Ｐｒｉｏｒを荒い事後分布を用いた複数視野のマッチングにおいて得られた、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの共通部分の面積、又は体積から得られる各パラメータの分布から得てもよい。
【0018】
本発明の第１の態様において、高精度事後分布算出部は、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする際、３視野の観測の組み合わせのすべてについてサンプリングを繰り返してもよい。
【0019】
本発明の第２の態様は、掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定方法であって、掃天観測によって得られた、複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出ステップと、掃天観測のうち、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定ステップと、同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出ステップとを備えることを要旨とする。
【0020】
本発明の第３の態様は、軌道推定プログラムであって、掃天観測によって得られた、複数の人工天体のそれぞれに関するデータを記憶する記憶部と、演算部とを備え、前記掃天観測に基づき、MCMCを用いたBayes推定によって軌道推定を行う複数の人工天体の軌道推定装置において、軌道推定装置に、掃天観測によって得られた、複数の人工天体のそれぞれの軌道の短い円弧から軌道要素および天体サイズの低精度事後分布を得る低精度事後分布算出機能と、掃天観測による観測のうち、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する人工天体特定機能と、同一の人工天体についての３つ以上の短い円弧の観測を用いて軌道要素の高精度事後分布を得る高精度事後分布算出機能とを実現させることを要旨とする。
【発明の効果】
【0021】
本発明によれば、掃天観測による短い円弧の観測からでも、精度よく軌道要素を推定可能な人工天体の軌道推定装置、軌道推定方法、及び軌道推定プログラムを提供できる。
【図面の簡単な説明】
【0022】
【図1】宇宙ゴミにより人工衛星が軌道変更する様子を説明する模式図である。
【図2】Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓによって記述される人工天体の軌道を説明する模式図であり、図２（ａ）は楕円軌道面での地球と人工天体の位置とａ、およびｅとの関係、図２（ｂ）は３次元における楕円軌道の位置とω、Ｉ、およびΩとの関係図２（ｃ）は楕円軌道面でのｔにおける人工天体の位置とＭ、真近点角ν、および離心近点角Ｅとの関係を示す図である。
【図3】人工天体の観測位置について説明する模式図である。
【図4】本発明の実施形態に係る軌道推定装置の一例を示すブロック図である。
【図5】本実施形態に係る軌道推定方法の一例の全体の流れを説明するフローチャートである。
【図6】本実施形態に係る軌道推定方法における低精度事後分布の算出方法を説明するフローチャートである。
【図7】本発明の実施形態に係る軌道推定装置において計算される、αｉ、δｉと平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇとの関係の一例を示すグラフである。
【図8】本発明の実施形態に係る軌道推定装置におけるＬｏｇ Ｌｉｋｅｌｉｈｏｏｄの計算結果の一例を示す。
【図10】サンプリングされた接線要素の事後分布の一例を示すグラフである。
【図11】Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓの事後分布のグラフである。
【図12】ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布のグラフである。
【図14】サンプリングされた接線要素の事後分布の一例を示すグラフである。
【図15】Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓの事後分布のグラフ
【図16】ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布のグラフ
【図18】サンプリングされた接線要素の事後分布の一例を示すグラフである。
【図19】Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓの事後分布のグラフ
【図20】ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布のグラフ
【図21】人工天体、観測装置、および太陽光の角度の関係を説明する模式図である。
【図22】アルベドの事前分布を示すグラフである。
【図23】計算された事後分布を示すグラフである。
【図24】３視野の観測が全て同じ天体である場合の事後分布を示すグラフである。
【図25】１視野が異なる天体である場合の事後分布を示すグラフである。
【図26】３つの観測点を用いてｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのｐｒｉｏｒを作成した一例を示すグラフである。
【図27】ＭｕｌｔｉＮｅｓｔでサンプリングされたｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布の一例を示すグラフである。
【図28】ＭｕｌｔｉＮｅｓｔでサンプリングされたＫｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓの事後分布の一例を示すグラフである。
【発明を実施するための形態】
【0023】
次に、図面を参照して、本発明の実施形態を説明する。実施形態に係る図面の記載において、同一又は類似の部分には同一又は類似の符号を付している。但し、図面は模式的なものである。
【0024】
又、実施形態は、本発明の技術的思想を具体化するための装置や方法を例示するものであって、本発明の技術的思想は、各構成要素の構成や配置、レイアウト等を下記のものに特定するものでない。本発明の技術的思想は、特許請求の範囲に記載された請求項が規定する技術的範囲内において、種々の変更を加えることができる。
【0025】
（実施形態）
本実施形態に係る軌道推定装置は、掃天観測によって得られた、複数の人工天体の軌道の短い円弧のデータから、精度よく軌道要素を推定することが出来る。
本実施形態に係る軌道推定装置は、掃天観測によって得られた、複数の人工天体の軌道の短い円弧のデータから、精度よく軌道要素を推定することが出来る。
【0026】
本実施形態に係る軌道推定方法は、軌道推定装置において、本実施形態に係る軌道推定方法を実現する処理内容が記載された軌道推定プログラムを実行させることにより実現可能である。
【0027】
まず、人工天体の軌道を記述するＫｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓ、およびｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓ（天の赤道要素）について以下に説明する。Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓは元期ｔにおけるパラメータａ：半長軸半径、ｅ：離心率、ω：近点引数、Ｉ：軌道傾斜角、Ω：昇交点赤経、Ｍ：平均近点角、である。図２に、これらのパラメータによって記述される人工天体２０の軌道を示す。図２（ａ）は楕円軌道面での地球２１と人工天体２０の位置とａ、およびｅとの関係を示す図である。図２（ｂ）は３次元における楕円軌道の位置とω、Ｉ、およびΩとの関係を、近地点２２、遠地点２３、軌道面２４、赤道面２５、人工天体２０の周回方向２６とともに示す図である。図２（ｃ）は楕円軌道面でのｔにおける人工天体の位置とＭ、真近点角ν、および離心近点角Ｅとの関係を示す図である。
【0028】
Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓにより、地球と人工天体の２体問題を前提とした楕円軌道が定まる。一方で、Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓでは真円に近い軌道や赤道面と平行に近い軌道を表現することが難しい。そこで、これらの様な軌道を含め、地球を周回するあらゆる人工天体の軌道を表現することができるｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓも用いられる。ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓは、パラメータａ：半長軸半径、ｅｘ＝ｅｃｏｓ（ω＋Ω）、ｅｙ＝ｅｓｉｎ（ω＋Ω）、ｈｘ＝ｔａｎ（ｉ／２）ｃｏｓ（Ω）、ｈｙ＝ｔａｎ（ｉ／２）ｓｉｎ（Ω）、Ｉν＝Ｍ＋ω＋Ωである。
【0029】
実際には人工天体の軌道は地球重力の不均一性、月や太陽の重力、太陽の輻射圧、上層大気による空気抵抗などによる摂動を受けて変化する。光学や電波などの観測から軌道要素を求める際には、これらの摂動を扱うＳＤＰ４、ＳＧＰ４といった軌道伝播アルゴリズムを用いることで、摂動を考慮したＫｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓを求めることができる。
【0030】
本実施形態では、掃天観測の複数夜のデータから、軌道要素を推定する。掃天観測について、以下に説明する。
【0031】
掃天観測としては、例えば、広視野動画サーベイ観測が挙げられる。広視野動画サーベイ観測は、所定の範囲の視野を望遠鏡に設置された複数のイメージセンサを用いて可視光単色の動画撮影を行う観測方法である。具体的な例としては、約２０平方度の視野を８４枚のＣＭＯＳセンサによってカバーする。１視野について毎秒２フレーム、０．５秒の露光が１２フレーム実施される。即ち、１視野について６秒の動画が取得される。１視野の周りに望遠鏡を微動して撮像を繰り返すことにより、イメージセンサ間の間隙が埋められる。さらに望遠鏡の角度を変えて動画撮影を繰り返すことにより、一晩に一回以上全天を掃天する。
【0032】
広視野動画サーベイ観測によって取得された動画データに対して、１視野の各イメージセンサ内において移動天体の検出が行われる。この検出の際、動画データにはノイズや「ｓｅｅｉｎｇダンス（大気揺らぎによる星像のふらつき）」等が多数含まれるため、ランダムフォレスト分類器によって真に移動天体である確率が所定値以上のものが抽出される。１晩に数千個の移動天体が検出され、その大部分は既知の人工天体、または未知の人工天体であると考えられる。
【0033】
地球近接天体（Ｎｅａｒ Ｅａｒｔｈ Ｏｂｊｅｃｔ；ＮＥＯ）の典型的な見かけの速さは～３ａｒｃｓｅｃ［ｓ－１］である。見かけの速さ～５０ａｒｃｓｅｃ［ｓ－１］は、高度約１万ｋｍに相当する。広視野動画サーベイの動画データには、この速度よりも移動速度が速い移動天体が映り込んでいるが、この速さを下回る人工天体が、宇宙ゴミや人工衛星といった、ＮＥＯのコンタミネーションとして検出される。
【0034】
検出される人工天体の円弧は、１視野６秒間の観測においては、１度よりもはるかに小さく、例えば、高度１万ｋｍでは、６秒間の円弧は１／１２度である。
【0035】
広視野動画サーベイで得られる人工天体のデータは、１視野につき、所定のフレームの観測位置、即ち角度と、明るさである。人工天体の観測位置は、各フレームの撮像時刻ｔｉにおけるＲＡ（αｉ）、およびＤｅｃ（δｉ）、明るさは可視光等級ｍｉとして得られる。ここで、ｉは１からフレーム数の間の整数であり、例えば１２フレーム観測の場合は１から最大１２の間の整数である。
【0036】
図３に、望遠鏡３３と、撮像時刻ｔ０、ｔ１、およびｔ２における人工天体３０のそれぞれの観測位置を示す。人工天体３０は軌道３４に沿って運動している。撮像時刻ｔ０における人工天体３０の観測位置はＲＡ（α０）、およびＤｅｃ（δ０）で表される。撮像時刻ｔ１、およびｔ２についても同様である。
【0037】
本実施形態に係る軌道推定装置を図４を参照しながら説明する。図４に示すように、本実施形態に係る軌道推定装置４０は、各種の演算実行のためのＣＰＵ４１、処理用のプログラム、人工天体のデータを含むデータ等の記憶のための記憶部４２、Ｉ／Ｏ（インプット・アウトプットインターフェース）４３、表示部４４、入力部４５から構成される。
【0038】
軌道推定装置４０は、パーソナルコンピュータ（ＰＣ）、メインフレーム、ワークステーション、クラウドコンピューティングシステム等、種々の電子計算機（計算リソース）であってよい。Ｉ／Ｏ４３は通信（送受信）用のインターフェース、バッファ等である。
【0039】
表示部４４は、ディスプレイ等の表示装置であり、ＣＰＵ４１による演算結果等を表示する。入力部４５は、キーボード、マウス等の入力装置であり、ユーザからの入力を受け付け、ＣＰＵ４１に送信する装置によって実現される。
【0040】
図４のブロック図に、ＣＰＵ４１内の機能部を示す。ＣＰＵ４１の各機能部をソフトウェアにより実現する場合、ＣＰＵ４１は各機能を実現するソフトウェアであるプログラムの命令を実行することで実現する。詳細には、低精度事後分布算出部４１１、人工天体特定部４１２、高精度事後分布算出部４１３等を備える。
【0041】
低精度事後分布算出部４１１は、Ｉ／Ｏ４３において取得した掃天観測による観測データに基づき、掃天観測の１視野の観測から軌道要素と天体サイズについての低精度事後分布を算出する。
【0042】
人工天体特定部４１２は、１視野の観測からＭＣＭＣサンプリングによって得られた荒い事後分布を用いて、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する、すなわちマッチングを行う。
【0043】
高精度事後分布算出部４１３は、掃天観測の複数視野の観測から軌道要素についての高精度の事後分布を算出する。
【0044】
本発明の実施形態に係る軌道推定方法を以下に説明する。本実施形態に係る軌道推定方法は、広視野動画サーベイを含む掃天観測において得られる短い円弧についての光学観測が複数夜に渡って得られる場合に、人工天体を同定し軌道決定する方法である。図５に、本実施形態に係る軌道推定方法の一例の全体の流れを説明するフローチャートを示す。本実施形態に係る軌道推定方法は、図５に示すように、以下の４ステップで実施される。
【0045】
ステップＳ５０１において、掃天観測による観測データを取得する。
【0046】
ステップＳ５０２において、掃天観測の１視野の観測から軌道要素と天体サイズについての精度の低い事後分布を算出する。本実施形態においては、以下に定義する接線軌道要素をＭＣＭＣサンプリングすることによって軌道要素についての事後分布を得る。さらに、軌道要素の事後分布と観測された明るさから天体サイズについての事後分布を算出する。
【0047】
ステップＳ５０３において、１視野の観測からＭＣＭＣサンプリングによって得られた荒い事後分布を用いて、異なる観測について事後分布同士を比較することによって同一の人工天体についての観測を特定する、すなわちマッチングを行う。
【0048】
ステップＳ５０４において、掃天観測の複数視野の観測から軌道要素についての高精度の事後分布を算出する。
【0049】
【0050】
【0051】
（接線要素を算出）
ステップＳ６０１において、接線要素（ｏｓｃｕｌａｔｉｎｇ ｐａｒａｍｅｔｅｒｓ）を算出する。まず、１視野の観測から得られる時刻ｔｉ、及び位置αｉ、δｉから、以下の関係式を基に平均時刻ｔａｖｇにおける平均位置αａｖｇ、δａｖｇと平均速度ｖα，ａｖｇ、ｖδ，ａｖｇとを算出する。平均時刻ｔａｖｇ、及び平均位置αａｖｇ、δａｖｇは以下の式で与えられる。
ステップＳ６０１において、接線要素（ｏｓｃｕｌａｔｉｎｇ ｐａｒａｍｅｔｅｒｓ）を算出する。まず、１視野の観測から得られる時刻ｔｉ、及び位置αｉ、δｉから、以下の関係式を基に平均時刻ｔａｖｇにおける平均位置αａｖｇ、δａｖｇと平均速度ｖα，ａｖｇ、ｖδ，ａｖｇとを算出する。平均時刻ｔａｖｇ、及び平均位置αａｖｇ、δａｖｇは以下の式で与えられる。
【0052】
【数1】
【0053】
【数2】
【0054】
【数3】
【0055】
また、αｉ、δｉより、各時刻における速度ｖα，ｉ、ｖδ，iが
【0056】
【数4】
【0057】
【数5】
【0058】
で与えられ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇは、
【0059】
【数6】
【0060】
【数7】
【0061】
で与えられる。
【0062】
接線要素は、平均時刻ｔａｖｇにおける地心距離ｒ、地心距離ｒの時間微分ｄｒ／ｄｔ、平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの６つのパラメータとして定義される。ここで、ｒ及びｄｒ／ｄｔは地球中心座標（Ｇｅｏｃｅｎｔｒｉｃ Ｆｒａｍｅ）、他の４パラメータは観測者中心座標系（Ｔｏｐｏｃｅｎｔｒｉｃ Ｆｒａｍｅ）で定義されることに注意されたい。
【0063】
図７にαｉ、δｉと平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇとの関係の一例を示す。図７に示すグラフにおいて、グラフ中のグレースケールで示されるデータ７１はαｉ、δｉであり、グレースケールの濃淡は時刻ｔｉに対応する。横軸と縦軸はそれぞれ、αｉ、δｉの平均値からの差分、グレースケールの濃淡は平均時刻からの差分である。グラフの中央のデータ７２は平均位置αａｖｇ、δａｖｇであり、データ７２からデータ７３までのベクトルは、（ｖα，ａｖｇ，ｖδ，ａｖｇ）に０．２５秒を乗じたベクトルである。
【0064】
（接線要素から軌道要素を得る手順）
ステップＳ６０２において、軌道要素（Ｋｅｐｌｅｒｉａｎ Ｅｌｅｍｅｎｔｓ）が接線要素と同様に６つのパラメータから構成され、各パラメータが互いに変数変換の関係であることから、得られた接線要素を基に、軌道要素を算出する。以下にその手順を説明する。
ステップＳ６０２において、軌道要素（Ｋｅｐｌｅｒｉａｎ Ｅｌｅｍｅｎｔｓ）が接線要素と同様に６つのパラメータから構成され、各パラメータが互いに変数変換の関係であることから、得られた接線要素を基に、軌道要素を算出する。以下にその手順を説明する。
【0065】
地球中心での天体の位置を計算する。地心距離ｒ、観測者中心でのＲＡαＧｅｏ、ＤｅｃδＧｅｏとすると、地球中心から観測者までのベクトルは
【0066】
【数8】
【0067】
で与えられ、観測者から天体までのベクトルは
【0068】
【数9】
【0069】
で与えられ、地球中心から天体までのベクトルは
【0070】
【数10】
【0071】
で与えられる。ここで、ベクトルＳをdについて解くことが出来、ベクトルｓを決定できる。これにより、Ｓすなわち地球中心座標系における天体の位置を求めることができる。
【0072】
接線要素ｒ、ｄｒ／ｄｔ、αａｖｇ、δａｖｇ、ｖα，ａｖｇ、ｖδ，ａｖｇが与えられると、地球中心座標系での天体の位置を求めることができ、さらにこれらの微小変位から天体の速度を求めることができる。すなわち、地球中心座標での位置ベクトル
【0073】
【数11】
【0074】
および速度ベクトル
【0075】
【数12】
【0076】
が得られる。ひとたび位置ベクトルｐＧｅｏ、および速度ベクトルｖＧｅｏが得られると、角運動量ベクトルを用いて地球―天体の２体問題における軌道要素を求めることができる。
【0077】
（ＭＣＭＣにおけるＬｏｇ Ｌｉｋｅｌｉｈｏｏｄの定義）
ステップＳ６０３において、接線軌道要素をＭＣＭＣ（マルコフ連鎖モンテカルロ法、Ｍａｒｃｏｖ Ｃｈａｉｎ Ｍｏｎｔｅ Ｃａｌｒｏ）サンプリングする際のＬｏｇ Ｌｉｋｅｌｉｈｏｏｄを定義する。なお、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄとは、最尤推定であり、統計学において、与えられたデータからそれが従う確率分布の母数を点推定する方法である。また、ＭＣＭＣとは、求める確率分布を均衡分布として持つマルコフ連鎖を作成することによって確率分布のサンプリングを行う種々のアルゴリズムの総称である。
ステップＳ６０３において、接線軌道要素をＭＣＭＣ（マルコフ連鎖モンテカルロ法、Ｍａｒｃｏｖ Ｃｈａｉｎ Ｍｏｎｔｅ Ｃａｌｒｏ）サンプリングする際のＬｏｇ Ｌｉｋｅｌｉｈｏｏｄを定義する。なお、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄとは、最尤推定であり、統計学において、与えられたデータからそれが従う確率分布の母数を点推定する方法である。また、ＭＣＭＣとは、求める確率分布を均衡分布として持つマルコフ連鎖を作成することによって確率分布のサンプリングを行う種々のアルゴリズムの総称である。
【0078】
まず、各サンプルについて、接線軌道要素から得られる軌道要素を用いて、観測時刻ｔｉにおけるＲＡ（αｍ，ｉ）、Ｄｅｃδｍ，ｉを計算する。観測のαｉ、δｉは、それぞれ、モデルのαｍ，ｉ、δｍ，ｉの周りに観測の位置決定精度の範囲で正規分布すると仮定する。これらの差分は、
【0079】
【数13】
【0080】
【数14】
【0081】
で表される。ここで、ＲＡ方向の差分についてはＣＯＳ（δ）を乗じる必要があることに注意されたい。
【0082】
さらに、観測誤差をσとすると、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄは
【0083】
【数15】
【0084】
で表される。ただし、ｆ（ｘ，σ）は、平均０、分散σ２の正規分布のｘにおける確率密度関数の自然対数を与える。
【0085】
接線軌道要素から得られる軌道要素は摂動が考慮されていない一方で、これらの軌道要素からαｍ，ｉ、δｍ，ｉを計算する際には摂動を考慮し、ＳＧＰ４等のアルゴリズムを用いることもできる。摂動考慮の有無による差異は、ここで考慮されるタイムスケールの範囲においては、ｇｌｏｂａｌな位置のｏｆｆｓｅｔとして現れる。したがって、上式をそのまま適用すると、ｇｌｏｂａｌ ｏｆｆｓｅｔによってＬＬが非常に小さくなってしまう。そこで、常識にｇｌｏｂａｌ ｏｆｆｓｅｔを吸収する項を追加する。
【0086】
まず、Δαｉ、Δδｉをｇｌｏｂａｌ成分とｌｏｃａｌ成分とに分ける。ｇｌｏｂａｌ成分Δαｇｌｏｂ、Δδｇｌｏｂはｉについて平均を取り、以下の式
【0087】
【数16】
【0088】
【数17】
【0089】
に基づき計算する。ｌｏｃａｌ成分Δαｌｏｃ，ｉ、Δδｌｏｃ，ｉはｇｌｏｂａｌ成分からの変位であって、以下の式
【0090】
【数18】
【0091】
【数19】
【0092】
で与えられるとする。ｇｌｏｂａｌ成分は、摂動考慮の有無による典型的な差異σｇｌｏｂを標準偏差とする正規分布、ｌｏｃａｌ成分は観測の位置決定精度を標準偏差とする正規分布に従うと仮定すると、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄは
【0093】
【数20】
【0094】
で表される。右辺第３項と４項が摂動考慮の有無に起因するｇｌｏｂａｌ ｏｆｆｓｅｔを吸収する項と解釈できる。
【0095】
実際の観測の際の一例として、観測の位置決定精度は、～０．３ａｒｃｓｅｃである。また、σｇｌｏｂの値は典型値として～２０ａｒｃｍｉｎ程度と推定される。以下では、軌道要素から人工天体の位置の計算を行う際、ＳＧＰ４によって摂動を考慮するものとする。
【0096】
（ＭＣＭＣの初期値の設定）
ステップＳ６０４において、ＭＣＭＣの初期値を設定する。接線要素のサンプリングにおけるＭＣＭＣの初期値は、平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇについては、観測から得られる値をそのまま用いることができる。一方、接線要素ｒ、ｄｒ／ｄｔについては、初期値を外部的に与える必要がある。
ステップＳ６０４において、ＭＣＭＣの初期値を設定する。接線要素のサンプリングにおけるＭＣＭＣの初期値は、平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇについては、観測から得られる値をそのまま用いることができる。一方、接線要素ｒ、ｄｒ／ｄｔについては、初期値を外部的に与える必要がある。
【0097】
接線要素ｒ、ｄｒ／ｄｔの初期値は、これらのｇｒｉｄ上のいくつかの点において軌道要素を計算し、閉じた軌道を与えるものから採用する。人工天体については、これらの取りうる範囲が限られている。例えば、接線要素ｒは高度１００ｋｍ程度から月の軌道高度３８万ｋｍ程度、ｄｒ／ｄｔは±５ｋｍｓ－１程度が現実的に取り得る範囲である。この範囲で適当な数（例えば２０程度）のｇｒｉｄを切った上、他の接線要素については観測から得られる値を使用し、軌道要素の計算を行う。閉じた軌道を与える点、またはＬＬが大きな点をＭＣＭＣの初期値とすることができる。
【0098】
軌道要素の事後分布を接線要素のサンプリングによって実現することの利点は、観測において、ＬＬが接線要素に対して凸関数となることである。ＬＬが凸関数となるため、ＭＣＭＣによってロバストなサンプリングが可能である。図８に、ｒ、ｄｒ／ｄｔのグリッド上の点においてＬＬを計算した例を示す。図８のグラフの、横軸がｒ、縦軸がｄｒ／ｄｔ、ドットは計算を行った点、グレースケールの濃度はＬＬの値を示し、グラフ中の、濃いグレーによって示される領域８０が、ＬＬの大きな領域である。
【0099】
（事前分布の設定）
ステップＳ６０５において、後述するＭＣＭＣサンプリングを行うにあたり、事前分布を設定することができる。ベイズ推定においては、事前分布を決定することから開始する。事前分布は、過去の文献、予備データ、経験等によりおよその分布が得られている場合には、そのような分布を今回の事前分布として用いることができる。事前分布とは、前から持っている情報のことを指し、尤度は取得した情報である。
ステップＳ６０５において、後述するＭＣＭＣサンプリングを行うにあたり、事前分布を設定することができる。ベイズ推定においては、事前分布を決定することから開始する。事前分布は、過去の文献、予備データ、経験等によりおよその分布が得られている場合には、そのような分布を今回の事前分布として用いることができる。事前分布とは、前から持っている情報のことを指し、尤度は取得した情報である。
【0100】
まず、接線要素のサンプリングにおいては、接線要素ｒ、ｄｒ／ｄｔの２次元分布、および平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇそれぞれの１次元分布を用いることができる。ただし、事後分布への事前分布の影響は小さい。接線要素ｒ、ｄｒ／ｄｔの事前分布については、既知の人工天体等の軌道要素の分布から計算されたもの等を用いることができる。平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇについては、例えば１２ｆｒａｍｅの観測から得られる標準偏差を用いて正規分布を用いることができる。
【0101】
（接線要素のサンプリング）
ステップＳ６０６において、ＭＣＭＣで接線要素のサンプリングを行う際には、以下の（１）～（３）の３パターンでそれぞれサンプリングを行う。
（１）接線要素ｒの１パラメータでのサンプリング
（２）接線要素ｒ、ｄｒ／ｄｔの２パラメータでのサンプリング
（３）接線要素ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリング
（１）では、ｄｒ／ｄｔ＝０を仮定し、平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇを観測値で固定する。（２）では、サンプリングされない平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇを観測値で固定する。（３）では、サンプリングされない平均位置αａｖｇ、δａｖｇを観測値で固定する。（２）および（３）は離心率の大きな軌道を含む様々な場合のサンプリングに適している一方、対象天体の軌道が円軌道に近い場合に、離心率を確率的に過大評価しやすい。（１）も、対象天体の軌道が円軌道に近い場合に離心率を確率的に過大評価しやすい傾向を完全には排除できないが、ｄｒ／ｄｔ＝０であるため、接線要素ｒが真の値に近い場合、真の軌道要素に近い値を与え得る。
ステップＳ６０６において、ＭＣＭＣで接線要素のサンプリングを行う際には、以下の（１）～（３）の３パターンでそれぞれサンプリングを行う。
（１）接線要素ｒの１パラメータでのサンプリング
（２）接線要素ｒ、ｄｒ／ｄｔの２パラメータでのサンプリング
（３）接線要素ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリング
（１）では、ｄｒ／ｄｔ＝０を仮定し、平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇを観測値で固定する。（２）では、サンプリングされない平均位置αａｖｇ、δａｖｇ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇを観測値で固定する。（３）では、サンプリングされない平均位置αａｖｇ、δａｖｇを観測値で固定する。（２）および（３）は離心率の大きな軌道を含む様々な場合のサンプリングに適している一方、対象天体の軌道が円軌道に近い場合に、離心率を確率的に過大評価しやすい。（１）も、対象天体の軌道が円軌道に近い場合に離心率を確率的に過大評価しやすい傾向を完全には排除できないが、ｄｒ／ｄｔ＝０であるため、接線要素ｒが真の値に近い場合、真の軌道要素に近い値を与え得る。
【0102】
（１）接線要素ｒの１パラメータでのサンプリング
接線要素のうち、ｒの１パラメータに対してＭＣＭＣサンプリングを行う。上述の通り、得られた初期値から、メトロポリス・ヘイスティングス法（Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ，ＭＨ法）サンプリング、またはギブス（Ｇｉｂｂｓ）サンプリングを用い、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄの式（数１９）からサンプリングを行う。事後分布に収束するまでのステップは１００に満たない程度であり、ｂｕｒｎ－ｉｎにおいてはステップは２００あれば十分である。また、初期値を複数設定し、初期値のそれぞれを起点に並列にサンプリングすることで計算時間を短縮できる。
接線要素のうち、ｒの１パラメータに対してＭＣＭＣサンプリングを行う。上述の通り、得られた初期値から、メトロポリス・ヘイスティングス法（Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ，ＭＨ法）サンプリング、またはギブス（Ｇｉｂｂｓ）サンプリングを用い、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄの式（数１９）からサンプリングを行う。事後分布に収束するまでのステップは１００に満たない程度であり、ｂｕｒｎ－ｉｎにおいてはステップは２００あれば十分である。また、初期値を複数設定し、初期値のそれぞれを起点に並列にサンプリングすることで計算時間を短縮できる。
【0103】
提案分布の幅は、ｂｕｒｎ－ｉｎにおける採択確率が適当な幅に収まるように設定することができる。実際の観測の一例として、６秒１２フレームの観測データの場合、採択確率を２０から４０％程度となるように提案分布の幅を設定することで、安定して収束することが確認された。図９（ａ）、図９（ｂ）、図９（ｃ）に、ｂｕｒｎ－ｉｎにおけるｔｒａｃｅ、Ｌｉｋｅｌｉｈｏｏｄ、採択確率の一例を示す。図９（ｂ）に示すＬｉｋｅｌｉｈｏｏｄは最大値が１となるように規格化されている。図９（ａ）、図９（ｂ）、図９（ｃ）のそれぞれのデータのグレースケールは、色の濃さが、並列サンプリングにおける、すなわち異なる初期値から開始した各ｃｈａｉｎを示す。３０ステップ目頃には初期値の影響がなくなっており、図９（ｃ）に示すように、採択確率の減少も止まっている。
【0104】
図１０に、サンプリングされた接線要素の事後分布の一例を示す。破線Ａに囲まれたグラフは、各パラメータの１次元の事後分布であり、破線Ｂで囲まれたグラフは、それぞれのパラメータペアの２次元の事後分布である。破線Ｃで囲まれたグラフは、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄを２次元上で示したものである。破線Ｂ、および破線Ｃで囲まれたグラフにおける点は、初期値を示す。１パラメータのサンプリングのため、ｒ以外は固定されている。ｒの決定精度は数万ｋｍである。
【0105】
接線要素がサンプリングされると、それに対応する軌道要素を得ることができる。図１１に、Ｋｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓの事後分布のグラフを、図１２にｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布のグラフを示す。図１０と同様に、破線Ａに囲まれたグラフは、各パラメータの１次元の事後分布、破線Ｂで囲まれたグラフは、それぞれのパラメータペアの２次元の事後分布、破線Ｃで囲まれたグラフは、Ｌｏｇ Ｌｉｋｅｌｉｈｏｏｄを２次元上で示したものである。図中の点線は、既知の人工天体で位置と速度によってマッチしたものの値を示す。接線要素空間では、ｒ以外は１値であるが、軌道要素空間ではそれぞれの接線要素セットに対応した軌道要素の値をとるため、全てのパラメータが広がった分布を持つ。
【0106】
図１０～図１２に示す天体は既知の人工天体と位置と速度とがマッチし、その離心率がｅ＝０．００４３３と０に近い。円軌道に近いため、近地点引数ωが不定になり、したがって平均近点角Ｍも不定となる。図１１では、これらのパラメータの分布がマッチした既知天体の値、すなわち真値に近い値をカバーしていない。これは天体の軌道が円軌道に近いことに起因する。実際、図１２に示すように、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布では全てのパラメータの事後分布が既知天体の値をカバーしている。
【0107】
（２）ｒ、ｄｒ／ｄｔの２パラメータでのサンプリング
ｒ、ｄｒ／ｄｔの２パラメータでのサンプリングについても、１パラメータでのサンプリングと同様に行う。ｒ、ｄｒ／ｄｔの２パラメータでのサンプリングの結果の一例を図１３～図１６に示す。図１３～図１６の詳細は、１パラメータでのサンプリングにおける図９～図１２の詳細と同様である。
ｒ、ｄｒ／ｄｔの２パラメータでのサンプリングについても、１パラメータでのサンプリングと同様に行う。ｒ、ｄｒ／ｄｔの２パラメータでのサンプリングの結果の一例を図１３～図１６に示す。図１３～図１６の詳細は、１パラメータでのサンプリングにおける図９～図１２の詳細と同様である。
【0108】
（３）ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリング
ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリングについても、１パラメータでのサンプリングと同様に行う。ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリングの結果の一例を図１７～図２０に示す。図１７～図２０の詳細は、１パラメータでのサンプリングにおける図９～図１２の詳細と同様である。
ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリングについても、１パラメータでのサンプリングと同様に行う。ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリングの結果の一例を図１７～図２０に示す。図１７～図２０の詳細は、１パラメータでのサンプリングにおける図９～図１２の詳細と同様である。
【0109】
（天体サイズの事後分布の計算）
ステップＳ６０７において、接線要素のサンプリングののち、天体の明るさから天体のサイズの事後分布を求める。人工天体の明るさは、人工天体までの距離、人工天体のサイズ、アルベド、形状、回転等の条件に依存する。人工天体のアルベドをρ、断面積をＡ、観測者から人工天体までの距離をｄ、反射の指向性を表す関数をＦ（ψ）とすると、人工天体の観測等級は、
ステップＳ６０７において、接線要素のサンプリングののち、天体の明るさから天体のサイズの事後分布を求める。人工天体の明るさは、人工天体までの距離、人工天体のサイズ、アルベド、形状、回転等の条件に依存する。人工天体のアルベドをρ、断面積をＡ、観測者から人工天体までの距離をｄ、反射の指向性を表す関数をＦ（ψ）とすると、人工天体の観測等級は、
【0110】
【数21】
【0111】
と表される。ただし、ｍ．は太陽の視等級であり、可視光（ＳＤＳＳｇバンド）では―２６．４７である。ここで、人工天体の反射の指向性について、球形モデルで太陽光が当たっている方向にのみ等方散乱すると仮定すると、Ｆ（ψ）＝１／２πである。さらに、人工天体の典型的なサイズをsとし、Ａ＝πｓ２と定義する。図２１に示すように、人工天体２１０が球形のとき、観測装置２１１から見た人工天体への視線２１２と太陽への視線２１３のなす角をβとすると、人工天体２１０が太陽光２１４によって照らされている部分の角度（位相角）γはβと等しい。したがって、人工天体のうち太陽に照らされている面積は、
【0112】
【数22】
【0113】
である。上式のＡをＡｓｈｉｎｅに置き換えると、
【0114】
【数23】
【0115】
【数24】
【0116】
となる。したがって、人工天体のサイズは、
【0117】
【数25】
【0118】
で表される。
【0119】
上式より、距離ｄ、およびアルベドρが決まれば、観測等級ｍから天体サイズを求めることが出来る。距離ｄは接線要素すなわち軌道要素が決まれば、軌道計算を行うことで計算できる。アルベドρについては事前分布であると仮定する。アルベドρの事前分布は実際の観測において観測されたもののうち、既知の天体とマッチするものの明るさを用いて求める。同じ天体の明るさの変動のうち、距離と位相角以外は天体の回転、姿勢の変化などによるものであるが、これらを全てアルベドの項に押し付け、上式を用い、ρの変動幅を取得する。ρの平均値は典型値である０．１とすることで図２２に示すようにρの分布が得られる。このρの分布をそのまま事前分布として用いてもよく、このρの分布に正規分布を当てはめたものを事前分布としてもよい。このρの分布に正規分布を当てはめたものを図２２にグラフ中の点線として示す。以降では、このρの分布に正規分布を当てはめたものを事前分布とする。
【0120】
観測等級、サンプリングされた接線要素、およびアルベドの事前分布から、天体サイズの事後分布を得る。観測等級は１視野の平均ｍａｖｇとして
【0121】
【数26】
【0122】
を用い、式（数２４）より計算する。
【0123】
図２３に得られた事後分布を示す。天体サイズの事後分布は、接線要素と同様、１パラメータ、２パラメータ、および４パラメータのそれぞれの場合について得られる。
【0124】
（荒い事後分布を用いた複数視野のマッチング）
図５に示すステップＳ５０３において、１視野の観測から、ＭＣＭＣによって得られた事後分布を使用して、複数の視野同士のマッチングを行う。まず、任意の３視野の観測の組み合わせを選択する。この組み合わせにおいて、１、２、および４パラメータのそれぞれのＭＣＭＣによって得られた事後分布同士を比較する。ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓにサイズを加えた７つのパラメータの事後分布について、３視野の観測が重なる部分、すなわち事後分布の共通部分の面積、又は体積を計算し、得られた共通部分の面積、又は体積を一致度と呼ぶことにする。事後分布の共通部分の面積とは、７つのパラメータの１次元の事後分布の共通部分の面積である。また、事後分布の共通部分の体積とは、７つのパラメータの事後分布の７次元空間内での共通部分の体積である。一致度として、共通部分の面積、又は体積のいずれを用いても構わない。
図５に示すステップＳ５０３において、１視野の観測から、ＭＣＭＣによって得られた事後分布を使用して、複数の視野同士のマッチングを行う。まず、任意の３視野の観測の組み合わせを選択する。この組み合わせにおいて、１、２、および４パラメータのそれぞれのＭＣＭＣによって得られた事後分布同士を比較する。ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓにサイズを加えた７つのパラメータの事後分布について、３視野の観測が重なる部分、すなわち事後分布の共通部分の面積、又は体積を計算し、得られた共通部分の面積、又は体積を一致度と呼ぶことにする。事後分布の共通部分の面積とは、７つのパラメータの１次元の事後分布の共通部分の面積である。また、事後分布の共通部分の体積とは、７つのパラメータの事後分布の７次元空間内での共通部分の体積である。一致度として、共通部分の面積、又は体積のいずれを用いても構わない。
【0125】
事後分布の一致度は、３視野の観測がすべて同じ天体についてのものであれば大きく、そうでない場合は小さくなる。一致度が適当な閾値を超えた場合、高精度の事後分布の取得に進む。７つのパラメータの中でも軌道面を決定するｈｘ、ｈｙにおいて、３視野の観測が全て同じ天体についてのものである場合とそうでない場合の差が大きい傾向にある。図２４に３視野の観測が全て同じ天体である場合、図２５に１視野が異なる場合の、それぞれの事後分布と共通部分の例を示す。図２４（ａ）～図２４（ｆ）、および図２５（ａ）～図２５（ｆ）はｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓ、図２４（ｇ）、および図２５（ｇ）はサイズを加えた７つのパラメータの事後分布の例である。各グラフのヒストグラムは左側の縦軸で表される確率密度分布であり、グレースケールの色の濃さは異なる視野であることを示す。共通部分はドットで示しており、右側の縦軸で表される。縦の点線はマッチした既知天体の値を示す。
【0126】
（高精度の事後分布の取得）
ステップＳ５０４において、ＲＡ、Ｄｅｃの観測値を用いてｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする。このサンプリングはＮｅｓｔｅｄ Ｓａｍｐｌｉｎｇの一種であるＭｕｌｔｉＮｅｓｔを用いて行う。ｌｏｇ ｌｉｋｅｌｉｈｏｏｄが凸関数とならないため、接線要素をサンプリングするメリットは得られない。そこで、ここではｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓそのものをサンプリングする。
ステップＳ５０４において、ＲＡ、Ｄｅｃの観測値を用いてｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする。このサンプリングはＮｅｓｔｅｄ Ｓａｍｐｌｉｎｇの一種であるＭｕｌｔｉＮｅｓｔを用いて行う。ｌｏｇ ｌｉｋｅｌｉｈｏｏｄが凸関数とならないため、接線要素をサンプリングするメリットは得られない。そこで、ここではｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓそのものをサンプリングする。
【0127】
ＭｕｌｔｉＮｅｓｔはｌｏｇ ｌｉｋｅｌｉｈｏｏｄが非凸の場合やピークが多数ある場合においても、事後分布への収束が比較的ロバストであるメリットがある。この長所を生かし、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓサンプリングを比較的ロバストに実施できる。Ｍｕｌｔｉ ｐｅａｋへの対応はロバストな軌道推定に欠かせない。なぜなら、例えば、全ての視野でそれなりに当てはまっている場合と、ある視野（１２フレーム）のみで当てはまりが良い場合にｌｏｇ ｌｉｋｅｌｉｈｏｏｄが大きくなる、すなわち、ｍｕｌｔｉ ｐｅａｋとなるためである。さらに、視野同士の観測時間が、例えば１日程度、離れている場合に、パラメータを少し変化させたときにｌｏｇ ｌｉｋｅｌｉｈｏｏｄが大きく変化する、すなわちｍｕｌｔｉ ｐｅａｋ同士の間隙が深く、ＭＣＭＣのようなランダムウォークサンプリングでは、ｍｕｌｔｉ ｐｅａｋ間を移動することが難しい。
【0128】
また、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓそのものをサンプリングすることで、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓからモデルのＲＡ，Ｄｅｃ計算にＳＤＰ４等の摂動を扱うアルゴリズムを用いることが出来る。これにより、最終的に得られる軌道要素の事後分布に摂動の効果を反映させることが出来る。
【0129】
（Ｐｒｉｏｒの設定）
ＭｕｌｔｉＮｅｓｔを用いるにあたり、サンプリングを行うパラメータのｐｒｉｏｒ（事前分布）を計算する。Ｐｒｉｏｒは、荒い事後分布を用いた複数視野のマッチングにおいて得られた、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの共通部分の面積、又は体積から得られる各パラメータの分布から得てもよく、又は、３視野の観測から得てもよい。
ＭｕｌｔｉＮｅｓｔを用いるにあたり、サンプリングを行うパラメータのｐｒｉｏｒ（事前分布）を計算する。Ｐｒｉｏｒは、荒い事後分布を用いた複数視野のマッチングにおいて得られた、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの共通部分の面積、又は体積から得られる各パラメータの分布から得てもよく、又は、３視野の観測から得てもよい。
【0130】
３視野の観測から、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのｐｒｉｏｒを得る場合、視野ｊの観測αｊ，ｉ、δｊ，ｉから代表点αｊ，ｏ、δｊ，ｏを抽出する。これは各視野の最初や最後の観測、平均などでかまわない。この段階で３つの観測点（ｊ＝１，２，３）が得られる。
【0131】
さらに、これら３つの観測点を用いて、以下の（１）～（３）の３パターンでそれぞれＭＣＭＣで接線要素のサンプリングを行ったのと同様の手順で接線要素のサンプリングを行ってもよい。
（１）接線要素ｒの１パラメータでのサンプリング
（２）接線要素ｒ、ｄｒ／ｄｔの２パラメータでのサンプリング
（３）接線要素ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリング
ここで、接線要素を用いることができるのは、各視野の観測点を一つに絞っているため、ｍｕｌｔｉ ｐｅａｋ性が軽減されるためである。
（１）接線要素ｒの１パラメータでのサンプリング
（２）接線要素ｒ、ｄｒ／ｄｔの２パラメータでのサンプリング
（３）接線要素ｒ、ｄｒ／ｄｔ、平均速度ｖα，ａｖｇ、ｖδ，ａｖｇの４パラメータでのサンプリング
ここで、接線要素を用いることができるのは、各視野の観測点を一つに絞っているため、ｍｕｌｔｉ ｐｅａｋ性が軽減されるためである。
【0132】
接線要素のサンプリングののち、さらに、１、２、および４パラメータのサンプリングから得られたｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布を足し合わせてもよい。この事後分布（確率密度関数）をＭｕｌｔｉＮｅｓｔに与えるｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのｐｒｉｏｒとする。実際には、累積確率密度関数の逆関数（ｐｅｒｃｅｎｔｉｌｅ ｐｏｉｎｔ ｆｕｎｃｔｉｏｎ，ｐｐｆ）を与える。
【0133】
図２６に、３つの観測点を用いてｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓのｐｒｉｏｒを作成した一例を示す。図２６に示すグラフのうち、左側のグラフはいずれも、確率密度関数であり、グラフ中のヒストグラム２６１は接線要素のサンプリングから得られたもので、ステップ２６２は数値的な補正を施したものである。図２６に示すグラフのうち、右側のグラフはいずれも、確率密度関数から累積確率密度関数（ｃｕｍｕｌａｔｉｖｅ ｄｅｎｓｉｔｙ ｆｕｎｃｔｉｏｎ）を計算したものである。
【0134】
３視野の観測αｊ，ｉ、δｊ，ｉすべてを用いて、ｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓをサンプリングする。ここでは、軌道伝搬にＳＤＰ４も用いるため、ｌｏｇ ｌｉｋｅｌｉｈｏｏｄとしては、摂動によるｏｆｆｓｅｔを考慮しない式を採用する。ＭｕｌｔｉＮｅｓｔでサンプリングすることでｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布を得る。図２７に、ＭｕｌｔｉＮｅｓｔでサンプリングされたｅｑｕｉｎｏｃｔｉａｌ ｅｌｅｍｅｎｔｓの事後分布の一例を示す。また、図２８に、ＭｕｌｔｉＮｅｓｔでサンプリングされたＫｅｐｌｅｒｉａｎ ｅｌｅｍｅｎｔｓの事後分布の一例を示す。例えばａについては、分布の幅は１０ｍオーダーとなっており、高精度で事後分布が求められていることがわかる。
【0135】
以上、本発明はここでは記載していない様々な実施形態等を含むことは勿論である。したがって、本発明の技術的範囲は上記の説明から妥当な特許請求の範囲に係る発明特定事項によってのみ定められるものである。
【符号の説明】
【0136】
１ 人工衛星
２ 宇宙ゴミ
３、４、５ 軌道
６ 予報円
２０、３０ 人工天体
２１ 地球
２２ 近地点
２３ 遠地点
２４ 軌道面
２５ 赤道面
２６ 周回方向
３３ 望遠鏡
３４ 軌道
４０ 軌道推定装置
４１ ＣＰＵ
４２ 記憶部
４３ Ｉ／Ｏ（インプット・アウトプットインターフェース）
４４ 表示部
４５ 入力部
４１１ 低精度事後分布算出部
４１２ マッチング部
４１３ 高精度事後分布算出部
７１、７２、７３ データ
２１０ 人工天体
２１１ 観測装置
２１２、２１３ 視線
２１４ 太陽光
２ 宇宙ゴミ
３、４、５ 軌道
６ 予報円
２０、３０ 人工天体
２１ 地球
２２ 近地点
２３ 遠地点
２４ 軌道面
２５ 赤道面
２６ 周回方向
３３ 望遠鏡
３４ 軌道
４０ 軌道推定装置
４１ ＣＰＵ
４２ 記憶部
４３ Ｉ／Ｏ（インプット・アウトプットインターフェース）
４４ 表示部
４５ 入力部
４１１ 低精度事後分布算出部
４１２ マッチング部
４１３ 高精度事後分布算出部
７１、７２、７３ データ
２１０ 人工天体
２１１ 観測装置
２１２、２１３ 視線
２１４ 太陽光
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【図26】
【図27】
【図28】
