(57)【特許請求の範囲】
【請求項１】
人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する格納部と、
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いることにより、時刻ｔ０における状態量の変化量Δｘ’０を収束判定をすることなく計算し、時刻ｔ０における基準軌道状態量ｘ’０に、計算した前記変化量Δｘ’０を加算することにより、時刻ｔ０における前記人工衛星の状態量の推定値である推定値ｘ０を計算する計算部と
を備える軌道計算装置。
【請求項２】
前記計算部は、
前記格納部に格納された前記基準軌道情報と、計算された前記状態量の前記推定値とを用いて、前記推定値の計算に用いられた前記観測量の存在する時刻よりも未来の時刻における前記人工衛星の状態量であり、新たな観測量を用いることなく計算される前記人工衛星の状態量である予測値を計算する請求項１に記載の軌道計算装置。
【請求項３】
前記計算部は、
前記基準軌道情報と、前記状態量の前記予測値とを用いて、前記観測量の存在する前記時刻よりも前記未来の時刻における前記人工衛星の前記予測値を異なる予測値に変更するための状態量の変化量であって、前記未来の時刻と前記観測量の存在する前記時刻との間の時刻における状態量の変化量である中間変化量を計算する請求項２に記載の軌道計算装置。
【請求項４】
前記軌道計算装置は、
前記基準軌道情報を生成する基準軌道情報生成装置から前記基準軌道情報を受信する受信部を備え、
前記格納部は、
前記受信部が受信した前記基準軌道情報を格納する請求項１から請求項３のいずれか一項に記載の軌道計算装置。
【請求項５】
コンピュータに、
人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する処理、
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いることにより、時刻ｔ０における状態量の変化量Δｘ’０を収束判定をすることなく計算し、時刻ｔ０における基準軌道状態量ｘ’０に、計算した前記変化量Δｘ’０を加算することにより、時刻ｔ０における前記人工衛星の状態量の推定値である推定値ｘ０を計算する処理
を実行させるための軌道計算プログラム。
【請求項６】
人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列と、重み行列とを含む基準軌道情報を格納する格納部と、
前記基準軌道情報を用いることにより基準軌道のヤコビアンを求め、前記人工衛星に対して観測された物理量である観測量とヤコビアンと重み行列とを用いて時刻ｔ０における状態量の変化量Δｘ’０を計算し、時刻ｔ０における基準軌道状態量ｘ’０に、計算した前記変化量Δｘ’０を加算することにより、時刻ｔ０における前記人工衛星の状態量の推定値である推定値ｘ０を計算する計算部と
を備える軌道計算装置。
【請求項７】
コンピュータに、
人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列と、重み行列とを含む基準軌道情報を格納する処理、
前記基準軌道情報を用いることにより基準軌道のヤコビアンを求め、前記人工衛星に対して観測された物理量である観測量とヤコビアンと重み行列とを用いて時刻ｔ０における状態量の変化量Δｘ’０を計算し、時刻ｔ０における基準軌道状態量ｘ’０に、計算した前記変化量Δｘ’０を加算することにより、時刻ｔ０における前記人工衛星の状態量の推定値である推定値ｘ０を計算する処理
を実行させるための軌道計算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、衛星が定められた周期で同一の位置を通過する基準軌道の近傍を飛行する場合において、観測により取得されたデータに基づいて衛星の軌道を推定・予測・制御する演算処理を行う装置に関する。
【背景技術】
【０００２】
衛星の利用分野の拡大に伴い、衛星の軌道を推定・予測・制御する性能に対して要求される精度が高くなっている。宇宙から地球の状態をモニタする地球観測衛星においては、観測した地球上の対象位置を特定したり、衛星に作用する加速度を測定したりするために、衛星が飛行した位置を高い精度で推定する必要がある。複数の衛星が定められた間隔で飛行し観測したデータを合成することにより地球上の３次元的な情報を得ることができるため、衛星が定められた位置を飛行するように軌道を予測・制御する技術が求められている。既に市民生活へ浸透している技術であるＧＰＳをはじめとする全地球航法衛星システムにおいては、衛星の電波を受信するユーザが自らの位置を特定するために衛星の軌道情報を利用する必要があり、衛星から配信される軌道情報の精度が測位精度に直接影響する。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特許第３５５０５２０号公報
【特許文献２】特開２０１３−３０９３号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
衛星の軌道を高い精度で推定する方法として、特許文献１および特許文献２で利用されているバッチ最小二乗推定法が用いられることが多い。この方法を用いる場合には、衛星の軌道力学モデルに基づく数値積分等の計算を繰り返し実施する必要があるため計算機の処理負荷が高く、推定に要する処理時間が長くなるという課題があった。
【課題を解決するための手段】
【０００５】
この発明の軌道計算装置は、人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する格納部と、
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いて、前記人工衛星の状態量の推定値を計算する計算部と
を備えることを特徴とする。
【発明の効果】
【０００６】
本発明の軌道計算装置では、基準軌道情報をあらかじめ計算し、格納部にデータとして記録し、記録されたデータを利用して、観測量から状態量の推定値を計算する。よって、処理が簡略化され、計算機である軌道計算装置の処理負荷を低減できる。
【図面の簡単な説明】
【０００７】
【図１】実施の形態１の図で、一般的なバッチ最小二乗推定法の処理のフローチャート。
【図２】実施の形態１の図で、基準軌道を説明する図。
【図３】実施の形態１の図で、状態量の推定値ｘ０の計算方法を模式的に示す図。
【図４】実施の形態１の図で、軌道計算装置１０，基準軌道情報生成装置２０のブロック図。
【図５】実施の形態１の図で、軌道計算装置１０，基準軌道情報生成装置２０のハードウェア構成図。
【図６】実施の形態１の図で、計算部１２の処理を示すフローチャート。
【図７】実施の形態１の図で、基準軌道情報１１Ａの取得経路を示す図。
【図８】実施の形態２の図で、推定値ｘ０から予測値ｘｃを計算する場合を示す図。
【図９】実施の形態３の図で、予測値ｘｃを変更する場合を説明する図。
【発明を実施するための形態】
【０００８】
実施の形態１．
はじめに、一般的なバッチ最小二乗推定法について説明する。
【０００９】
時刻ｔにおける衛星の状態量をｘ（ｔ）と表し、次の式１のように定義する。
【数１】
0006640615.tif 000002
ただし、式１のｐ（ｔ）、ｖ（ｔ）は、それぞれ時刻ｔにおける衛星の位置、速度であり、ｄ（ｔ）は衛星の加速度モデルにおいて推定すべき力学パラメータである。
状態量ｘ（ｔ）は時刻ｔ０における状態量ｘ０＝ｘ（ｔ０）から軌道力学モデルを用いて計算することができる。
状態量ｘ（ｔ）に対する観測モデルをｙ（ｔ）とする。
【００１０】
いま、時刻ｔｋ（ｋ＝１，２，．．．，ｎ）において観測量ｚ（ｔ）が得られたとき、時刻ｔ０における状態量の推定値はｘ０に次の式２の変化量Δｘ０を加算することで得られる。
【数２】
0006640615.tif 000003
ただし、式２において、ｙｙ、ｚｚ、Ｊは次の式３のように表され、Ｗはｚｚの各成分に対する重みを表す行列である。
【数３】
0006640615.tif 000004
【００１１】
式３のＪｋ（ｋ＝１，２，．．．，ｎ）はヤコビアンであり、次の式４により計算される。
【数４】
0006640615.tif 000005
【００１２】
式４のＴｉｊは状態遷移行列と呼ばれ、状態量の時刻ｔｉにおける微小変化δｘｉに対して時刻ｔｊに生じる微小変化δｘｊの関係δｘｊ＝Ｔｉｊδｘｉを表す行列である。
ここで、Ｔ０ｋとＴ０（ｋ＋１）には次の式５の関係がある。
【数５】
0006640615.tif 000006
従って、Ｔ０１とＴｋ（ｋ＋１）（ｋ＝１，２，．．．，ｎ−１）を計算することによりＴ０ｋ（ｋ＝１，２，．．．，ｎ）を得ることができる。
【００１３】
以上のようにして、時刻ｔ０における状態量の推定値が（ｘ０＋Δｘ０）として得られる。その推定値を用いて軌道力学モデルに基づく数値積分等により時刻ｔにおける状態量ｘ（ｔ）を推定することができる。
【００１４】
状態量の推定値を高い精度で得るためには、推定された状態量ｘ（ｔ）を用いて同様の計算を繰り返し、Δｘ０の各成分の絶対値があらかじめ定めた閾値より小さくなると計算が収束したと判定し、推定処理を終了して状態量の推定値を得る。
【００１５】
一般的には、時刻ｔ０における状態量ｘ０から時刻ｔｋ（ｋ＝１，２，．．．，ｎ）における状態量と状態遷移行列を計算するには軌道力学モデルに基づく数値積分等の方法が用いられ、計算機の処理負荷が高くなることがこの方式の課題となっている。
【００１６】
図１は、一般的なバッチ最小二乗推定法を示すフローチャートである。次に、図１を用いて一般的なバッチ最小二乗推定法の処理フローについて説明する。
【００１７】
ステップＳ１１において、推定処理の開始時に、時刻ｔ０における状態量（初期値）を与える。この時点では状態量の変化量は零であり、ステップＳ１２において、状態量（初期値）を状態量（推定値）とする。
【００１８】
ステップＳ１２の状態量（推定値）から、軌道力学モデルに基づいた数値積分等を行い、ステップＳ１３で観測モデルｙｙを計算し、ステップＳ１４でヤコビアンＪを計算する。
【００１９】
ステップＳ１５において、観測により取得された観測量ｚｚからステップＳ１３で計算された観測モデルｙｙを差し引く。なお観測量ｚｚとは人工衛星に対して観測された物理量である。
ステップＳ１６において、それに対してヤコビアンＪと重み行列Ｗとを用いて、状態量ｘ０の変化量Δｘ０を計算する。
【００２０】
ステップＳ１７において、ステップＳ１６で計算された状態量の変化量Δｘ０の各成分ｐ、ｖ、ｄの絶対値と閾値とを比較することにより収束判定を行い、閾値より小さい場合には収束と判定し推定処理を終了する。
【００２１】
ステップＳ１７の収束判定において状態量の変化量Δｘ０の各成分ｐ、ｖ、ｄの絶対値が閾値より大きい場合には、ステップＳ１２の状態量（推定値）に状態量の変化量Δｘ０を加算し、それを改めてステップＳ１２の状態量（推定値）として同様の処理を繰り返す。
【００２２】
以下、実施の形態１に係る発明の具体的な実施の形態を説明する。
【００２３】
図２は、基準軌道を説明する概略図である。人工衛星が周期Ｐで同一の位置を通過する軌道は基準軌道と称される。図２は、人工衛星が周期Ｐで同一の位置を通過する基準軌道４０を示している。人工衛星３０は、周期Ｐで位置５１を通過し、同様に周期Ｐで位置５２を通過し、同様に周期Ｐで位置５３を通過する。同一の位置とは、位置５１、位置５２、位置５３である。
衛星があらかじめ定められた基準軌道の近傍を飛行する場合には、以下のように推定処理を簡略化することができる。図２では人工衛星３０の実際の軌道３１を破線で示している。基準軌道４０は理想化された軌道である。
【００２４】
人工衛星３０が基準軌道４０を飛行する場合には、状態量は周期Ｐで同一の値になる。  このため、
時刻ｔにおける状態量をｘ’（ｔ）、
観測モデルをｙ’（ｔ）と表すと、
ｘ’（ｔ＋Ｐ）＝ｘ’（ｔ）、
ｙ’（ｔ＋Ｐ）＝ｙ’（ｔ）が成り立つ。
なおｘ’は基準軌道状態量である。
従って、
ｔ（ｋ＋ｍ）＝ｔｋ＋Ｐのとき、
基準軌道４０におけるＹｋ、Ｔｋ（ｋ＋１）をそれぞれＹ’ｋ、Ｔ’ｋ（ｋ＋１）と表すと、
Ｙ’（ｋ＋ｍ）＝Ｙ’ｋ、
Ｔ’（ｋ＋ｍ）（ｋ＋ｍ＋１）＝Ｔ’ｋ（ｋ＋１）
が成り立つ。
なおＴ’は基準軌道状態遷移行列である。
【００２５】
人工衛星３０が基準軌道４０の近傍を飛行する場合には、
ｘ（ｔ）≒ｘ’（ｔ）
と近似することができる。
また同様に、
ｙ（ｔｋ）≒ｙ’（ｔｋ）、
Ｔｋ（ｋ＋１）≒Ｔ’ｋ（ｋ＋１）
と近似することができる。
【００２６】
時刻ｔ０における状態量の推定値ｘ０は、基準軌道４０における状態量ｘ’０に次の式６のように、変化量Δｘ’０を加算することにより計算する。
つまり、ｘ０＝ｘ’０＋Δｘ’０      （式６）
で計算することができる。
式６のΔｘ’０は、以下の式７、式８、式９、式１０で計算する。
【数６】
0006640615.tif 000007
【００２７】
ｘ’（ｔｋ）、ｙ’（ｔｋ）、Ｙ’ｋ、Ｔ’（ｋ−１）ｋは、あらかじめ定められた基準軌道４０から計算できる。このため、その全て又は一部をあらかじめ計算して計算機の記憶装置にデータとして予め記録しておくことが可能である。また、これらの値は状態量の推定値に依存しないため、繰り返し計算による収束処理を行う必要がない。
【００２８】
図３は、上記の式６を模式的に示す図である。上記のように、基準軌道４０から計算できるデータを計算機の記憶装置に予め記憶しておき、実施の形態１では、このデータを用いて時刻ｔ０における状態量の推定値ｘ０を計算する。図３において、四角の破線で囲んだデータが、記憶装置に予め記憶されているデータである。これらをまとめて基準軌道情報１１Ａと呼ぶ。
【００２９】
観測量を取得する時刻ｔｓにおける
ｘ’（ｔｓ）、ｙ’（ｔｓ）、Ｙ’ｓ、Ｔ’（ｓ−１）ｓは、
あらかじめ計算された時刻ｔｋ（ｋ＝１，２，．．．，ｎ）に対して周期Ｐの整数（ｎｓ）倍を加算しｔｓ＝ｔｋ＋ｎｓＰが成り立つ時刻ｔｋのデータから得ることができる。
【００３０】
もし、
ｔｓ＝ｔｋ＋ｎｓＰ
が成り立つ時刻ｔｋが存在せず
ｔｓ＝ｔｋ＋Δｔｋ＋ｎｓＰ
となる場合には、｜Δｔｋ｜が最も小さくなる時刻ｔｋを含む１つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。
【００３１】
以上のようにして、記憶装置に記録したデータである基準軌道情報１１Ａを用いて推定処理を行うことにより、計算機の処理負荷を低減することができる。
【００３２】
以下に、記憶装置に基準軌道４０から得られたｘ’（ｔｋ）、ｙ’（ｔｋ）、Ｙ’ｋ、Ｔ’（ｋ−１）ｋ等を含む基準軌道情報１１Ａを予め記憶しておき、基準軌道情報１１Ａのデータ用いて人工衛星３０の軌道の推定値を計算する軌道計算装置１０を説明する。
【００３３】
（＊＊＊構成の説明＊＊＊）
図４は、軌道計算装置１０及び基準軌道情報生成装置２０の機能ブロック図を示す。
図５は、軌道計算装置１０及び基準軌道情報生成装置２０のハードウェア構成を示す。
軌道計算装置１０、基準軌道情報生成装置２０はいずれもコンピュータである。
【００３４】
まず図４の軌道計算装置１０を説明する。軌道計算装置１０は、格納部１１、計算部１２、受信部１３を備えている。格納部１１は、基準軌道情報１１Ａを格納している。基準軌道情報１１Ａとは上記のように、基準軌道４０における人工衛星３０の位置ｐ（ｔ）、速度ｖ（ｔ）、力学パラメータｐ（ｔ）とを表す状態量である基準軌道状態量ｘ’（ｔ）と、基準軌道状態量ｘ’（ｔ）の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列Ｔ’ｉｊとを含む情報である。計算部１２は、基準軌道情報１１Ａと、人工衛星３０に対して観測された物理量である観測量ｚｚとを用いて、人工衛星３０の状態量の推定値ｘ（ｔ０）を計算する。
受信部１３は、後述のように、基準軌道情報生成装置２０の生成した基準軌道情報１１Ａを受信する。格納部１１は受信部１３が受信した基準軌道情報１１Ａを格納する。これは格納部１１に基準軌道情報１１Ａがない場合の他、現在の基準軌道情報１１Ａを新しい基準軌道情報１１Ａに更新する場合である。
【００３５】
次に図４の基準軌道情報生成装置２０を説明する。基準軌道情報生成装置２０は、格納部２１、生成部２２、送信部２３を備えている。格納部２１は、基準軌道情報１１Ａを生成するための元データであるオリジナル情報２１Ａを格納している。生成部２２は、オリジナル情報２１Ａを用いて基準軌道情報１１Ａを生成する。格納部２１が格納している基準軌道情報１１Ａは生成部２２が生成したものである。送信部２３は基準軌道情報１１Ａを軌道計算装置１０へ送信する。
【００３６】
図５に示すように、軌道計算装置１０、基準軌道情報生成装置２０はコンピュータである。
【００３７】
以下では軌道計算装置１０を主体に説明するが、この説明は同じコンピュータである
基準軌道情報生成装置２０にも当てはまる。軌道計算装置１０は、プロセッサ８１、メモリ８２、通信装置８３を備える。プロセッサ８１は、プログラムを実行する。メモリ８２は格納部１１（基準軌道情報生成装置２０では格納部２１）を実現する。メモリ８２には、軌道計算装置１０の場合は基準軌道情報１１Ａが記憶されており、基準軌道情報生成装置２０の場合はオリジナル情報２１Ａ、基準軌道情報１１Ａが記憶されている。またメモリ８２には、図４に示す「計算部１２」（基準軌道情報生成装置２０では「生成部２２」）の機能を実現するプログラムが記憶されている。そして、プロセッサ８１がプログラムを実行して、計算部１２（基準軌道情報生成装置２０では生成部２２）の動作を実行する。通信装置８３は受信部１３（基準軌道情報生成装置２０では送信部２３）を実現する。プロセッサ８１は、プロセッシングを行うＩＣ（Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）である。プロセッサ８１は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、ＤＳＰ（Ｄｉｇｉｔａｌ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｏｒ）等である。図５に示すメモリ８２は、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）等である。メモリ８２には、計算部１２（基準軌道情報生成装置２０では生成部２２）を実現するプログラムの他にＯＳ（Ｏｐｅｒａｔｉｎｇ  Ｓｙｓｔｅｍ）も記憶されている。そして、ＯＳの少なくとも一部がプロセッサ８１により実行される。図５では１つのプロセッサが図示されているが複数のプロセッサを備えていてもよい。計算部１２（基準軌道情報生成装置２０では生成部２２）の処理の結果を示す情報やデータや信号値や変数値が、メモリ８２、又は、プロセッサ８１内のレジスタ又はキャッシュメモリに記憶される。
【００３８】
計算部１２あるいは生成部２２の「部」を、「回路」又は「工程」又は「手順」又は「処理」に読み替えてもよい。また、計算部１２あるいは生成部２２は、ロジックＩＣ（Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）、ＧＡ（Ｇａｔｅ  Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｐｅｃｉｆｉｃ  Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）といった電子回路により実現されてもよい。なお、プロセッサ及び上記の電子回路を総称してプロセッシングサーキットリーともいう。
【００３９】
（＊＊＊動作の説明＊＊＊）
図６は、計算部１２の処理を示すフローチャートである。図６を用いて、軌道計算装置１０の計算部１２が格納部１１に格納されている基準軌道情報１１Ａを用いて実行する、「簡略化されたバッチ最小二乗推定法」の処理フローについて説明する。
【００４０】
ステップＳ２０において、基準軌道情報生成装置２０は、軌道計算装置１０による推定処理の開始前に、基準軌道４０における人工衛星３０の状態量ｘ’、状態遷移行列Ｔ’等を計算し、軌道計算装置１０に送信する。軌道計算装置１０は基準軌道情報１１Ａを格納部１１に記録する。なお、軌道計算装置１０が基準軌道情報生成装置２０を内蔵する構成でも構わない。つまり、軌道計算装置１０が基準軌道情報生成装置２０でもある構成でもよい。
【００４１】
ステップＳ２１において、計算部１２は、推定処理の開始時に、格納部１１に記録された基準軌道４０における状態量ｘ’、状態遷移行列Ｔ’等のデータ読み出しを行い、基準軌道４０の観測モデルｙ’ｙを求め（ステップＳ２２）、また基準軌道４０のヤコビアンＪ’を求める（ステップＳ２３）。
【００４２】
ステップＳ２４において、計算部１２は、観測により取得された観測量ｚｚから基準軌道４０の観測モデルｙ’ｙを差し引く。ステップＳ２５において、計算部１２は、その（ｚｚ−ｙ’ｙ）に対して、基準軌道４０のヤコビアンＪ’と重み行列Ｗ’を用いて、状態量の変化量Δｘ’０の計算を行う。
【００４３】
ステップＳ２６において、計算された状態量の変化量Δｘ’０を、格納部１１に記録されている基準軌道４０における時刻ｔ０の状態量ｘ’０に加算し、それを状態量の推定値ｘ０（図３）として推定処理を終了する。
【００４４】
図３、図６に示すように、計算部１２は、基準軌道情報１１Ａを用いて、時刻ｔ０における、状態量の変化量Δｘ’０を計算し、時刻ｔ０における基準軌道状態量ｘ’０に、計算した変化量Δｘ’０を加算することにより、時刻ｔ０の推定値ｘ０を計算する。
【００４５】
図７は、軌道計算装置１０の基準軌道情報１１Ａの取得経路を示す図である。図７の（ａ）は、軌道計算装置１０が地上側に存在し、軌道計算装置１０が地上側に存在する基準軌道情報生成装置２０から基準軌道情報１１Ａを取得する場合を示している。図７の（ｂ）は、軌道計算装置１０が人工衛星３０に搭載されており、軌道計算装置１０は地上側に存在する基準軌道情報生成装置２０から基準軌道情報１１Ａを取得する場合を示している。なお、軌道計算装置１０が基準軌道情報生成装置２０を内蔵しても良いことは上記で述べたとおりである。
【００４６】
（＊＊＊実施の形態１の効果＊＊＊）
軌道計算装置１０では、基準軌道情報１１Ａをあらかじめ計算し、メモリにデータとして記録し、記録されたデータを利用したバッチ最小二乗推定法により、観測量から状態量の推定処理を行う。よって、処理が簡略化され計算機の処理負荷を低減できる。
【００４７】
また、軌道計算装置１０では、図１のステップＳ１７の判定処理が不要であるので、状態量の推定値を迅速に得ることができる。
【００４８】
実施の形態２．
図８は、実施の形態２を説明するための図である。実施の形態２では、実施の形態１のステップＳ２６で計算された時刻ｔ０の状態量の推定値ｘ０を用いて、計算部１２が、時刻ｔ０よりも未来の時刻である時刻ｔｃにおける状態量の予測値ｘｃを計算する。推定値が観測量ｚｚから計算されるのに対して、「予測値」とは、観測量ｚｚを用いることなく計算される状態量の予測結果である。なお実施の形態２は計算部１２による処理である。
【００４９】
図８に示すように、時刻ｔ０における状態量の推定値ｘ０が得られると、それに基づいて時刻ｔ０よりも未来の時刻ｔｃにおける状態量の予測値ｘｃを、次の式１１に従って計算部１２が計算する。
【数７】
0006640615.tif 000008
【００５０】
時刻ｔｃにおけるｘ’（ｔｃ）、Ｔ’（ｋ−１）ｃは、あらかじめ計算された時刻ｔｋ（ｋ＝１，２，．．．，ｎ）に対して周期Ｐの整数（ｎｃ）倍を加算しｔｃ＝ｔｋ＋ｎｃＰが成り立つ時刻ｔｋのデータから得ることができる。
【００５１】
もし、ｔｃ＝ｔｋ＋ｎｃＰが成り立つ時刻ｔｋが存在せずｔｃ＝ｔｋ＋Δｔｋ＋ｎｐＰとなる場合には、｜Δｔｋ｜が最も小さくなる時刻ｔｋを含む１つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。
【００５２】
以上のように計算部１２は、実施の形態２によれば、格納部１１に記録されたデータと観測量とから推定された衛星の状態量ｘ０を用いて、観測量ｚｚの存在する時刻ｔ０に対して未来の時刻ｔｃにおける人工衛星３０の状態量ｘｃを簡便に予測することができる。
【００５３】
以上のように、計算部１２は、格納部１１に格納された基準軌道情報１１Ａと、計算された状態量の推定値ｘ０とを用いて、新たな観測量を用いることなく計算される人工衛星３０の状態量である予測値ｘｃを計算する。予測値ｘｃは、推定値ｘ０の計算に用いられた観測量ｚｚの存在する時刻よりも未来の時刻ｔｃにおける人工衛星３０の状態量である。
【００５４】
実施の形態３．
図９は、実施の形態３を説明するための図である。実施の形態２では予測値ｘｃを計算した。実施の形態３では、計算部１２は、時刻ｔｃでの予測値ｘｃを予測値ｘｔにしたい場合の、時刻ｔｆ（ｔ０＜ｔｆ＜ｔｃ）における状態量の変化量Δｘｆを計算する場合を説明する。実施の形態３も計算部１２による処理である。
【００５５】
時刻ｔ０における状態量の推定値ｘ０と時刻ｔｃにおける状態量の予測値ｘｃが得られたとき、時刻ｔｃにおける状態量を目標値ｘｔに変更するためには、時刻ｔｆ（ｔ０＜ｔｆ＜ｔｃ）において次の式１２の変化量Δｘｆを時刻ｔｆの状態量ｘｆに加える。変化量Δｘｆは中間変化量である。
【数８】
0006640615.tif 000009
【００５６】
時刻ｔｆにおけるＴ’ｆ（Ｌ＋１）は、あらかじめ計算された時刻ｔＬ（Ｌ＝１，２，．．．，ｎ）に対して周期Ｐの整数（ｎｆ）倍を加算しｔｆ＝ｔＬ＋ｎｆＰが成り立つ時刻ｔＬのデータから得ることができる。
【００５７】
もしｔｆ＝ｔＬ＋ｎｆＰが成り立つ時刻ｔＬが存在せずｔｆ＝ｔＬ＋ΔｔＬ＋ｎｆＰとなる場合には、｜ΔｔＬ｜が最も小さくなる時刻ｔＬを含む１つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。
【００５８】
実際に人工衛星３０が飛行する軌道を制御する際には、時刻ｔｃにおける状態量の目標値として与えられるのは位置ｐｔあるいは位置ｐｔと速度ｖｔであり、一方、時刻ｔｆにおいて状態量に加えることができるのは速度の変化量Δｖｆである。
【００５９】
実施の形態３によれば、格納部１１に記録されたデータと、未来の時刻において予測された衛星の状態量を用いて、その時刻における衛星の状態量を簡便に制御することができる。
【００６０】
なお、図９の場合は時刻ｔｆにおける状態量の変化、つまり、状態量（ｐ、ｖ、ｄ）の変化を計算したが、状態量のうちｐ、ｖを用いて以下のように計算してもよい。
【００６１】
時刻ｔｃにおける状態量の目標値として位置ｐｔが与えられたとき、時刻ｔｆにおいて状態量に加える速度の変化量Δｖｆは次の式１３により計算できる。
【数９】
0006640615.tif 000010
ただし、Ｔ’ｆｃ＿ｐ／ｖはＴ’ｆｃの部分行列であり、次の式１４により表される。
【数１０】
0006640615.tif 000011
【００６２】
時刻ｔｃにおける状態量の目標値として位置ｐｔと速度ｖｔが与えられたとき、２回以上の時刻で状態量に速度の変化量を加える必要がある。時刻ｔｆｊ（ｊ＝１，２，．．．，Ｎ）において状態量に加える速度の変化量Δｖｆｊは次の式１５により計算できる。
【数１１】
0006640615.tif 000012
ただし、行列の右肩の記号＃は一般逆行列を意味する。Ｔ’ｆｊｃ＿ｐｖ／ｖはＴ’ｆｊｃの部分行列であり、次の式１６により表される。
【数１２】
0006640615.tif 000013
【００６３】
時刻ｔｆｊにおけるＴ’ｆｊｃについても、時刻ｔｆにおけるＴ’ｆｃと同様の方法で計算することができる。
【００６４】
以上のように、計算部１２は、基準軌道情報１１Ａと、状態量の予測値ｘｃとを用いて、中間変化量Δｘｆを計算する。中間変化量Δｘｆは、観測量ｚｚの存在する時刻ｔ０よりも未来の時刻ｔｃにおける人工衛星３０の予測値ｘｃを異なる予測値ｘｔに変更するための状態量の変化量であり、未来の時刻ｔｃと観測量ｚｚの存在する時刻ｔ０との間の時刻ｔｆにおける状態量の変化量である。
【符号の説明】
【００６５】
１０  軌道計算装置、１１  格納部、１２  計算部、１３  受信部、２０  基準軌道情報生成装置、２１  格納部、２２  生成部、２３  送信部、３０  人工衛星、４０  基準軌道、８１  プロセッサ、８２  メモリ、８３  通信装置。
【図１】
0006640615.tif 000014
【図２】
0006640615.tif 000015
【図３】
0006640615.tif 000016
【図４】
0006640615.tif 000017
【図５】
0006640615.tif 000018
【図６】
0006640615.tif 000019
【図７】
0006640615.tif 000020
【図８】
0006640615.tif 000021
【図９】
0006640615.tif 000022