(57)【特許請求の範囲】
【請求項１】
衛星の軌道を、前記衛星を観測することにより得られた物理量から計算する軌道計算装置において、
前記物理量に基づいて前記衛星の軌道を第１の手法によって計算する第１の軌道計算手段と、
前記物理量に基づいて前記衛星の軌道を前記第１の手法とは異なる第２の手法によって計算する第２の軌道計算手段と、
前記第１の軌道計算手段によって得られた結果を、前記第２の軌道計算手段に供給する供給手段と、
前記第１および第２の軌道計算手段によって得られた計算結果を同時に表示出力する出力手段と、
を有することを特徴とする軌道計算装置。
【請求項２】
前記第１の軌道計算手段は、バッチ推定により軌道を計算し、
前記第２の軌道計算手段は、前記供給手段から供給された情報を参照してリアルタイム推定により軌道を計算することを特徴とする請求項１記載の軌道計算装置。
【請求項３】
前記第１の軌道計算手段は、軌道の予測誤差が所定の閾値を越えた場合に前記バッチ推定処理による計算を行うことを特徴とする請求項２記載の軌道計算装置。
【請求項４】
衛星の軌道を、衛星を観測することにより得られた物理量から計算する軌道計算方法において、
第１の軌道計算手段が、前記物理量に基づいて前記衛星の軌道を第１の手法によって計算する第１の軌道計算ステップと、
第２の軌道計算手段が、前記物理量に基づいて前記衛星の軌道を前記第１の手法とは異なる第２の手法によって計算する第２の軌道計算ステップと、
供給手段が、前記第１の軌道計算ステップによって得られた結果を、前記第２の軌道計算ステップに供給する供給ステップと、
出力手段が、前記第１および第２の軌道計算ステップによって得られた計算結果を同時に表示出力する出力ステップと、
を有することを特徴とする軌道計算方法。
【請求項５】
衛星の軌道を、前記衛星を観測することにより得られた物理量から計算する処理をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体において、
コンピュータを、
前記物理量に基づいて前記衛星の軌道を第１の手法によって計算する第１の軌道計算手段、
前記物理量に基づいて前記衛星の軌道を前記第１の手法とは異なる第２の手法によって計算する第２の軌道計算手段、
前記第１の軌道計算手段によって得られた結果を、前記第２の軌道計算手段に供給する供給手段、
前記第１および第２の軌道計算手段によって得られた計算結果を同時に表示出力する出力手段、
として機能させるプログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は軌道計算装置および軌道計算方法に関し、特に、宇宙空間を所定の軌道を描きながら飛翔する衛星の軌道を、衛星を観測することにより得られた物理量から計算する軌道計算装置および軌道計算方法に関する。
【０００２】
【従来の技術】
人工衛星等のように宇宙空間に軌道を描きながら飛翔する物体の状態（ある時刻における位置ベクトルおよび速度ベクトル）を観測に基づいて計算する方法には、バッチ推定方式とリアルタイム推定方式とがある。
【０００３】
バッチ推定方式は、観測データを蓄積し、一括処理する方式である。このようなバッチ推定方式では、計算モデルの精密化が可能であるとともに、大量の観測データを扱うことから、精度の高い結果を得ることができるという特徴がある。
【０００４】
一方、リアルタイム推定方式は、観測データを逐次的に処理する方式であり、観測データが取得される度に処理を行うことから、迅速に計算結果を得ることができるという特徴がある。
【０００５】
【発明が解決しようとする課題】
しかしながら、バッチ推定方式では、一定時間毎に処理を行う関係上、即時性がないという問題点があった。
【０００６】
また、リアルタイム推定方式では、短い時間間隔（例えば、３秒）で取り込まれる観測データに対して逐次処理を施す必要があることから、モデルを簡略化せざるを得ず、結果的に精度の低下を招くという問題点があった。また、リアルタイム推定方式には過収束という現象が発生しやすく、バッチ推定方式のように高精度の結果を得ることができないという問題点があった。ここで、過収束とは、推定の確からしさを示す推定結果共分散行列を見ただけでは十分精度が高いように見えて、実は誤った結果に収束しているという現象である。
【０００７】
本発明はこのような点に鑑みてなされたものであり、高精度でしかも即時的な軌道の推定を行うことが可能な軌道計算システムおよびそのようなシステムを構成する軌道計算装置および軌道計算方法を提供することを目的とする。
【０００８】
【課題を解決するための手段】
本発明では上記課題を解決するために、図１に示す、人工衛星１の軌道を、前記人工衛星１を観測することにより得られた物理量から計算する軌道計算装置（計算装置４）において、前記物理量に基づいて前記人工衛星１の軌道を第１の手法によって計算する第１の軌道計算手段４ｂと、前記物理量に基づいて前記人工衛星１の軌道を前記第１の手法とは異なる第２の手法によって計算する第２の軌道計算手段４ｃと、前記第１の軌道計算手段４ｂによって得られた結果を、前記第２の軌道計算手段４ｃに供給する供給手段４ｄと、前記第１および第２の軌道計算手段によって得られた計算結果を同時に表示出力する出力手段４ｅと、を有することを特徴とする軌道計算装置（計算装置４）が提供される。
【０００９】
ここで、第１の軌道計算手段４ｂは、物理量に基づいて人工衛星１の軌道を第１の手法によって計算する。第２の軌道計算手段４ｃは、物理量に基づいて人工衛星１の軌道を第１の手法とは異なる第２の手法によって計算する。供給手段４ｄは、第１の軌道計算手段４ｂによって得られた結果を、第２の軌道計算手段４ｃに供給する。出力手段４ｅは、第１および第２の軌道計算手段によって得られた計算結果を同時に表示出力する。
【００１０】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して説明する。
図１は、本発明に係る軌道計算システムの動作原理を説明する原理図である。この図において、人工衛星１は、例えば、周回衛星であり地球の周りを所定の周期で周回する。
【００１１】
パラボラアンテナ２は、人工衛星１との間に電波リンクを有しており、人工衛星１との間で情報を授受する。
観測装置３は、測定手段３ａおよび送信手段３ｂによって構成されており、人工衛星１とパラボラアンテナ２の間の電波リンクを利用して観測データを取得し、計算装置４に対して得られた観測データを送信する。
【００１２】
ここで、測定手段３ａは、人工衛星１とパラボラアンテナ２の間の電波リンクを利用して、人工衛星１の運動に係る物理量を測定する。
送信手段３ｂは、測定手段３ａの測定結果のデータ（観測データ）を、計算装置４に対して送信する。
【００１３】
計算装置４は、受信手段４ａ、第１の軌道計算手段４ｂ、第２の軌道計算手段４ｃ、供給手段４ｄ、および、出力手段４ｅによって構成されており、観測装置３から供給された観測データから、人工衛星１の軌道を計算して表示装置５に表示する。
【００１４】
ここで、受信手段４ａは、観測装置３から送信された観測データを受信する。第１の軌道計算手段４ｂは、第１の計算手法であるバッチ推定処理により人工衛星１の軌道を計算する。
【００１５】
第２の軌道計算手段４ｃは、第２の計算手法であるリアルタイム推定処理により人工衛星１の軌道を計算する。
供給手段４ｄは、第１の軌道計算手段４ｂによって得られた結果を、第２の軌道計算手段４ｃに供給する。
【００１６】
出力手段４ｅは、第１または第２の軌道計算手段４ｂ，４ｃによって得られた計算結果を表示装置５に対して出力する。
表示装置５は、例えば、ＣＲＴ（Ｃａｔｈｏｄｅ  Ｒａｙ  Ｔｕｂｅ）モニタによって構成されており、計算装置４から供給された計算結果を表示出力する。
【００１７】
次に、以上の原理図の動作について説明する。
観測装置３の測定手段３ａは、人工衛星１とパラボラアンテナ２の間の電波リンクを利用して、電波の伝送時間、電波の周波数の変化量、および、パラボラアンテナ２の仰角等を観測データとして取得し、送信手段３ｂに供給する。
【００１８】
送信手段３ｂは、測定手段３ａから供給された観測データを、計算装置４に対して送信する。
計算装置４の受信手段４ａは、観測装置３から送信されてきた観測データを受信し、第１の軌道計算手段４ｂと第２の軌道計算手段４ｃにそれぞれ供給する。
【００１９】
第１の軌道計算手段４ｂは、観測データを一定量だけ蓄積した後、蓄積された観測データと、人工衛星１の運動方程式から理論的に導出された計算データとの誤差が最小となるように人工衛星１の軌道を計算する。従って、第１の軌道計算手段４ｂからは、所定の時間毎に計算結果が出力される。なお、第１の軌道計算手段４ｂによって計算された結果は、供給手段４ｄと出力手段４ｅにそれぞれ供給される。
【００２０】
供給手段４ｄは、第１の軌道計算手段４ｂから所定の時間毎に出力される計算結果を第２の軌道計算手段４ｃに対して供給する。
第２の軌道計算手段４ｃは、受信手段４ａによって受信された観測データと、供給手段４ｄによって供給された第１の軌道計算手段４ｂの計算結果とを参照して、リアルタイム推定処理により人工衛星１の軌道を算出して出力する。なお、第２の軌道計算手段４ｃは、観測データに対して逐次的に処理を行うので、観測データが受信される周期と同一の周期で計算結果を出力する。
【００２１】
ここで、第２の軌道計算手段４ｃは、観測データを逐次的に処理することを要求されるので、第１の軌道計算手段４ｂよりも簡易な方法を用いて計算を行っており、計算結果が実際とは異なる解に収束する「過収束」が発生する場合がある。過収束に陥った場合には、推定された結果と、実際の測定結果がよく一致することから、そのような状態に陥っていることを知り得ない場合が多く、万一、知り得たとしても過収束の状態から離脱することが困難であった。
【００２２】
本発明においては、前述のように、第２の軌道計算手段４ｃが、第１の軌道計算手段４ｂの計算結果を参照して軌道計算を行うようにしたので、過収束が発生することを防止できる。即ち、本発明では、即時性は高いが精度が低い第２の計算手法（リアルタイム推定）を、即時性は低いが精度が高い第１の計算手法（バッチ推定）によって算出された結果により補間し、より精度が高い計算結果を得ることが可能となる。
【００２３】
このようにして、第１の軌道計算手段４ｂおよび第２の軌道計算手段４ｃの計算結果は、出力手段４ｅを介して表示装置５に表示出力される。このとき、第２の軌道計算手段４ｃによって算出された結果のみでなく、第１の軌道計算手段４ｂによる計算結果を同一の画面に重ねて表示することにより、両者の誤差が一目瞭然となるので、計算結果に異常が生じた場合には直ちにこれを知ることが可能となる。
【００２４】
次に、図２を参照して、本発明の実施の形態の構成例について説明する。
この図において、人工衛星１０は、例えば、周回衛星であり、地上局に向けて測位信号を送信している。
【００２５】
パラボラアンテナ２０，２１は、人工衛星１０から送信された測位信号を地球上の異なる位置で受信し、観測装置３０に対して供給する。
観測装置３０は、観測データ取得部３０ａおよび送信部３０ｂによって構成されており、人工衛星１０から送信されてくる測位信号を利用して人工衛星１０の観測データを生成し、ネットワーク４０を介して計算装置５０に供給する。
【００２６】
ここで、観測データ取得部３０ａは、前述のように、人工衛星１０から送信されてくる測位信号を利用して電波の伝搬時間、電波の周波数の変化量、および、パラボラアンテナ２０の仰角等を観測データとして取得する。
【００２７】
送信部３０ｂは、観測データ取得部３０ａによって取得された観測データを、ネットワーク４０を介して計算装置５０に送信する。
ネットワーク４０は、例えば、ディジタル回線等によって構成されており、観測装置３０によって取得された観測データを計算装置５０まで伝送する。
【００２８】
計算装置５０は、受信部５０ａ、記憶部５０ｂ、バッチ推定処理部５０ｃ、基準情報生成部５０ｄ、リアルタイム推定処理部５０ｅ、表示処理部５０ｆ、および、制御部５０ｇによって構成されており、観測装置３０から送信された観測データから、人工衛星１０の軌道を計算して表示装置６０に結果を表示する。
【００２９】
ここで、受信部５０ａは、ネットワーク４０を介して送信されてきた観測データを受信する。
記憶部５０ｂは、受信部５０ａによって受信された観測データを所定量だけ記憶し、バッチ推定処理部５０ｃに対して供給する。
【００３０】
バッチ推定処理部５０ｃは、人工衛星１０が軌道上で受ける様々な外乱（摂動）を考慮に入れたニュートンの運動方程式に対して、各々の観測時刻において人工衛星１０が受ける加速度をモデル化した運動モデルを適用することにより人工衛星１０の理論的な軌道情報（人工衛星の位置ベクトルおよび速度ベクトル等）を算出する。そして、得られた軌道情報に対して、その位置において観測データが受ける誤差を示す観測モデルを適用することにより、理論的な観測データを導出する。最後に、計算によって得られた理論的な観測データと、実際の観測データとを比較し、これらの差の二乗和が最小となるように人工衛星１０軌道を算定する。
【００３１】
基準情報生成部５０ｄは、バッチ推定処理部５０ｃによって算出された軌道情報から、各観測時刻における人工衛星１０の軌道その他を示す基準情報を生成し、リアルタイム推定処理部５０ｅと表示処理部５０ｆにそれぞれ供給する。ここで、基準情報とは、以下のデータを示す。
（Ａ）人工衛星１０の赤道面座標系における位置ベクトルおよび速度ベクトル
（Ｂ）地上局（パラボラアンテナ２０，２１）の赤道面座標系における位置ベクトル
（Ｃ）人工衛星１０の赤道面座標系における状態遷移行列
（Ｄ）観測データ取得可能時間帯情報
なお、状態遷移行列は、運動方程式に従うパラメータＸ（ｔ）があった場合に、以下の式で表すことができる。
【００３２】
【数１】
Φ（ｔ’，ｔ）＝∂Ｘ（ｔ’）／∂Ｘ（ｔ）  ・・・（１）
いまの場合、Ｘ（ｔ）は、人工衛星１０の位置ベクトルおよび速度ベクトルであり、Φの物理的な意味としては、時刻ｔで位置ベクトルまたは速度ベクトルに微小変化が生じた場合に、ｔ’における位置ベクトルまたは速度ベクトルがどれだけ変化するかを示している。換言すると、ΦはＸの誤差を伝播する式である。
【００３３】
また、観測データ取得可能時間帯情報とは、人工衛星１０が観測可能な時間帯を示す情報である。即ち、本実施の形態では、人工衛星１０は周回衛星であるので、この人工衛星１０が地上から観測可能な時間帯を示す情報である。
【００３４】
リアルタイム推定処理部５０ｅは、カルマンフィルタを用いたリアル推定処理により人工衛星１０の軌道を計算する。このカルマンフィルタによる処理は、観測更新と時間更新と呼ばれる２つの処理によって成り立っている。観測更新では、取得された観測データと、カルマンフィルタによって予測された観測量との差分に応じたフィードバックをかけて人工衛星１０の軌道情報を修正する。更新された軌道情報は、時間更新処理によって次の観測データが取得される時刻まで伝播され、そこでまた新たな更新処理を受けることになる。このような処理が繰り返されることにより、人工衛星１０の軌道情報が逐次推定される。なお、本実施の形態では、リアルタイム推定処理部５０ｅは、バッチ推定処理部５０ｃから供給された推定値を初期値とするとともに、基準情報生成部５０ｄから供給された基準情報を参照して前述の観測更新と時間更新とを行うことにより推定精度を向上させている。
【００３５】
表示処理部５０ｆは、基準情報生成部５０ｄと、リアルタイム推定処理部５０ｅとによって算出された軌道情報を表示装置６０に供給して表示させる。
制御部５０ｇは、装置の各部を制御する。
【００３６】
表示装置６０は、例えば、ＣＲＴモニタ等によって構成されており、表示処理部５０ｆから供給された画像信号を表示出力する。
次に、以上の実施の形態の動作について説明する。
【００３７】
人工衛星１０から送信された測位信号は、地上局であるパラボラアンテナ２０，２１によって受信され、観測装置３０に対して供給される。
観測装置３０は、観測データ取得部３０ａによって、パラボラアンテナ２０，２１から供給された測位信号を観測データとして取得し、送信部３０ｂに供給する。
【００３８】
送信部３０ｂは、観測データ取得部３０ａから供給された観測データを、ネットワーク４０を介して計算装置５０に対して送信する。
計算装置５０では、観測装置３０から送信された観測データを、受信部５０ａにより受信する。
【００３９】
受信部５０ａによって受信された観測データは、記憶部５０ｂに供給されてそこに記憶され、一定量の観測データが蓄積された場合には、バッチ推定処理部５０ｃに対して蓄積された観測データが供給される。
【００４０】
バッチ推定処理部５０ｃは、運動方程式に対して運動モデルを適用して人工衛星１０の軌道情報の理論値を算出する。次に、得られた軌道情報の理論値に対して観測モデルを適用し、各観測データが取得された時刻における観測データの理論値を算出する。続いて、観測データの理論値と、実際の観測データとの差分の二乗和が最小となる軌道情報を算出し、これを人工衛星１０の最も確からしい軌道情報と推定する。
【００４１】
このようにして得られた軌道情報は、基準情報生成部５０ｄに供給されるとともに、初期値情報としてリアルタイム推定処理部５０ｅに供給される。
基準情報生成部５０ｄでは、バッチ推定処理部５０ｃによる処理結果を取得し、これに基づいて前述の（Ａ）〜（Ｄ）からなる基準情報を生成し、リアルタイム推定処理部５０ｅに供給する。
【００４２】
リアルタイム推定処理部５０ｅは、バッチ推定処理部５０ｃから供給された初期値を参照して、リアルタイム推定処理を実行する。即ち、リアルタイム推定処理部５０ｅは、バッチ推定処理部５０ｃから供給された初期値を用いてカルマンフィルタを設定する。そして、このようにして設定されたカルマンフィルタに対して、基準情報の状態遷移行列を用いて時間更新処理を行うとともに、軌道情報と状態遷移行列とを用いて観測更新処理を行う。
【００４３】
図３は、リアルタイム推定処理とバッチ推定処理の関係を示す図である。この図に示すように、観測データ（図３（Ａ）参照）の入力が開始されると、これらのデータは記憶部５０ｂに蓄積される（図３（Ｃ）参照）。そして、一定量の観測データが蓄積されると、バッチ推定処理（図３（Ｄ）参照）が起動され、初期値と基準情報とがリアルタイム推定処理部５０ｅに供給される。その結果、リアルタイム推定処理部５０ｅは、これらの情報と、観測データとを参照して、リアルタイム推定処理を実行し、人工衛星１０の軌道情報を算出する。
【００４４】
なお、リアルタイム推定処理部５０ｅは、図４に示すように、バッチ推定処理によって導出された基準情報の軌道情報を「基準軌道」とし、各時刻における基準軌道と、リアルタイム推定による推定値とのずれΔｒｋを計算結果として出力する。
【００４５】
表示処理部５０ｆは、基準情報生成部５０ｄとリアルタイム推定処理部５０ｅによる推定結果のそれぞれを評価対象パラメータに変換し、表示装置６０に供給して表示させる。なお、評価対象パラメータとは、推定によって得られた軌道情報である速度ベクトルと位置ベクトルはそれぞれ３次元のベクトルであることから、これらを直感的に理解可能な値（例えば、人工衛星１０の直下点経度）に変換したものをいう。
【００４６】
図５は、以上の処理の結果、表示装置６０に表示される画面の表示例である。この表示例では、横軸は時間を示し、また、縦軸は直下点経度を示している。また、横軸において三角形「△」が付されている位置では、基準情報の更新がなされた（バッチ推定処理による処理結果が出力された）ことを示している。更に、太線はバッチ推定処理の結果から求められた直下点経度であり、破線はバッチ推定処理の誤差範囲を示しており、細線はリアルタイム推定処理から求められた直下点経度を示している。
【００４７】
このような表示画面を参照することにより、リアルタイム推定処理によって得られた結果と、バッチ推定処理によって得られた結果とを比較考量することが可能となるので、これらによる推定結果が乖離していない場合には、双方の推定結果がともに確からしいことを知ることができる。
【００４８】
次に、図６〜図８を参照して、図２に示す実施の形態において実行される処理の詳細について説明する。
先ず、図６は、バッチ推定処理部５０ｃにおいて実行される初期設定処理の一例を説明するフローチャートである。なお、この処理は、図３に示す最初の「バッチ推定処理」に対応している。このフローチャートが開始されると、以下の処理が実行される。
［Ｓ１］バッチ推定処理部５０ｃは、記憶部５０ｂに一定量の観測データが蓄積されたか否かを判定し、蓄積された場合にはステップＳ２に進み、それ以外の場合にはステップＳ１に戻って同様の処理を繰り返す。
［Ｓ２］バッチ推定処理部５０ｃは、蓄積された観測データを用いてバッチ推定処理を開始する。
［Ｓ３］バッチ推定処理部５０ｃは、リアルタイム推定の初期値を算出し、リアルタイム推定処理部５０ｅに対して供給する。
［Ｓ４］基準情報生成部５０ｄは、バッチ推定処理部５０ｃによる軌道の推定結果を基準軌道に設定する。
［Ｓ５］基準情報生成部５０ｄは、状態遷移行列を所定の将来時刻（次回のバッチ推定による更新予定時刻）まで伝播する。
【００４９】
この際、ある一定間隔の時刻に対して伝播を行って複数の状態遷移行列を予め生成しておき、リアルタイム推定処理において実際に用いる時刻に対しては、これらの情報を適宜補間することにより、任意の時刻の状態遷移行列を得る。
［Ｓ６］基準情報生成部５０ｄは、地上局であるパラボラアンテナ２０，２１の位置情報と、可視情報を算出する。
【００５０】
ここで、パラボラアンテナ２０，２１の位置情報とは、慣性座標系である赤道面座標系におけるパラボラアンテナ２０，２１の位置と速度ベクトルを示している。また、可視情報は、前述のように人工衛星１０が測定可能な時間帯を示す情報である。基準情報生成部５０ｄは、これらの情報を生成する。
【００５１】
次に、図７を参照して、通常状態におけるバッチ推定処理（図３に示す第２番目以降のバッチ推定処理）の詳細について説明する。このフローチャートが開始されると、以下の処理が実行されることになる。
［Ｓ２０］記憶部５０ｂは、観測データを蓄積する。
［Ｓ２１］バッチ推定処理部５０ｃは、予測誤差ΔＲを算定する。
【００５２】
ここで、予測誤差ΔＲは以下の式で表される。
【００５３】
【数２】
ΔＲ＝（３ｎ／２）Δａ×Ｔ  ・・・（２）
ここで、ΔＲは、正確には人工衛星１０の進行方向における位置の誤差である。また、ｎは人工衛星１０の周回周期の逆数である平均運動を、Δａは軌道長半径誤差を、Ｔは前回のバッチ推定処理からの経過時間を示している。
［Ｓ２２］バッチ推定処理部５０ｃは、予測誤差ΔＲが運用上要求される制度ΔＲＬよりも大きい場合にはステップＳ２３に進み、それ以外の場合にはステップＳ２０に戻って同様の処理を繰り返す。
【００５４】
即ち、バッチ推定処理では、予測値とともにその予測誤差範囲が得られるが、誤差範囲が時間とともに増大するパラメータがあり、その典型例が前述の軌道長半径の誤差による人工衛星１０の軌道上進行方向の位置誤差ΔＲである。このΔＲは前述のように時間Ｔとともに増大するので、ΔＲの値が運用上要求される制度ΔＲＬを越えた場合には、新たなバッチ推定処理を行って精度の低下を補う。［Ｓ２３］バッチ推定処理部５０ｃは、バッチ推定処理を実行する。
［Ｓ２４］バッチ推定処理部５０ｃは、バッチ推定処理の初期値を算出し、リアルタイム推定処理部５０ｅに供給する。
［Ｓ２５］基準情報生成部５０ｄは、バッチ推定処理部５０ｃによる軌道の予測結果を基準軌道に設定する。
［Ｓ２６］基準情報生成部５０ｄは、状態遷移行列を所定の時刻（リアルタイム推定の開始時刻）まで伝播する。
［Ｓ２７］基準情報生成部５０ｄは、地上局であるパラボラアンテナ２０，２１の位置情報の初期値と、可視情報の初期値を算出する。
［Ｓ２８］基準情報生成部５０ｄは、基準情報を更新する。即ち、前回のバッチ推定処理によって算出された基準情報を、今回の処理によって生成された新たな基準情報により更新する。
［Ｓ２９］処理を継続する場合には、ステップＳ２０に戻って同様の処理を繰り返し、それ以外の場合には処理を終了する。
【００５５】
次に、図８を参照して、リアルタイム推定処理の詳細について説明する。この処理が開始されると、以下の処理が実行される。
［Ｓ４０］リアルタイム推定処理部５０ｅは、人工衛星１０から送信された測位信号がパラボラアンテナ２０，２１によって受信された時刻ｔを、実際の時刻ｔ’に補正する。
【００５６】
即ち、人工衛星１０から送信され、パラボラアンテナ２０，２１によって受信された測位信号には、受信時刻を示すタグが付与されるが、実際に受信されてからタグが付与される迄の間にタイムラグがあるのでこれを修正する。
［Ｓ４１］リアルタイム推定処理部５０ｅは、基準情報生成部５０ｄから状態遷移行列を取得する。
【００５７】
なお、必要に応じて補間処理等を行い、目的となる時刻における状態遷移行列を生成する。
［Ｓ４２］リアルタイム推定処理部５０ｅは、カルマンフィルタを初期値（前回の値）からｔ’まで時間更新する。
【００５８】
ここで、通常のリアルタイム処理の時間更新処理では積分により所定の時刻まで時間更新を行うが、本実施の形態では、バッチ推定処理によって生成された状態遷移行列をカルマンフィルタに対して乗算することにより時間更新を行う。従って、処理に時間を要する積分を省略することができるので、処理時間を短縮することが可能となる。また、状態遷移行列は、リアルタイム処理に比較して精度の高いバッチ推定処理によって生成されたものであるので、リアルタイム推定処理の精度を向上させることも可能となる。
［Ｓ４３］リアルタイム推定処理部５０ｅは、基準情報生成部５０ｄから地上局位置情報と、可視情報とを取得し、これらから時刻ｔ’における地上局位置情報と、時刻ｔ’における可視情報とを算出する。
【００５９】
なお、この算出方法としては、例えば、エルミート補間処理を用いることができる。
［Ｓ４４］リアルタイム推定処理部５０ｅは、ステップＳ４３において取得した地上局位置情報と可視情報とから、人工衛星１０の観測が可能であるか否かを判定し、観測可能である場合にはステップＳ４５に進み、それ以外の場合にはステップＳ４９に進む。
［Ｓ４５］リアルタイム推定処理部５０ｅは、人工衛星１０が測位信号を発信した時刻を推定する処理を実行する。
【００６０】
即ち、測位データが発信されてから、受信されるまでの間に人工衛星１０が移動するので、その移動量を考慮して測位信号が発信された時刻をライトタイムイクエーション（Ｌｉｇｈｔ  Ｔｉｍｅ  Ｅｑｕａｔｉｏｎ）により算出する。なお、その際には、基準情報生成部５０ｄによって生成された基準情報を用いる。
［Ｓ４６］リアルタイム推定処理部５０ｅは、人工衛星１０のクロックが有する誤差（クロックバイアス）を算出する。
【００６１】
なお、この誤差は、時刻を変数とする二次式によって与えられるので、この二次式を解くことにより、クロックバイアスを得る。
［Ｓ４７］リアルタイム推定処理部５０ｅは、ステップＳ４５において得られた発信時刻と、ステップＳ４６において得られたクロックバイアスとから、測位信号の実際の送信時刻を算出する。
［Ｓ４８］リアルタイム推定処理部５０ｅは、観測感度行列を算出する。
【００６２】
ここで、観測感度行列とは、観測される物理量ｈ（地上局と人工衛星間の距離、または、地上局からみた人工衛星の方向等）を推定パラメータＸ（いまの場合は人工衛星１０の位置ベクトルまたは速度ベクトル）の関数として、ｈ（Ｘ）として表したとすると、観測感度行列Ｈは以下の式により与えられる。
【００６３】
【数３】
Ｈ＝∂ｈ（Ｘ）／∂Ｘ  ・・・（３）
なお、この物理的な意味は次のように解される。即ち、Ｈが大きいということは、推定パラメータＸのわずかな変化に対して観測される物理量ｈが大きく変化することである。換言すると、Ｈが大きい観測量を使用すると、Ｘを精度よく決めることができる。逆にＨが小さい観測量ではＸの変化に対してｈの変化が鈍いのでＸが決まりにくい。
［Ｓ４９］リアルタイム推定処理部５０ｅは、基準情報生成部５０ｄから供給された基準情報に含まれている軌道情報と情報遷移行列、および、ステップＳ４８において算出された観測感度行列を用いて、カルマンフィルタの観測更新処理を実行する。
［Ｓ５０］処理を継続する場合にはステップＳ４０に戻って同様の処理を繰り返し、それ以外の場合には処理を終了する。
【００６４】
以上に説明したように、本発明の実施の形態では、バッチ推定処理とリアルタイム推定処理とを並行して実行するとともに、バッチ推定処理により得られた結果を、リアルタイム推定処理に反映させることができるので、リアルタイム推定処理の精度を向上させることが可能となる。
【００６５】
また、バッチ推定処理に基づく状態遷移行列を、リアルタイム推定処理の時間更新処理に利用するようにしたので、時間更新処理をより高精度にしかも迅速に実行することが可能となる。
【００６６】
更に、リアルタイム推定処理の実運用精度をモニタすることが可能となる。これについて次に詳述する。即ち、軌道決定精度評価に使用されるパラメータは通常以下の２種類である。
（ａ）共分散行列
（ｂ）観測データの実測値と理論値との差
ここで、リアルタイム推定処理単独では、その場で評価できるのは（ａ）と（ｂ）のトレンドだけである。しかし、（ａ）は過収束の状態になった場合はもとより、絶対精度を評価するパラメータにはなり得ない。また、（ｂ）については、実測値と理論値の差は小さいにも拘わらず、実際には誤差が含まれていることがあり得る。
【００６７】
本来、「実運用的精度」として要求されるのは、人工衛星の位置が真値からどれだけずれているかということであるので、リアルタイム推定処理とバッチ推定処理の結果を比較することができる本実施の形態における方法が最もふさわしいと言える。また、リアルタイム推定処理を実施しながら、かつ、結果の実運用的精度評価を同時に実施できるので、従来のように、問題が生じた場合に事後的に解析作業を行うといった手間を省略することが可能となる。
【００６８】
なお、以上の実施の形態においては、周回衛星を観測の対象としたが、静止衛星を観測の対象とすることも可能であることはいうまでもない。
最後に、上記の処理機能は、コンピュータによって実現することができる。その場合、軌道計算装置が有すべき機能の処理内容は、コンピュータで読み取り可能な記録媒体に記録されたプログラムに記述されており、このプログラムをコンピュータで実行することにより、上記処理がコンピュータで実現される。コンピュータで読み取り可能な記録媒体としては、磁気記録装置や半導体メモリ等がある。市場へ流通させる場合には、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｋ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やフロッピーディスク等の可搬型記録媒体にプログラムを格納して流通させたり、ネットワークを介して接続されたコンピュータの記憶装置に格納しておき、ネットワークを通じて他のコンピュータに転送することもできる。コンピュータで実行する際には、コンピュータ内のハードディスク装置等にプログラムを格納しておき、メインメモリにロードして実行する。
【００６９】
【発明の効果】
以上説明したように本発明では、衛星の軌道を、衛星を観測することにより得られた物理量から計算する軌道計算装置において、物理量に基づいて衛星の軌道を第１の手法によって計算する第１の軌道計算手段と、物理量に基づいて衛星の軌道を第１の手法とは異なる第２の手法によって計算する第２の軌道計算手段と、第１の軌道計算手段によって得られた結果を、第２の軌道計算手段に供給する供給手段と、第１および第２の軌道計算手段によって得られた計算結果を同時に表示出力する出力手段と、を有するようにしたので、計算によって得られた軌道の精度を簡単に確認することが可能となる。
【図面の簡単な説明】
【図１】本発明の動作原理を説明する原理図である。
【図２】本発明の実施の形態の構成例を示す図である。
【図３】リアルタイム推定処理とバッチ推定処理の関係を示す図である。
【図４】基準軌道を利用したリアルタイム推定の概念図である。
【図５】図２に示す表示装置に表示される表示画面の一例である。
【図６】図２に示すバッチ推定処理部において実行される初期設定処理の一例を説明するためのフローチャートである。
【図７】通常状態におけるバッチ推定処理の一例を説明するためのフローチャートである。
【図８】リアルタイム推定処理の一例を説明するフローチャートである。
【符号の説明】
１  人工衛星
２  パラボラアンテナ
３  観測装置
３ａ  測定手段
３ｂ  送信手段
４  計算装置
４ａ  受信手段
４ｂ  第１の軌道計算手段
４ｃ  第２の軌道計算手段
４ｄ  供給手段
４ｅ  出力手段
５  表示装置
１０  人工衛星
２０，２１  パラボラアンテナ
３０  観測装置
３０ａ  観測データ取得手段
３０ｂ  送信手段
４０  ネットワーク
５０  計算装置
５０ａ  受信部
５０ｂ  記憶部
５０ｃ  バッチ推定処理部
５０ｄ  基準情報生成部
５０ｅ  リアルタイム推定処理部
５０ｆ  表示処理部
５０ｇ  制御部
６０  表示装置
【図１】
0003550520.tif 000002
【図２】
0003550520.tif 000003
【図３】
0003550520.tif 000004
【図４】
0003550520.tif 000005
【図５】
0003550520.tif 000006
【図６】
0003550520.tif 000007
【図７】
0003550520.tif 000008
【図８】
0003550520.tif 000009