JP6401976B2 - 宇宙物体同定方法、宇宙物体同定装置及び宇宙物体同定プログラム - Google Patents

Description

本発明は、宇宙物体を同定する技術に関する。

人類が宇宙活動を開始した１９５０年代以降、宇宙空間に打上げられるロケット及び人工衛星は年々増加している。使用済みロケット及び人工衛星は、スペースデブリ（宇宙ごみ）となって宇宙空間に留まる。近年これら宇宙物体の急激な増加が大きな問題となっており、宇宙物体をレーダまたは光学望遠鏡で観測してその軌道を把握し、運用中の人工衛星との接近回避予測に用いるということが行われるようになってきた。

宇宙物体の軌道把握については米国が最大の観測システムを有しており、現在約17000個の宇宙物体情報が公開され、その軌道情報が日々更新されている。宇宙物体情報はカタログデータと言われ、以下の項目が包含される。
・識別番号、識別名称
・発生日（打上日など）
・レーダ反射断面積 RCS（Rader Cross Section）
・軌道情報
軌道情報はケプラー要素に近い形式で表現され、ＴＬＥ（Two Line Element）と呼ばれる。

日本国内ではＴＬＥに基づく観測が実施されている。光学観測ではサーベイ観測が行われている。サーベイ観測では、得られた宇宙物体画像を目視でＴＬＥ情報と比較し、どのカタログ番号の物体であるか判別する。この処理を“同定”という。光学観測の場合、観測できるのは主に静止軌道付近、すなわち地上から見るとあまり動かない宇宙物体であるので、目視による同定もそれほど困難性はない。

一方、レーダは低高度（数100km〜1000km程度）の宇宙物体観測に使用される。この場合、地上から宇宙物体が見える時間は軌道１周回約100分中10分内外である。レーダ観測方式には追尾観測とサーベイ観測（フェンスモード）があり、現在、国内では主に追尾観測を実施している。

特開２０１１−１５７０３０号公報 特開２０１３−１２９３０７号公報

しかしながら、より多くの宇宙物体、特にカタログ化されていない宇宙物体まで観測するにはサーベイ観測が必要であり、その場合は得られた観測データをＴＬＥと照合して同定する必要がある。ＴＬＥ情報は日々更新されるといっても17000個全てが毎日更新されているわけではなく、放置されているＴＬＥ情報は時間経過と共に精度が劣化する。特に低高度物体については大気抵抗力の不確定さによる誤差が無視できない。したがって、目視で観測データとＴＬＥ情報の照合を行うのは極めて困難である。

本発明は、一側面として、識別された宇宙物体の宇宙物体の軌道情報の誤差と観測誤差とを考慮して、観測された宇宙物体の同定を行う技術を提供することを目的とする。

１つの態様では、宇宙物体同定方法は、コンピュータが次の処理を行う。すなわち、コンピュータが、所定地点で観測された、宇宙空間に存在する対象物体の位置に関する時系列の観測情報を取得する。コンピュータが、格納部から取得した軌道情報と誤差情報とに基づいて、識別された物体の軌道予測における進行方向の位置のマイナス側の位置誤差が最大となる軌道伝播による、観測の開始時刻における識別された物体の第１位置ベクトルと、識別された物体の軌道予測における進行方向の位置のプラス側の位置誤差が最大となる軌道伝播による、観測の終了時刻における識別された物体の第２位置ベクトルとを算出する。コンピュータが、第１位置ベクトルと、第２位置ベクトルとに基づいて、観測情報が観測された観測時間内で、識別された物体が軌道伝播する軌道伝播範囲を算出する。軌道情報は、宇宙空間に存在する物体であってそれぞれが識別された物体の軌道に関する情報である。誤差情報は、軌道の誤差に関する情報である。コンピュータは、観測時間内における、識別された物体毎の軌道伝播範囲に含まれる対象物体の数に基づいて、対象物体を同定する。

本発明の一側面によれば、識別された宇宙物体の宇宙物体の軌道情報の誤差と観測誤差とを考慮して、観測された宇宙物体の同定を行うことができる。

本実施形態における宇宙物体同定装置のブロック図を示す。 本実施形態における宇宙物体同定装置の一例を示す。 本実施形態におけるＴＬＥ軌道誤差テーブルの一例を示す。 本実施形態における全体の処理フローを示す。 本実施形態における観測設備性能によるフィルタリング処理（Ｓ１１）フローを示す。 本実施形態における観測位置ベクトル計算処理（Ｓ１２）フローを示す。 本実施形態におけるレンジデータの補間を説明するための図である。 本実施形態における赤道面座標系と地表面座標系を説明するための図である。 本実施形態における幾何学的可視判定によるフィルタリング処理（Ｓ１３）フローを示す。 本実施形態における幾何学的可視判定によるフィルタリング処理（Ｓ１３）フローを示す。 本実施形態におけるエレベーション角による可視判定を説明するための図である。 本実施形態における位置ベクトルの予測存在領域を説明するための図である。 本実施形態における観測値と予測値のマッチング処理（Ｓ１４）フローを示す。 本実施形態における観測値と予測値のマッチング処理（Ｓ１４）フローを示す。 本実施形態における予測位置ベクトル「方向」の予測存在領域を説明するための図である。 本実施形態における、２つの基準方向単位ベクトルで定義する座標系を説明するための図である。 本実施形態における、４つの単位ベクトルによって形成される領域の包含条件を説明するための図である。 本実施形態における同定結果ファイルの一例を示す。 本実施形態における同定結果を表示する画面の一例を示す。 本実施形態におけるプログラムを実行するコンピュータのハードウェア環境の構成ブロック図の一例である。

図１は、本実施形態における宇宙物体同定装置のブロック図を示す。宇宙物体同定装置１は、取得部２、格納部６、算出部３、同定部４を含む。

格納部６には、軌道情報７と、誤差情報８とが格納されている。軌道情報７は、宇宙空間に存在する物体であってそれぞれが識別された物体の軌道に関する情報である。誤差情報８は、該軌道の誤差に関する情報である。格納部６の一例として、記憶装置２０が挙げられる。軌道情報７の一例として、カタログＤＢ２１が挙げられる。誤差情報８の一例として、ＴＬＥ軌道誤差テーブル２４が挙げられる。

取得部２は、所定地点で観測された、宇宙空間に存在する対象物体の位置に関する時系列の観測情報を取得する。取得部２の一例として、制御部１２が挙げられる。観測情報の一例として、観測データ２２が挙げられる。

算出部３は、格納部６から取得した軌道情報７と誤差情報８とに基づいて、観測情報が観測された観測時間内で、識別された物体が軌道伝播する軌道伝播範囲を算出する。算出部３の一例として、Ｓ２およびＳ３（Ｓ１３，Ｓ１４−１〜Ｓ１４−８）の処理を行う制御部１２が挙げられる。

同定部４は、観測時間内における、識別された物体毎の軌道伝播範囲に含まれる対象物体の数に基づいて、対象物体を同定する。同定部４の一例として、Ｓ３，Ｓ１４の処理を行う制御部１２が挙げられる。

このように構成することにより、識別された宇宙物体の宇宙物体の軌道情報の誤差と観測誤差とを考慮して、観測された宇宙物体の同定を行うことができる。

算出部３は、軌道情報７と誤差情報８とに基づいて、少なくとも、第１位置ベクトルと第２位置ベクトルとを算出する。第１位置ベクトルは、識別された物体の軌道予測における進行方向の位置のマイナス側の位置誤差が最大となる軌道伝播による観測開始時刻での前記識別された物体の位置ベクトルである。第２位置ベクトルは、識別された物体の軌道予測における進行方向の位置のプラス側の位置誤差が最大となる軌道伝播による観測終了時刻での識別された物体の位置ベクトルである。算出部３は、誤差情報に基づいて、第１位置ベクトルと第２位置ベクトルが成す面の傾斜誤差を算出する。算出部３は、第１位置ベクトルと、第２位置ベクトルと、傾斜誤差とに基づいて、識別された物体の軌道伝播範囲を算出する。

このように構成することにより、観測時間における識別された物体の軌道伝播範囲を算出することができる。

同定部５は、対象物体は、複数の識別された物体のそれぞれの軌道伝播範囲において、計測された対象物体の数が最も多い軌道伝播範囲に対応する識別された物体であると同定する。

このように構成することにより、計測された対象物体の数に応じて、宇宙物体の同定を行うことができる。

以下では、本実施形態について、さらに詳述する。本実施形態は、以下の事項（Ｉ）（ＩＩ）を前提にする。なお、以下では、宇宙物体を「ＳＯ」と称する場合がある。
（Ｉ）ＳＯの同定では、その時保持されているカタログデータベース（米国の公開データを日々反映）に基づき、観測値とＴＬＥによる予測値が最も近いＳＯを特定する。
（ＩＩ）観測値は、レーダ等により所定時間追跡して観測された所定の宇宙物体についての時系列の観測データより得られるものであり、具体的には、以下に示すように、レンジと角度データとが得られているものとする。ただし、レンジと角度データの観測時刻が一致していることは前提としない。
レンジ： 観測局とＳＯと間の距離データ
角度データ：観測局からＳＯへの方向の方位角AZ、仰角EL

図２は、本実施形態における宇宙物体同定装置の一例を示す。宇宙物体同定装置１１は、制御部１２、記憶装置２０を含む。制御部１２により実行される同定処理は、観測設備性能によるフィルタリング処理（Ｓ１）、幾何学的可視判定によるフィルタリング処理（Ｓ２）、観測値と予測値のマッチング処理（Ｓ３）を含む。

観測設備性能によるフィルタリング処理（Ｓ１）では、観測局の対象の観測設備が観測可能な高度の上限、軌道傾斜角の下限、RCSの下限が予め指定されている。制御部１２は、これらの指定された情報に基づいて、カタログデータベース（ＤＢ）２１のカタログデータをサーチして、観測設備の性能上、観測できないＳＯを除外する。

幾何学的可視判定によるフィルタリング処理（Ｓ２）では、制御部１２は、ＴＬＥに基づく軌道伝播を行う。軌道伝播とは、ある時刻での宇宙物体の位置、速度等に基づいて、異なる時刻での宇宙物体の位置、速度等を予測することである。それから、制御部１２は、その軌道伝播に基づいて、観測時間帯において対象観測局の位置から幾何学的に可視の可能性があるか判定し、対象観測局の位置から幾何学的に可視できないＳＯを除外する。Ｓ２の処理では、軌道伝播にて予測される誤差が考慮される。

観測値と予測値のマッチング処理（Ｓ３）では、制御部１２は、Ｓ１及びＳ２のフィルタ処理で残ったＳＯに対して以下の処理を行う。制御部１２は、ＴＬＥによる予測誤差と観測誤差を考慮し、観測開始〜観測終了までの期間においてＳＯの位置ベクトルが存在し得る範囲を求める。制御部１２は、不明の宇宙物体についての時系列の観測データから計算される位置ベクトルがその存在し得る範囲に包含されるかを判定し、包含される観測データ数が最も多いＳＯを抽出する。制御部１２は、不明の宇宙物体は、その抽出されたSOであると同定し、同定結果を出力する。

ここで、観測値と予測値のマッチング処理（Ｓ３）では、以下の２つの理由により、赤道面座標における位置ベクトルに対して、処理が行われる。なお、赤道面座標とは、地心を原点とし、地球自転軸の北極方向を＋Ｚ軸（すなわちＸＹ面が赤道面）、春分点方向を＋Ｘ軸とした座標系であり、ＳＯの軌道を表現するのに使用される。
・ＴＬＥに基づく軌道予測の誤差において、優勢（dominant）な誤差は、進行方向位置誤差である。赤道面座標においては、この進行方向位置誤差は単純に軌道上位相の進み／遅れで表現することができる。観測データ（レンジ、角度）では、この進行方向位置誤差を単純に表現することが困難であり、場合によっては観測されていたＳＯが不可視になる場合もある。
・マッチングの対象は複数のＳＯのＴＬＥ（赤道面座標）と観測データであるが、地表面座標の観測データを赤道面座標に変換する方が座標変換の回数は少なくて済む。

記憶装置２０には、本実施形態で用いる計算に使用するデータが格納される。記憶装置２０には、カタログＤＢ２１、観測データ２２、観測設備条件情報２３、ＴＬＥ軌道誤差テーブル２４、プログラム実行条件情報２５、出力データ３０が格納される。出力データ３０は、フィルタ１後カタログデータ（ワークファイル）２６、観測位置ベクトルファイル（ワークファイル）２７、フィルタ２後予測位置データファイル（ワークファイル）２８、同定結果ファイル２９を含む。

カタログＤＢ２１は、日々、外部機関から源泉取得されて、更新されている全カタログデータを格納する。カタログデータは、以下のデータを含む。
・カタログ識別番号／名称
・ＴＬＥ軌道パラメータ
・RCS（Rader Cross Section，レーダ反射断面積）

観測データ２２は、観測局で取得した以下を含む時系列のデータであって、ファイルに格納されている。観測データ２２は、レーダ等により追跡して観測された不明の宇宙物体（対象ＳＯ）についての時系列の観測データである。
・観測時刻、レンジデータ （ｔi、ρi ） ｉ＝１，，ｎ
・観測時刻、角度データ （ｔj、AZj、ELj ） ｊ＝１，，ｍ
観測設備条件情報２３として、以下の情報が予め設定されている。
・局位置（測地緯度、経度、高度） φ、λ、ｈ
・観測可能範囲 高度上限 Ｈmax
軌道傾斜角下限 ｉmin
ＲＣＳ下限 RCSmin

ＴＬＥ軌道誤差テーブル２４は、図３に示すように、所定の軌道傾斜角(deg)の範囲毎の、位置誤差、速度誤差を格納する。位置誤差は、高度方向ΔH、進行方向ΔL、クロストラック方向ΔCで示される。速度誤差は、高度方向ΔVH、進行方向ΔVL、クロストラック方向ΔVCで示される。ＴＬＥ軌道誤差テーブル２４に格納される位置誤差および速度誤差の設定値は、例えば欧州宇宙機関（ＥＳＡ）の調査結果を参考にして、予め設定されている。

プログラム実行条件情報２５は、入力により指定される以下のプログラム実行条件である。
・幾何学的可視判定のＥＬ下限値 ELCUT
・大気抵抗による高度低下率の誤差率 ε（％）
・観測値と予測値マッチング時に考慮する誤差マージン
レンジ誤差 Δρobs
方向角度誤差 Δθobs
（誤差マージンには、必ずしも観測設備性能に基づく標準偏差やバイアスではなく、観測値と予測値マッチングの際にマージンとして考慮する値を設定することもできる。）

フィルタ１後カタログデータファイル２６は、観測設備性能によるフィルタリング処理（Ｓ１）が実行された結果、残ったＳＯのカタログ識別名、ＴＬＥ軌道パラメータを格納するワークファイルである。

観測位置ベクトルファイル２７は、観測データ（レンジ、角度）から計算した赤道面座標における位置ベクトルの時系列データを格納するワークファイルである。

フィルタ２後予測位置データファイルは、幾何学的可視判定によるフィルタリング処理（Ｓ２）を実行した結果、残ったＳＯについて、TLEに基づく以下の軌道予測値を格納するワークファイルである。
（i） 誤差考慮なし： 観測開始時刻の位置ベクトル、ケプラー要素
（ii） 誤差考慮なし： 観測終了時刻の位置ベクトル、ケプラー要素
（iii） ＋誤差考慮： 観測開始時刻の位置ベクトル、ケプラー要素
（iv） ＋誤差考慮： 観測終了時刻の位置ベクトル、ケプラー要素
（y） −誤差考慮： 観測開始時刻の位置ベクトル、ケプラー要素
（vi） −誤差考慮： 観測終了時刻の位置ベクトル、ケプラー要素
ここで＋／−誤差考慮とは、進行方向位置の予測誤差について＋／−最大の状態を意味する。また、ケプラー要素とは、ケプラーの法則に基づく軌道要素であり、ａ：軌道長半径、ｅ：離心率、ｉ：軌道傾斜角、Ω：昇交点赤経、ω：近地点引数、ｆ：真近点離角（代わりに、Ｍ：平均近点離角を使う場合が多い）の要素を含む。

同定結果ファイル２９には、カタログデータから得られたＳＯの存在し得る範囲に包含される観測データのＳＯの計測数Ｎinc＞０である対象ＳＯについて、カタログ識別名とＮincが出力される（Ｎincについて降順）。また、同定結果ファイル２９は、同定結果として表示画面への表示に使用される。Ｎincは、具体的には、ＴＬＥ予測値と誤差マージンから導出される「観測開始時刻〜観測終了時刻においてＳＯが存在しうる領域」に包含される、対象ＳＯについての観測データのＳＯの計測数である。同定結果ファイル２９には、「ＳＯが存在し得る領域」の境界値情報が併せて出力される。

図４は、本実施形態における全体の処理フローを示す。観測設備性能によるフィルタリング処理（Ｓ１１）は、不明ＳＯの観測データに依存しないので、カタログＤＢ２１のカタログデータの更新直後に実行されるものとする。不明ＳＯ観測直後に「観測位置ベクトル計算」処理（Ｓ１２）以降が実行される。以下、Ｓ１１〜Ｓ１４の処理について説明する。まず、制御部１２は、観測設備性能によるフィルタリング処理（Ｓ１１）を実行する。Ｓ１１の処理について図５を用いて説明する。

図５は、本実施形態における観測設備性能によるフィルタリング処理（Ｓ１１）フローを示す。ここでは、制御部１２は、カタログＤＢ２１からＴＬＥ軌道パラメータとＲＣＳを１つずつ読み込む。制御部１２は、ＴＬＥ軌道パラメータから、ＳＯの軌道で、地球に最も近づく点、すなわち近地点（ペリジ）の高度を計算する（Ｓ１１−１）。制御部１２は、式(1-1)に示すように算出したペリジ高度およびＴＬＥ軌道パラメータと、観測設備条件情報２３の観測可能範囲パラメータとを比較する（Ｓ１１−２〜Ｓ１１−４）。
式(1-1)の３つの条件を全て満たす場合（Ｓ１１−２で「Ｙｅｓ」、Ｓ１１−３で「Ｙｅｓ」、Ｓ１１−４で「Ｙｅｓ」）、制御部１２は、フィルタ１後カタログデータファイル２６へカタログ識別名とＴＬＥ軌道パラメータとを出力する（Ｓ１１−５）。制御部１２は、Ｓ１１−１〜Ｓ１１−５の処理を、カタログＤＢ２１に格納された全カタログ数分繰り返す。

次に、図４に示すように、制御部１２は、観測位置ベクトル計算処理（Ｓ１２）を行う。ここでは、制御部１２は、観測データ（レンジ、AZ、EL）２２から赤道面座標におけるＳＯ位置ベクトルを計算する。このとき、レンジとＡＺ／ＥＬデータの観測時刻が一致しているとは限らないことから、以下に説明するように、最初に同一時刻の（レンジ、AZ、EL）セットを揃える処理が行われる。Ｓ１２の処理について図６を用いて説明する。

図６は、本実施形態における観測位置ベクトル計算処理（Ｓ１２）フローを示す。まず、制御部１２は、観測データ２２に含まれるレンジデータが時刻昇順に並んでいるか、およびAZ/ELデータが時刻昇順に並んでいるかをチェックし、それぞれレンジ用時刻昇順ファイル、AZ/EL用時刻昇順ファイルに格納する（Ｓ１２−１）。

次に、制御部１２は、レンジデータの補間処理を行う（Ｓ１２−２〜Ｓ１２−９）。観測時刻はAZ/ELデータを基準にし、レンジデータはその基準にしたAZ/ELデータの時刻前後のデータで補間される。レンジデータの方を補間する理由は、球面座標のAZ/ELが、直接、線形補間するパラメータとして適さないからである。図７に示すように、AZ/ELデータの観測時刻とレンジデータの観測時刻から、レンジデータの線形補間が行われ、AZ/EL観測時刻におけるレンジデータが算出される。算出されたAZ/EL観測時刻におけるレンジデータを用いて、レンジベクトル作成時刻が算出される。ここで、Ｓ１２−２〜Ｓ１２−９で用いるパラメータの表記は、以下の通りである。
・TAE： AZ/EL観測時刻
・TR： レンジ観測時刻
・ρ： 時刻TRのレンジ
・T0： 観測時刻TRの１つ前のレンジ観測時刻
・ρ0： 時刻TR0のレンジ

制御部１２は、AZ/EL用時刻昇順ファイルからAZ/EL観測データを１つずつ読み込む。制御部１２は、読み込んだAZ/EL観測データを用いて、AZ/EL観測時刻、観測局からSOへの方向の方位角および仰角をそれぞれ、TAE、AZ、ELに設定する（Ｓ１２−２）。さらに、制御部１２は、レンジ用時刻昇順ファイルからレンジ観測データを１つ読み込む。制御部１２は、読み込んだレンジ観測データを用いて、観測時刻TR、時刻TRのレンジρを設定する（Ｓ１２−３）。

観測時刻TRが微小判定値δの範囲内でTAEと一致している場合は（Ｓ１２−４で「Ｙｅｓ」）、制御部１２は、ρm＝ρと設定し（Ｓ１２−９）、Ｓ１２−１０の処理へ進む。

観測時刻TRが微小判定値δの範囲内でTAEと一致していない場合（Ｓ１２−４で「Ｎｏ」）、制御部１２は、TR＞TAEとなる最初の時刻を順番に探す（Ｓ１２−５で「Ｙｅｓ」）。TR＞TAEとなる最初の時刻が見つからなかった場合（Ｓ１２−５で「Ｎｏ」）、制御部１２は、レンジ用時刻昇順ファイルから、次のレンジ観測データを１つ読込み（Ｓ１２−３）、Ｓ１２−４以降の処理を行う。

TR＞TAEとなる最初の時刻が見つかった場合（Ｓ１２−５で「Ｙｅｓ」）、制御部１２は、１つ前のレンジデータが存在するかを判定する（Ｓ１２−６）。１つ前のレンジデータがない場合は（Ｓ１２−６で「Ｎｏ」）、制御部１２は、AZ/EL用時刻昇順ファイルから次のAZ/EL観測データを読み込み（Ｓ１２−２）、Ｓ１２−３以降の処理を行う。

１つ前のレンジデータが存在する場合（Ｓ１２−６で「Ｙｅｓ」）、制御部１２は、式（2-1）により、Ｓ１２−３で読み込んだレンジデータと１つ前のレンジデータとの線形補間でTAEにおけるレンジを計算する（Ｓ１２−６で「Ｙｅｓ」、Ｓ１２−７）。
ρm＝ρ0＋（ρ−ρ0）（TAE−T0）／（TR−T0） (2-1)

Ｓ１２−８またはＳ１２−９の処理後、制御部１２は、（TAE、ρm、AZ、EL）を１つのセットとし、位置ベクトル計算処理に進む。

位置ベクトル計算処理では、まず、赤道面座標系XYZ-Iに対して、地上局位置を基準にして、地表面座標系XYZ-Sが図８のように定義される。XYZ-S座標系において、レンジと方向角度（ρm、AZ、EL）が与えられた場合、制御部１２は、式（2-2）により、レンジベクトルρを計算する（Ｓ１２−１０）。ここで、AZ角定義は東をＸ、北をＹ方向としている。以下では、ベクトルについては、パラメータに下線をつけて表記する場合がある。

制御部１２は、観測時刻TAEにおいて、局位置ベクトルＳとレンジベクトルρとをそれぞれ、赤道面座標系へ変換する（Ｓ１２−１１）。

制御部１２は、式（2-3）により、観測時刻TAEにおける観測位置ベクトルＲobsを計算し、観測位置ベクトルファイル２７に格納する（Ｓ１２−１２）。ここで、<I>は赤道面座標系における表現を示す。

制御部１２は、レンジ用時刻昇順ファイルの読み込みカウンタを１つ戻す（Ｓ１２−１４）。これは、レンジ用時刻昇順ファイルから次のAZ/EL観測データを読み込んだとき、最新の読み込み済みのレンジ観測時刻から処理を開始するためである。

次に、図４に示すように、制御部１２は、幾何学的可視判定によるフィルタリング処理（Ｓ１３）を実行する。制御部１２は、ＴＬＥ軌道パラメータに基づいて、観測設備性能によるフィルタリング処理（Ｓ１２）を通過した各ＳＯが、観測期間において対象局から可視の可能性があるかを判定する。制御部１２は、最初に、ＴＬＥ軌道伝播による誤差を考慮し、観測開始時刻および観測終了時刻における予測位置ベクトルを計算する。次に、制御部１２は、予測位置ベクトルと局位置ベクトルからレンジベクトルを計算し、観測設備性能によるフィルタリング処理（Ｓ１２）を通過した各ＳＯが、観測基地局から、幾何学的に可視であるか判定する。Ｓ１３の処理について図９を用いて説明する。

図９は、本実施形態における幾何学的可視判定によるフィルタリング処理（Ｓ１３）フローを示す。制御部１２は、観測位置ベクトルファイル２７から、観測開始時刻Ｔｓおよび観測終了時刻Ｔｅを抽出する（Ｓ１３−１）。制御部１２は、抽出した時刻Ｔｓ，Ｔｅそれぞれにおける局位置ベクトルＳｓ，Ｓｅを赤道面座標に変換する（Ｓ１３−２）。

制御部１２は、誤差を考慮した予測位置ベクトル計算を行う（Ｓ１３−３〜Ｓ１３−７）。中高度ＳＯの可視／不可視に最も大きく影響するのは、軌道予測時の進行方向位置誤差である。そこで、以下の（１）〜（６）の軌道予測位置を算出し、いずれか１点でも幾何学的に可視であるならば、そのＳＯは観測されたＳＯである可能性がある、と考える。
（１）誤差を考慮しない軌道伝播による観測開始時刻Ｔsでの予測位置
（２）誤差を考慮しない軌道伝播による観測終了時刻Ｔeでの予測位置
（３）進行方向位置誤差＋側最大となる軌道伝播による観測開始時刻Ｔsでの予測位置
（４）進行方向位置誤差＋側最大となる軌道伝播による観測終了時刻Ｔeでの予測位置
（５）進行方向位置誤差−側最大となる軌道伝播による観測開始時刻Ｔsでの予測位置
（６）進行方向位置誤差−側最大となる軌道伝播による観測終了時刻Ｔeでの予測位置

制御部１２は、フィルタ１後カタログデータファイル２６からカタログデータＴＬＥ０を読み込む（Ｓ１３−３）。ここで、カタログデータのＴＬＥ軌道パラメータが以下であるとする。
・平均運動 ｎ
・軌道長半径 ａ＝（μ／ｎ²）^1/3 μ：地心重力定数
・抗力項 Ｂ*
ＴＬＥ軌道パラメータを入力して軌道伝播による予測位置の算出には、一般に公開されている軌道計算用アルゴリズムのツール（ＳＧＰ４）を使用する。ここで、ＳＧＰ４は、Simplified General Perturbations Satellite Orbit Model 4の略称である。誤差を考慮しない軌道伝播を行う場合（（１）、（２））は、カタログデータの値がそのまま入力される。

進行方向位置±誤差の軌道伝播を行う場合（（３）〜（６））は、ＴＬＥ軌道パラメータのｎ、Ｂ*について誤差が考慮される。

まず、制御部１２は、ＴＬＥ軌道誤差テーブル２４から進行方向の速度誤差を読み出し、式(3-1)により平均運動誤差の大きさδｎを計算する（Ｓ１３−４）。
δｎ＝｜−（３／ａ）ΔVL｜[rad/sec] (3-1)
さらに、Ｂ*についてε％の誤差を考慮するが、制御部１２は、δｎによる進行方向誤差と符号が同じになるように設定する。

まとめると、進行方向位置±誤差の軌道伝播を行う時のＴＬＥ軌道パラメータ初期値設定は以下の表のようになる。

以上の軌道伝播により得られた（１）〜（６）のＳＯの予測位置ベクトルＲとその時のケプラー要素（ａ、ｅ、ｉ、Ω、ω、Ｍ）を以下のsuffixで表す。
Ｒij、 （ａ、ｅ、ｉ、Ω、ω、Ｍ）ij
ｉ＝０：誤差なし、 ｐ：＋誤差、 ｍ：−誤差
ｊ＝ｓ：観測開始時刻、 ｅ：観測終了時刻

これらの処理は、具体的には、制御部１２は、ＴＬＥ０を初期値として、ＳＧＰ４を用いて、観測開始時刻Ｔsから観測終了時刻ＴeまでのＳＯの予測位置ベクトルとケプラー要素（（１）（２））とを計算する（Ｓ１３−５）。

制御部１２は、ＴＬＥ０に＋誤差を考慮した軌道要素を初期値として、ＳＧＰ４と上記表の式(3-2)を用いて、観測開始時刻Ｔsから観測終了時刻ＴeまでのＳＯの予測位置ベクトルとケプラー要素（（３）（４））とを計算する（Ｓ１３−６）。

制御部１２は、ＴＬＥ０に−誤差を考慮した軌道要素を初期値として、ＳＧＰ４と上記表の式(3-3)を用いて、観測開始時刻Ｔsから観測終了時刻ＴeまでのＳＯの予測位置ベクトルとケプラー要素（（５）（６））とを計算する（Ｓ１３−７）。

次に、制御部１２は、エレベーション角による可視判定を行う（Ｓ１３−８〜Ｓ１３−１４）。Ｓ１３−１，Ｓ１３−２で算出したように、観測開始／終了時刻における赤道面座標の局位置ベクトルは各ＳＯに依存しないので、予め算出しておく。

ここでは地球の扁平を無視し、図１０のようにレンジベクトルρと局位置ベクトルＳから算出される仰角（EL）とエレベーションカット角ELCUTとの比較により可視が判定される。
ρ＝Ｒ−Ｓ (3-4)
EL [deg] ＝９０−cos^-1（Ｓ・ρ／｜Ｓ｜｜ρ｜）×180／π (3-5)
制御部１２は、（１）〜（６）のＲijについてELを計算する（Ｓ１３−８，Ｓ１３−１０，Ｓ１３−１２）。計算したＥＬが１つでもEL＞ELCUTとなれば（Ｓ１３−９で「Ｙｅｓ」、Ｓ１３−１１で「Ｙｅｓ」またはＳ１３−１３で「Ｙｅｓ」）、制御部１２は、次の処理を行う。すなわち、制御部１２は、フィルタ２後予測位置データファイル２８へ、カタログ識別名、（１）〜（６）のＳＯの予測位置ベクトルＲij、ケプラー要素（ａ、ｅ、ｉ、Ω、ω、Ｍ）ijを出力する（Ｓ１３−１４）。

制御部１２は、フィルタ１後カタログデータ数分、Ｓ１３−３〜Ｓ１３−１４の処理を繰り返す。

次に、制御部１２は、図４に示すように、観測値と予測値のマッチング処理（Ｓ１４）を行う。幾何学的可視判定によるフィルタリング処理（Ｓ１３）を通過した各ＳＯ（予測ＳＯ候補）について、観測開始時刻Ｔｓから観測終了時刻Ｔｅまでにおける予測値からその予測ＳＯ候補が存在しうる領域（予測存在領域）が設定される。予測存在領域を、図１１を用いて説明する。ここで予測存在領域は「軌道１周回分を超えない」ことが前提となる。

制御部１２は、観測位置ベクトル計算処理（Ｓ１２）で算出された対象ＳＯの位置ベクトルの時系列データで示される対象ＳＯが、予測ＳＯ候補の予測存在領域に包含されているかどうかを判定し、包含されているデータ個数をカウントする。Ｓ１４の処理について図１２を用いて説明する。

図１２は、本実施形態における観測値と予測値のマッチング処理（Ｓ１４）フローを示す。以下では、まず予測ＳＯ候補の予測存在領域の算出方法について記述し、次に位置ベクトル包含判定について記述する。予測ＳＯ候補の存在領域の設定と、対象ＳＯが予測ＳＯ候補の予測存在領域に包含されているかの判定とは、地心距離と位置ベクトル方向に分けて行われる。地心距離とは、地球の中心からＳＯまでの距離をいう。

制御部１２は、予測ＳＯ候補の位置ベクトルの予測存在領域の設定を行う（Ｓ１４−１〜Ｓ１４−８）。ここでは、予測ＳＯ候補の位置ベクトルの予測存在領域の設定（Ｓ１４−２〜Ｓ１４−７）と、予測存在領域の位置ベクトル方向角度の存在範囲の設定（Ｓ１４−８）とが行われる。

制御部１２は、まず、フィルタ２後予測位置データファイル２８から、１つの予測ＳＯ候補のデータを取得し、以下６つの位置ベクトルとケプラー要素を読み込む（Ｓ１４−１）。
Ｒij、（ａ、ｅ、ｉ、Ω、ω、Ｍ）ij
ｉ＝０：誤差なし、 ｐ：＋誤差、 ｍ：−誤差
ｊ＝ｓ：観測開始時刻、 ｅ：観測終了時刻

次に、制御部１２は、予測ＳＯ候補の地心距離の存在範囲を設定する（Ｓ１４−２〜Ｓ１４−７）。ここでは、制御部１２は、６つの位置ベクトルＲijからそれぞれ地心距離Ｒijを算出し、式(4-1)により予測ＳＯ候補の地心距離の最大値と最小値を求める（Ｓ１４−２）。

次に、制御部１２は、観測時間帯において、可視できる観測範囲に、地球から最も遠距離の点である遠地点（アポジ）または地球から最も近距離の点である近地点（ペリジ）が含まれるかを判定するため、次の処理を行う。すなわち、制御部１２は、ＭmsとＭpeを、−１８０deg＜・・・≦１８０degの範囲のパラメータで表す。

Ｍms≦０かつ０≦Ｍpeの場合（Ｓ１４−３で「Ｙｅｓ」）、可視できる観測範囲に、ペリジが含まれるので、制御部１２は、式(4-2)により各軌道のペリジ高度を算出する。

さらに、制御部１２は、式(4-3)により、Ｒminを以下の値で置き換える（Ｓ１４−４）。
Ｒmin＝MIN（Ｐ0、Ｐp、Ｐm） (4-3)

一方、Ｍms＞０かつ０＞Ｍpeの場合（Ｓ１４−５で「Ｙｅｓ」）、可視できる観測範囲に、アポジが含まれるので、制御部１２は、式(4-4)により、各軌道のアポジ高度を算出する。

さらに、制御部１２は、式(4-5)により、Ｒmaxを以下の値で置き換える（Ｓ１４−６）。
Ｒmax＝MAX（Ａ0、Ａp、Ａm） (4-5)

さらに、ＴＬＥ軌道誤差テーブル２４の設定値およびマッチング誤差のマージンを考慮すると、制御部１２は、式(4-6)を用いて、予測ＳＯ候補の地心距離の存在範囲の上限および下限を計算する（Ｓ１４−７）。

次に、制御部１２は、予測ＳＯ候補の位置ベクトル方向角度の存在範囲を設定する（Ｓ１４−８）。予測ＳＯ候補の位置ベクトル方向は、２つのベクトルＲmsとＲpeに挟まれた範囲にあると考えられる。２つのベクトルＲms，Ｒpeが成す面に対して、方向誤差マージンΔθC、ΔθLを加えることにより、図１３のような位置ベクトル方向の予測存在領域を定義する。また、この予測存在領域の４隅の方向単位ベクトルをP1、P2、Q1、Q2と定義する。

以下、この予測存在領域の４隅の方向単位ベクトルP1、P2、Q1、Q2の計算方法を示す。まず、制御部１２は、ＲmsとＲpeとから、基準となる方向単位ベクトルを計算する。

制御部１２は、式(4-8)に示すように、ＴＬＥ軌道誤差の位置誤差進行方向成分と、マッチング誤差マージンにより、方向誤差マージンΔθLを計算する。
ΔθL＝ΔL／Ｒmin［rad］＋Δθobs (4-8)
ここで、Ｒminは(4-1)〜(4-3)式から計算される最小地心距離である。

制御部１２は、式(4-9)に示すように、ＴＬＥ軌道誤差の位置誤差クロストラック方向成分、速度誤差クロストラック方向成分による軌道傾斜角誤差、およびマッチング誤差マージンにより、方向誤差マージンΔθCを計算する。
ΔθC＝ΔC／Ｒmin［rad］＋ΔVC／Ｖ[rad]＋Δθobs (4-9)
Ｖ＝√（μ／ａ0s）： 平均運動

ここで、基準方向単位ベクトルUS、UEを用いて図１４のような座標系ＸＹＺsを定義する（図８の地表面座標系とは異なることに注意する。）。
・Ｘs方向： US方向
・Ｚs方向： USとUEが成す面に垂直
・Ｙs方向： Ｚs、Ｘsに対して右手系で定義

制御部１２は、各軸方向の単位ベクトルを以下で計算する。
XS＝US
ZS＝US×UE
YS＝ZS×XS

基準方向単位ベクトルUSがＸＹＺs座標に変換されると、ベクトルP1は、USをＹ軸回りに−ΔθC、次にＺ軸回りに−ΔθL回転させて得られる。ベクトルQ1は、Ｙ軸回り回転を＋ΔθCとすることにより得られる。

同様に、視線ベクトルUEをＸs方向として座標系を定義することにより、ベクトルP2、Q2が得られる。

まとめると、基準方向単位ベクトルで定義する座標系とベクトルの回転方向の関係は、以下の表のようになる。
ベクトルP1、P2、P3、P4を計算する式は、以下のようになる。なお、制御部１２は、式(4-10)により、予めZSを計算しておく。
ZS＝US×UE (4-10)
ここで、Ry（θ）、Rz（θ）は、Ｙ軸回り、Ｚ軸回りの「ベクトルの回転行列」（座標変換行列ではない）で、以下の式で表される。
Φは、赤道面座標からＸＹＺs座標への変換行列で、以下の式で表される。

次に、制御部１２は、観測位置ベクトルファイル２７から読み出した対象ＳＯの観測位置ベクトルが、予測ＳＯ候補の予測存在領域に包含されるか否かの判定を行う（Ｓ１４−９〜Ｓ１４−１６）。まず、制御部１２は、予測ＳＯ候補の予測存在領域に包含されている対象ＳＯの個数Ｎincを０で初期化する（Ｓ１４−９）。

制御部１２は、観測位置ベクトルファイル２７から対象ＳＯの位置ベクトルＲobsを１つずつ読込む。制御部１２は、式(4-15)より、位置ベクトルＲobsの地心距離Ｒobsと、位置ベクトルＲobsの方向単位ベクトルUを計算する（Ｓ１４−１０）。

制御部１２は、式(4-16)により地心距離Ｒobsが予測ＳＯ候補の予測存在範囲内に存在するか否かを判定する（Ｓ１４−１１）。
ＲＣmin≦Ｒobs かつ Ｒobs≦ＲＣmax (4-16)

地心距離Ｒobsが予測ＳＯ候補の予測存在範囲内に存在すると判定された場合（Ｓ１４−１１で「Ｙｅｓ」）、制御部１２は、位置ベクトルＲobsの方向が、予測ＳＯ候補の予測存在範囲の境界ベクトル（P1、Q1、Q2、P2）が作る範囲内に存在するか否かを判定する（Ｓ１４−１２）。ここで、制御部１２は、図１５（Ａ）に示すように、境界ベクトル（P1、Q1、Q2、P2）をそれぞれ、端点方向単位ベクトル（U1、U2、U3、U4）と置く。

制御部１２は、式(4-17)により、２つのベクトル（U1、U2）、（U2、U3）、（U3、U4）、（U4、U1）のそれぞれが成す面の法線ベクトルを求める。
図１５（Ｂ）では、一例として、判定対象ベクトルＵに対して、ベクトルU1，U2の成す面（P1−Q1）の法線ベクトルH12と、ベクトルU3，U4の成す面（P2−Q2）の法線ベクトルH34とが示されている。

位置ベクトルＲobsの方向単位ベクトルＵが、境界ベクトル（P1、Q1、Q2、P2）によって形成される範囲内に存在する条件は、図１５より、式(4-18)で表される。

位置ベクトルＲobsの方向単位ベクトル方向が境界ベクトル（P1、Q1、Q2、P2）によって形成される範囲内に存在すると判定された場合（Ｓ１４−１２で「Ｙｅｓ」）、制御部１２は、当該観測位置ベクトルは予測存在範囲に包含されると判定する。この場合、制御部１２は、対象ＳＯの個数Ｎincを＋１カウントする。

観測位置ベクトルファイル２７に格納された観測データ数分、Ｓ１４−１０〜Ｓ１４−１３の処理を繰り返す。

Ｎincのカウントが終了し、Ｎinc＞０の場合、制御部１２は、予測ＳＯ候補のカタログ識別名、Ｎinc、当該カタログＳＯの予測存在領域情報をワークファイルに出力する（Ｓ１４−１４）。当該カタログＳＯの予測存在領域情報は、例えば、ＲＣmin、ＲＣmax、端点方向単位ベクトル（U1、U2、U3、U4）等に関する情報を含む。

制御部１２は、Ｓ１４−１〜Ｓ１４−１４の処理をフィルタ２後予測位置データファイル２８に格納されている全ての予測ＳＯ候補について実施する。

制御部１２は、ワークファイルのデータをＮincについて降順にソートし（Ｓ１４−１５）、Ｎincの大きい順にカタログ名とＮinc、および予測存在領域情報を最終結果として同定結果ファイル２９に出力する。その結果、Ｎincの最も大きいカタログ名が同定結果ファイル２９の最上位に出力される。これにより、制御部１２は、対象ＳＯは、Ｎincの最も大きいカタログ名で特定されるＳＯであると同定する（Ｓ１４−１６）。なお、本実施形態では、制御部１２は、ワークファイルのデータをＮincについて降順にソートしたが、昇順でソートしてもよい。また、制御部１２は、ワークファイルのデータから、Ｎincの最も大きいカタログ名とＮinc、および予測存在領域情報を最終結果として同定結果ファイル２９に出力してもよい。

図１６は、本実施形態における同定結果ファイルの一例を示す。図１６（Ａ）では、対象ＳＯについて、４６個の観測データ（AZ/ELおよびレンジ）が取得されている。その４６個の観測データのうち３６個がカタログ番号’00877 640538’の宇宙物体の予測存在領域に包含されたことが示されている。これにより、その対象ＳＯは、カタログ番号’00877 640538’の宇宙物体であると同定される。

図１６（Ｂ）では、対象ＳＯについて、４０個の観測データ（AZ/ELおよびレンジ）が取得されている。その４０個の観測データのうち４０個がカタログ番号’37782 11043B’の宇宙物体の予測存在領域に包含されたことが示されている。また、その４０個の観測データのうち２個がカタログ番号’03611 66025E’の宇宙物体の予測存在領域に包含されたことが示されている。この場合、その対象ＳＯは、の予測存在領域に包含された数の多いカタログ番号’37782 11043B’の宇宙物体であると同定される。

図１６（Ｃ）では、対象ＳＯについて、５６個の観測データ（AZ/ELおよびレンジ）が取得されたが、予測存在領域とマッチする宇宙物体は見つからなかったことが示されている。この場合、その対象ＳＯは、ＴＬＥが公開されていない宇宙物体の可能性がある。

次に、制御部１２は、図４に示すように、同定結果の表示処理（Ｓ１５）を行う。Ｓ１５の処理について図１７を用いて説明する。

図１７は、本実施形態における同定結果を表示する画面の一例を示す。制御部１２は、Ｓ１４で作成した同定結果ファイル２９の内容を画面に表示する。このとき、制御部１２は、同定結果ファイル２９の予測存在領域情報に基づいて、対象ＳＯの観測値と、予測ＳＯ候補の予測地心距離のマッチング状況を、予測ＳＯ候補毎に色分け等を行って異なる表示形態で描画する。

図１７の例では、図１６（Ｂ）の同定結果ファイル２９の内容に基づいて、対象ＳＯの観測値と、予測ＳＯ候補の予測地心距離のマッチング状況を、予測ＳＯ候補毎に、色分けして表示されている。このとき、予測存在領域は、同定結果ファイル２９の予測存在領域情報に基づいて、描画されている。

なお、制御部１２は、対象ＳＯがＮincの最も大きいカタログ名で特定されるＳＯであると同定された場合、そのことを強調表示してもよい。

図１８は、本実施形態におけるプログラムを実行するコンピュータのハードウェア環境の構成ブロック図の一例である。コンピュータ４０は、宇宙物体同定装置１、１１として機能する。コンピュータ４０は、ＣＰＵ４２、ＲＯＭ４３、ＲＡＭ４６、通信Ｉ／Ｆ４４、記憶装置４７、出力Ｉ／Ｆ４１、入力Ｉ／Ｆ４５、読み取り装置４８、バス４９、出力機器５１、入力機器５２によって構成されている。

ここで、ＣＰＵは、中央演算装置を示す。ＲＯＭは、リードオンリメモリを示す。ＲＡＭは、ランダムアクセスメモリを示す。Ｉ／Ｆは、インターフェースを示す。バス４９には、ＣＰＵ４２、ＲＯＭ４３、ＲＡＭ４６、通信Ｉ／Ｆ４４、記憶装置４７、出力Ｉ／Ｆ４１、入力Ｉ／Ｆ４５、及び読み取り装置４８が接続されている。読み取り装置４８は、可搬型記録媒体を読み出す装置である。出力機器５１は、出力Ｉ／Ｆ４１に接続されている。入力機器５２は、入力Ｉ／Ｆ４５に接続にされている。

記憶装置４７としては、ハードディスク、フラッシュメモリ、磁気ディスクなど様々な形式の記憶装置を使用することができる。記憶装置４７またはＲＯＭ４３には、ＣＰＵ４２を取得部２、算出部３、同定部４として機能させる本実施形態に係るプログラム、具体的には、観測設備性能によるフィルタリング処理、幾何学的可視判定によるフィルタリング処理、観測値と予測値のマッチング処理を行うプログラムが格納されている。また、記憶装置４７またはＲＯＭ４３には、カタログＤＢ２１、観測データ２２、観測設備条件情報２３、ＴＬＥ軌道誤差テーブル２４、プログラム実行条件情報２５、出力データ３０、ワークファイル等が格納される。ＲＡＭ４６には、情報が一時的に記憶される。

ＣＰＵ４２は、記憶装置４７またはＲＯＭ４３から本実施形態に係るプログラムを読み出し、当該プログラムを実行する。

上記実施形態で説明した処理を実現するプログラムは、プログラム提供者側から通信ネットワーク５０、および通信Ｉ／Ｆ４４を介して、例えば記憶装置４７に格納されてもよい。また、上記実施形態で説明した処理を実現するプログラムは、市販され、流通している可搬型記憶媒体に格納されていてもよい。この場合、この可搬型記憶媒体は読み取り装置４８にセットされて、ＣＰＵ４２によってそのプログラムが読み出されて、実行されてもよい。可搬型記憶媒体としてはＣＤ−ＲＯＭ、フレキシブルディスク、光ディスク、光磁気ディスク、ＩＣカード、ＵＳＢメモリ装置など様々な形式の記憶媒体を使用することができる。このような記憶媒体に格納されたプログラムが読み取り装置４８によって読み取られる。

また、入力機器５２には、キーボード、マウス、電子カメラ、ウェブカメラ、マイク、スキャナ、センサ、タブレットなどを用いることが可能である。また、出力機器５１には、ディスプレイ、プリンタ、スピーカなどを用いることが可能である。また、ネットワーク５０は、インターネット、ＬＡＮ、ＷＡＮ、専用線、有線、無線等の通信網であってよい。

なお、本発明は、以上に述べた実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲内で種々の構成または実施形態を取ることができる。

１ 宇宙物体同定装置
２ 取得部
３ 算出部
４ 同定部
６ 格納部
７ 軌道情報
８ 誤差情報
１１ 宇宙物体同定装置
１２ 制御部
２０ 記憶装置
２１ カタログＤＢ
２２ 観測データ
２３ 観測設備条件情報
２４ ＴＬＥ軌道誤差テーブル
２５ プログラム実行条件情報
２６ フィルタ１後カタログデータ
２７ 観測位置ベクトルファイル
２８ フィルタ２後予測位置データファイル
２９ 同定結果ファイル
３０ 出力データ

Claims (5)

1. コンピュータが、
   所定地点で観測された、宇宙空間に存在する対象物体の位置に関する時系列の観測情報を取得し、
   宇宙空間に存在する物体であってそれぞれが識別された該物体の軌道に関する軌道情報と、該軌道の誤差に関する誤差情報とを格納する格納部から取得した前記軌道情報と前記誤差情報とに基づいて、前記識別された物体の軌道予測における進行方向の位置のマイナス側の位置誤差が最大となる軌道伝播による、前記観測の開始時刻における前記識別された物体の第１位置ベクトルと、前記識別された物体の軌道予測における進行方向の位置のプラス側の位置誤差が最大となる軌道伝播による、前記観測の終了時刻における前記識別された物体の第２位置ベクトルとを算出し、
   前記第１位置ベクトルと、前記第２位置ベクトルとに基づいて、前記観測情報が観測された観測時間内で、前記識別された物体が軌道伝播する軌道伝播範囲を算出し、
   前記観測時間内における、前記識別された物体毎の前記軌道伝播範囲に含まれる前記対象物体の数に基づいて、前記対象物体を同定する
   ことを特徴とする宇宙物体同定方法。
2. 前記軌道伝播範囲の算出において、
   前記誤差情報に基づいて、前記第１位置ベクトルと前記第２位置ベクトルが成す面の傾斜誤差を算出し、
   前記第１位置ベクトルと、前記第２位置ベクトルと、前記傾斜誤差とに基づいて、前記識別された物体の前記軌道伝播範囲を算出する
   ことを特徴とする請求項１に記載の宇宙物体同定方法。
3. 前記同定において、前記対象物体は、複数の前記識別された物体のそれぞれの前記軌道伝播範囲において、前記数が最も多い前記軌道伝播範囲に対応する該識別された物体であると同定する
   ことを特徴とする請求項１または２に記載の宇宙物体同定方法。
4. 所定地点で観測された、宇宙空間に存在する対象物体の位置に関する時系列の観測情報を取得する取得部と、
   宇宙空間に存在する物体であってそれぞれが識別された該物体の軌道に関する軌道情報と、該軌道の誤差に関する誤差情報とを格納する格納部と、
   前記格納部から取得した前記軌道情報と前記誤差情報とに基づいて、前記識別された物体の軌道予測における進行方向の位置のマイナス側の位置誤差が最大となる軌道伝播による、前記観測の開始時刻における前記識別された物体の第１位置ベクトルと、前記識別された物体の軌道予測における進行方向の位置のプラス側の位置誤差が最大となる軌道伝播による、前記観測の終了時刻における前記識別された物体の第２位置ベクトルとを算出し、前記第１位置ベクトルと、前記第２位置ベクトルとに基づいて、前記観測情報が観測された観測時間内で、前記識別された物体が軌道伝播する軌道伝播範囲を算出する算出部と、
   前記観測時間内における、前記識別された物体毎の前記軌道伝播範囲に含まれる前記対象物体の数に基づいて、前記対象物体を同定する同定部と、
   を備えることを特徴とする宇宙物体同定装置。
5. コンピュータに、
   所定地点で観測された、宇宙空間に存在する対象物体の位置に関する時系列の観測情報を取得し、
   宇宙空間に存在する物体であってそれぞれが識別された該物体の軌道に関する軌道情報と、該軌道の誤差に関する誤差情報とを格納する格納部から取得した前記軌道情報と前記誤差情報とに基づいて、前記識別された物体の軌道予測における進行方向の位置のマイナス側の位置誤差が最大となる軌道伝播による、前記観測の開始時刻における前記識別された物体の第１位置ベクトルと、前記識別された物体の軌道予測における進行方向の位置のプラス側の位置誤差が最大となる軌道伝播による、前記観測の終了時刻における前記識別された物体の第２位置ベクトルとを算出し、
   前記第１位置ベクトルと、前記第２位置ベクトルとに基づいて、前記観測情報が観測された観測時間内で、前記識別された物体が軌道伝播する軌道伝播範囲を算出し、
   前記観測時間内における、前記識別された物体毎の前記軌道伝播範囲に含まれる前記対象物体の数に基づいて、前記対象物体を同定する
   処理を実行させることを特徴とする宇宙物体同定プログラム。