JP2017161420A - 軌道計算装置及び軌道計算プログラム - Google Patents

Abstract

【課題】人工衛星の軌道の推定において、計算機の処理負荷を抑え、推定に要する処理時間の短い計算装置を提供する。
【解決手段】軌道計算装置１０は、基準軌道情報１１Ａを格納する格納部１１と、基準軌道情報１１Ａを用いて人工衛星の軌道の推定値を計算する計算部１２を備える。基準軌道情報１１Ａは、人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む情報である。計算部１２は、基準軌道情報１１Ａと、人工衛星の観測量とを用いて、状態量の推定値を計算する。
【選択図】図４

Description

本発明は、衛星が定められた周期で同一の位置を通過する基準軌道の近傍を飛行する場合において、観測により取得されたデータに基づいて衛星の軌道を推定・予測・制御する演算処理を行う装置に関する。

衛星の利用分野の拡大に伴い、衛星の軌道を推定・予測・制御する性能に対して要求される精度が高くなっている。宇宙から地球の状態をモニタする地球観測衛星においては、観測した地球上の対象位置を特定したり、衛星に作用する加速度を測定したりするために、衛星が飛行した位置を高い精度で推定する必要がある。複数の衛星が定められた間隔で飛行し観測したデータを合成することにより地球上の３次元的な情報を得ることができるため、衛星が定められた位置を飛行するように軌道を予測・制御する技術が求められている。既に市民生活へ浸透している技術であるＧＰＳをはじめとする全地球航法衛星システムにおいては、衛星の電波を受信するユーザが自らの位置を特定するために衛星の軌道情報を利用する必要があり、衛星から配信される軌道情報の精度が測位精度に直接影響する。

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

衛星の軌道を高い精度で推定する方法として、特許文献１および特許文献２で利用されているバッチ最小二乗推定法が用いられることが多い。この方法を用いる場合には、衛星の軌道力学モデルに基づく数値積分等の計算を繰り返し実施する必要があるため計算機の処理負荷が高く、推定に要する処理時間が長くなるという課題があった。

この発明の軌道計算装置は、人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する格納部と、
前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いて、前記人工衛星の状態量の推定値を計算する計算部と
を備えることを特徴とする。

本発明の軌道計算装置では、基準軌道情報をあらかじめ計算し、格納部にデータとして記録し、記録されたデータを利用して、観測量から状態量の推定値を計算する。よって、処理が簡略化され、計算機である軌道計算装置の処理負荷を低減できる。

実施の形態１の図で、一般的なバッチ最小二乗推定法の処理のフローチャート。 実施の形態１の図で、基準軌道を説明する図。 実施の形態１の図で、状態量の推定値ｘ_０の計算方法を模式的に示す図。 実施の形態１の図で、軌道計算装置１０，基準軌道情報生成装置２０のブロック図。 実施の形態１の図で、軌道計算装置１０，基準軌道情報生成装置２０のハードウェア構成図。 実施の形態１の図で、計算部１２の処理を示すフローチャート。 実施の形態１の図で、基準軌道情報１１Ａの取得経路を示す図。 実施の形態２の図で、推定値ｘ_０から予測値ｘ_ｃを計算する場合を示す図。 実施の形態３の図で、予測値ｘ_ｃを変更する場合を説明する図。

実施の形態１．
はじめに、一般的なバッチ最小二乗推定法について説明する。

時刻ｔにおける衛星の状態量をｘ（ｔ）と表し、次の式１のように定義する。

ただし、式１のｐ（ｔ）、ｖ（ｔ）は、それぞれ時刻ｔにおける衛星の位置、速度であり、ｄ（ｔ）は衛星の加速度モデルにおいて推定すべき力学パラメータである。
状態量ｘ（ｔ）は時刻ｔ_０における状態量ｘ_０＝ｘ（ｔ_０）から軌道力学モデルを用いて計算することができる。
状態量ｘ（ｔ）に対する観測モデルをｙ（ｔ）とする。

いま、時刻ｔ_ｋ（ｋ＝１，２，．．．，ｎ）において観測量ｚ（ｔ）が得られたとき、時刻ｔ_０における状態量の推定値はｘ_０に次の式２の変化量Δｘ_０を加算することで得られる。

ただし、式２において、ｙ_ｙ、ｚ_ｚ、Ｊは次の式３のように表され、Ｗはｚ_ｚの各成分に対する重みを表す行列である。

式３のＪ_ｋ（ｋ＝１，２，．．．，ｎ）はヤコビアンであり、次の式４により計算される。

式４のＴ_ｉｊは状態遷移行列と呼ばれ、状態量の時刻ｔ_ｉにおける微小変化δｘ_ｉに対して時刻ｔ_ｊに生じる微小変化δｘ_ｊの関係δｘ_ｊ＝Ｔ_ｉｊδｘ_ｉを表す行列である。
ここで、Ｔ_０ｋとＴ_{０（ｋ＋１）}には次の式５の関係がある。

従って、Ｔ_０１とＴ_{ｋ（ｋ＋１）}（ｋ＝１，２，．．．，ｎ−１）を計算することによりＴ_０ｋ（ｋ＝１，２，．．．，ｎ）を得ることができる。

以上のようにして、時刻ｔ_０における状態量の推定値が（ｘ_０＋Δｘ_０）として得られる。その推定値を用いて軌道力学モデルに基づく数値積分等により時刻ｔにおける状態量ｘ（ｔ）を推定することができる。

状態量の推定値を高い精度で得るためには、推定された状態量ｘ（ｔ）を用いて同様の計算を繰り返し、Δｘ_０の各成分の絶対値があらかじめ定めた閾値より小さくなると計算が収束したと判定し、推定処理を終了して状態量の推定値を得る。

一般的には、時刻ｔ_０における状態量ｘ_０から時刻ｔ_ｋ（ｋ＝１，２，．．．，ｎ）における状態量と状態遷移行列を計算するには軌道力学モデルに基づく数値積分等の方法が用いられ、計算機の処理負荷が高くなることがこの方式の課題となっている。

図１は、一般的なバッチ最小二乗推定法を示すフローチャートである。次に、図１を用いて一般的なバッチ最小二乗推定法の処理フローについて説明する。

ステップＳ１１において、推定処理の開始時に、時刻ｔ_０における状態量（初期値）を与える。この時点では状態量の変化量は零であり、ステップＳ１２において、状態量（初期値）を状態量（推定値）とする。

ステップＳ１２の状態量（推定値）から、軌道力学モデルに基づいた数値積分等を行い、ステップＳ１３で観測モデルｙ_ｙを計算し、ステップＳ１４でヤコビアンＪを計算する。

ステップＳ１５において、観測により取得された観測量ｚ_ｚからステップＳ１３で計算された観測モデルｙ_ｙを差し引く。なお観測量ｚ_ｚとは人工衛星に対して観測された物理量である。
ステップＳ１６において、それに対してヤコビアンＪと重み行列Ｗとを用いて、状態量ｘ_０の変化量Δｘ_０を計算する。

ステップＳ１７において、ステップＳ１６で計算された状態量の変化量Δｘ_０の各成分ｐ、ｖ、ｄの絶対値と閾値とを比較することにより収束判定を行い、閾値より小さい場合には収束と判定し推定処理を終了する。

ステップＳ１７の収束判定において状態量の変化量Δｘ_０の各成分ｐ、ｖ、ｄの絶対値が閾値より大きい場合には、ステップＳ１２の状態量（推定値）に状態量の変化量Δｘ_０を加算し、それを改めてステップＳ１２の状態量（推定値）として同様の処理を繰り返す。

以下、実施の形態１に係る発明の具体的な実施の形態を説明する。

図２は、基準軌道を説明する概略図である。人工衛星が周期Ｐで同一の位置を通過する軌道は基準軌道と称される。図２は、人工衛星が周期Ｐで同一の位置を通過する基準軌道４０を示している。人工衛星３０は、周期Ｐで位置５１を通過し、同様に周期Ｐで位置５２を通過し、同様に周期Ｐで位置５３を通過する。同一の位置とは、位置５１、位置５２、位置５３である。
衛星があらかじめ定められた基準軌道の近傍を飛行する場合には、以下のように推定処理を簡略化することができる。図２では人工衛星３０の実際の軌道３１を破線で示している。基準軌道４０は理想化された軌道である。

人工衛星３０が基準軌道４０を飛行する場合には、状態量は周期Ｐで同一の値になる。 このため、
時刻ｔにおける状態量をｘ’（ｔ）、
観測モデルをｙ’（ｔ）と表すと、
ｘ’（ｔ＋Ｐ）＝ｘ’（ｔ）、
ｙ’（ｔ＋Ｐ）＝ｙ’（ｔ）が成り立つ。
なおｘ’は基準軌道状態量である。
従って、
ｔ_{（ｋ＋ｍ）}＝ｔ_ｋ＋Ｐのとき、
基準軌道４０におけるＹ_ｋ、Ｔ_{ｋ（ｋ＋１）}をそれぞれＹ’_ｋ、Ｔ’_{ｋ（ｋ＋１）}と表すと、
Ｙ’_{（ｋ＋ｍ）}＝Ｙ’_ｋ、
Ｔ’_{（ｋ＋ｍ）（ｋ＋ｍ＋１）}＝Ｔ’_{ｋ（ｋ＋１）}
が成り立つ。
なおＴ’は基準軌道状態遷移行列である。

人工衛星３０が基準軌道４０の近傍を飛行する場合には、
ｘ（ｔ）≒ｘ’（ｔ）
と近似することができる。
また同様に、
ｙ（ｔ_ｋ）≒ｙ’（ｔ_ｋ）、
Ｔ_{ｋ（ｋ＋１）}≒Ｔ’_{ｋ（ｋ＋１）}
と近似することができる。

時刻ｔ_０における状態量の推定値ｘ_０は、基準軌道４０における状態量ｘ’_０に次の式６のように、変化量Δｘ’_０を加算することにより計算する。
つまり、ｘ_０＝ｘ’_０＋Δｘ’_０ （式６）
で計算することができる。
式６のΔｘ’_０は、以下の式７、式８、式９、式１０で計算する。

ｘ’（ｔ_ｋ）、ｙ’（ｔ_ｋ）、Ｙ’_ｋ、Ｔ’_{（ｋ−１）ｋ}は、あらかじめ定められた基準軌道４０から計算できる。このため、その全て又は一部をあらかじめ計算して計算機の記憶装置にデータとして予め記録しておくことが可能である。また、これらの値は状態量の推定値に依存しないため、繰り返し計算による収束処理を行う必要がない。

図３は、上記の式６を模式的に示す図である。上記のように、基準軌道４０から計算できるデータを計算機の記憶装置に予め記憶しておき、実施の形態１では、このデータを用いて時刻ｔ_０における状態量の推定値ｘ_０を計算する。図３において、四角の破線で囲んだデータが、記憶装置に予め記憶されているデータである。これらをまとめて基準軌道情報１１Ａと呼ぶ。

観測量を取得する時刻ｔ_ｓにおける
ｘ’（ｔ_ｓ）、ｙ’（ｔ_ｓ）、Ｙ’_ｓ、Ｔ’_{（ｓ−１）ｓ}は、
あらかじめ計算された時刻ｔ_ｋ（ｋ＝１，２，．．．，ｎ）に対して周期Ｐの整数（ｎ_ｓ）倍を加算しｔ_ｓ＝ｔ_ｋ＋ｎ_ｓＰが成り立つ時刻ｔ_ｋのデータから得ることができる。

もし、
ｔ_ｓ＝ｔ_ｋ＋ｎ_ｓＰ
が成り立つ時刻ｔｋが存在せず
ｔ_ｓ＝ｔ_ｋ＋Δｔ_ｋ＋ｎ_ｓＰ
となる場合には、｜Δｔ_ｋ｜が最も小さくなる時刻ｔ_ｋを含む１つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。

以上のようにして、記憶装置に記録したデータである基準軌道情報１１Ａを用いて推定処理を行うことにより、計算機の処理負荷を低減することができる。

以下に、記憶装置に基準軌道４０から得られたｘ’（ｔ_ｋ）、ｙ’（ｔ_ｋ）、Ｙ’_ｋ、Ｔ’_{（ｋ−１）ｋ}等を含む基準軌道情報１１Ａを予め記憶しておき、基準軌道情報１１Ａのデータ用いて人工衛星３０の軌道の推定値を計算する軌道計算装置１０を説明する。

（＊＊＊構成の説明＊＊＊）
図４は、軌道計算装置１０及び基準軌道情報生成装置２０の機能ブロック図を示す。
図５は、軌道計算装置１０及び基準軌道情報生成装置２０のハードウェア構成を示す。
軌道計算装置１０、基準軌道情報生成装置２０はいずれもコンピュータである。

まず図４の軌道計算装置１０を説明する。軌道計算装置１０は、格納部１１、計算部１２、受信部１３を備えている。格納部１１は、基準軌道情報１１Ａを格納している。基準軌道情報１１Ａとは上記のように、基準軌道４０における人工衛星３０の位置ｐ（ｔ）、速度ｖ（ｔ）、力学パラメータｐ（ｔ）とを表す状態量である基準軌道状態量ｘ’（ｔ）と、基準軌道状態量ｘ’（ｔ）の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列Ｔ’_ｉｊとを含む情報である。計算部１２は、基準軌道情報１１Ａと、人工衛星３０に対して観測された物理量である観測量ｚ_ｚとを用いて、人工衛星３０の状態量の推定値ｘ（ｔ_０）を計算する。
受信部１３は、後述のように、基準軌道情報生成装置２０の生成した基準軌道情報１１Ａを受信する。格納部１１は受信部１３が受信した基準軌道情報１１Ａを格納する。これは格納部１１に基準軌道情報１１Ａがない場合の他、現在の基準軌道情報１１Ａを新しい基準軌道情報１１Ａに更新する場合である。

次に図４の基準軌道情報生成装置２０を説明する。基準軌道情報生成装置２０は、格納部２１、生成部２２、送信部２３を備えている。格納部２１は、基準軌道情報１１Ａを生成するための元データであるオリジナル情報２１Ａを格納している。生成部２２は、オリジナル情報２１Ａを用いて基準軌道情報１１Ａを生成する。格納部２１が格納している基準軌道情報１１Ａは生成部２２が生成したものである。送信部２３は基準軌道情報１１Ａを軌道計算装置１０へ送信する。

図５に示すように、軌道計算装置１０、基準軌道情報生成装置２０はコンピュータである。

以下では軌道計算装置１０を主体に説明するが、この説明は同じコンピュータである
基準軌道情報生成装置２０にも当てはまる。軌道計算装置１０は、プロセッサ８１、メモリ８２、通信装置８３を備える。プロセッサ８１は、プログラムを実行する。メモリ８２は格納部１１（基準軌道情報生成装置２０では格納部２１）を実現する。メモリ８２には、軌道計算装置１０の場合は基準軌道情報１１Ａが記憶されており、基準軌道情報生成装置２０の場合はオリジナル情報２１Ａ、基準軌道情報１１Ａが記憶されている。またメモリ８２には、図４に示す「計算部１２」（基準軌道情報生成装置２０では「生成部２２」）の機能を実現するプログラムが記憶されている。そして、プロセッサ８１がプログラムを実行して、計算部１２（基準軌道情報生成装置２０では生成部２２）の動作を実行する。通信装置８３は受信部１３（基準軌道情報生成装置２０では送信部２３）を実現する。プロセッサ８１は、プロセッシングを行うＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）である。プロセッサ８１は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）等である。図５に示すメモリ８２は、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）等である。メモリ８２には、計算部１２（基準軌道情報生成装置２０では生成部２２）を実現するプログラムの他にＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）も記憶されている。そして、ＯＳの少なくとも一部がプロセッサ８１により実行される。図５では１つのプロセッサが図示されているが複数のプロセッサを備えていてもよい。計算部１２（基準軌道情報生成装置２０では生成部２２）の処理の結果を示す情報やデータや信号値や変数値が、メモリ８２、又は、プロセッサ８１内のレジスタ又はキャッシュメモリに記憶される。

計算部１２あるいは生成部２２の「部」を、「回路」又は「工程」又は「手順」又は「処理」に読み替えてもよい。また、計算部１２あるいは生成部２２は、ロジックＩＣ（Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＧＡ（Ｇａｔｅ Ａｒｒａｙ）、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）といった電子回路により実現されてもよい。なお、プロセッサ及び上記の電子回路を総称してプロセッシングサーキットリーともいう。

（＊＊＊動作の説明＊＊＊）
図６は、計算部１２の処理を示すフローチャートである。図６を用いて、軌道計算装置１０の計算部１２が格納部１１に格納されている基準軌道情報１１Ａを用いて実行する、「簡略化されたバッチ最小二乗推定法」の処理フローについて説明する。

ステップＳ２０において、基準軌道情報生成装置２０は、軌道計算装置１０による推定処理の開始前に、基準軌道４０における人工衛星３０の状態量ｘ’、状態遷移行列Ｔ’等を計算し、軌道計算装置１０に送信する。軌道計算装置１０は基準軌道情報１１Ａを格納部１１に記録する。なお、軌道計算装置１０が基準軌道情報生成装置２０を内蔵する構成でも構わない。つまり、軌道計算装置１０が基準軌道情報生成装置２０でもある構成でもよい。

ステップＳ２１において、計算部１２は、推定処理の開始時に、格納部１１に記録された基準軌道４０における状態量ｘ’、状態遷移行列Ｔ’等のデータ読み出しを行い、基準軌道４０の観測モデルｙ’_ｙを求め（ステップＳ２２）、また基準軌道４０のヤコビアンＪ’を求める（ステップＳ２３）。

ステップＳ２４において、計算部１２は、観測により取得された観測量ｚ_ｚから基準軌道４０の観測モデルｙ’_ｙを差し引く。ステップＳ２５において、計算部１２は、その（ｚ_ｚ−ｙ’_ｙ）に対して、基準軌道４０のヤコビアンＪ’と重み行列Ｗ’を用いて、状態量の変化量Δｘ’_０の計算を行う。

ステップＳ２６において、計算された状態量の変化量Δｘ’_０を、格納部１１に記録されている基準軌道４０における時刻ｔ_０の状態量ｘ’_０に加算し、それを状態量の推定値ｘ_０（図３）として推定処理を終了する。

図３、図６に示すように、計算部１２は、基準軌道情報１１Ａを用いて、時刻ｔ_０における、状態量の変化量Δｘ’_０を計算し、時刻ｔ_０における基準軌道状態量ｘ’_０に、計算した変化量Δｘ’_０を加算することにより、時刻ｔ_０の推定値ｘ_０を計算する。

図７は、軌道計算装置１０の基準軌道情報１１Ａの取得経路を示す図である。図７の（ａ）は、軌道計算装置１０が地上側に存在し、軌道計算装置１０が地上側に存在する基準軌道情報生成装置２０から基準軌道情報１１Ａを取得する場合を示している。図７の（ｂ）は、軌道計算装置１０が人工衛星３０に搭載されており、軌道計算装置１０は地上側に存在する基準軌道情報生成装置２０から基準軌道情報１１Ａを取得する場合を示している。なお、軌道計算装置１０が基準軌道情報生成装置２０を内蔵しても良いことは上記で述べたとおりである。

（＊＊＊実施の形態１の効果＊＊＊）
軌道計算装置１０では、基準軌道情報１１Ａをあらかじめ計算し、メモリにデータとして記録し、記録されたデータを利用したバッチ最小二乗推定法により、観測量から状態量の推定処理を行う。よって、処理が簡略化され計算機の処理負荷を低減できる。

また、軌道計算装置１０では、図１のステップＳ１７の判定処理が不要であるので、状態量の推定値を迅速に得ることができる。

実施の形態２．
図８は、実施の形態２を説明するための図である。実施の形態２では、実施の形態１のステップＳ２６で計算された時刻ｔ_０の状態量の推定値ｘ_０を用いて、計算部１２が、時刻ｔ_０よりも未来の時刻である時刻ｔ_ｃにおける状態量の予測値ｘ_ｃを計算する。推定値が観測量ｚ_ｚから計算されるのに対して、「予測値」とは、観測量ｚ_ｚを用いることなく計算される状態量の予測結果である。なお実施の形態２は計算部１２による処理である。

図８に示すように、時刻ｔ_０における状態量の推定値ｘ_０が得られると、それに基づいて時刻ｔ_０よりも未来の時刻ｔ_ｃにおける状態量の予測値ｘ_ｃを、次の式１１に従って計算部１２が計算する。

時刻ｔ_ｃにおけるｘ’（ｔ_ｃ）、Ｔ’_{（ｋ−１）ｃ}は、あらかじめ計算された時刻ｔ_ｋ（ｋ＝１，２，．．．，ｎ）に対して周期Ｐの整数（ｎ_ｃ）倍を加算しｔ_ｃ＝ｔ_ｋ＋ｎ_ｃＰが成り立つ時刻ｔ_ｋのデータから得ることができる。

もし、ｔ_ｃ＝ｔ_ｋ＋ｎ_ｃＰが成り立つ時刻ｔ_ｋが存在せずｔ_ｃ＝ｔ_ｋ＋Δｔ_ｋ＋ｎ_ｐＰとなる場合には、｜Δｔ_ｋ｜が最も小さくなる時刻ｔ_ｋを含む１つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。

以上のように計算部１２は、実施の形態２によれば、格納部１１に記録されたデータと観測量とから推定された衛星の状態量ｘ_０を用いて、観測量ｚ_ｚの存在する時刻ｔ_０に対して未来の時刻ｔ_ｃにおける人工衛星３０の状態量ｘ_ｃを簡便に予測することができる。

以上のように、計算部１２は、格納部１１に格納された基準軌道情報１１Ａと、計算された状態量の推定値ｘ_０とを用いて、新たな観測量を用いることなく計算される人工衛星３０の状態量である予測値ｘ_ｃを計算する。予測値ｘ_ｃは、推定値ｘ_０の計算に用いられた観測量ｚ_ｚの存在する時刻よりも未来の時刻ｔ_ｃにおける人工衛星３０の状態量である。

実施の形態３．
図９は、実施の形態３を説明するための図である。実施の形態２では予測値ｘ_ｃを計算した。実施の形態３では、計算部１２は、時刻ｔ_ｃでの予測値ｘ_ｃを予測値ｘ_ｔにしたい場合の、時刻ｔ_ｆ（ｔ_０＜ｔ_ｆ＜ｔ_ｃ）における状態量の変化量Δｘ_ｆを計算する場合を説明する。実施の形態３も計算部１２による処理である。

時刻ｔ_０における状態量の推定値ｘ_０と時刻ｔ_ｃにおける状態量の予測値ｘ_ｃが得られたとき、時刻ｔ_ｃにおける状態量を目標値ｘ_ｔに変更するためには、時刻ｔ_ｆ（ｔ_０＜ｔ_ｆ＜ｔ_ｃ）において次の式１２の変化量Δｘ_ｆを時刻ｔ_ｆの状態量ｘ_ｆに加える。変化量Δｘ_ｆは中間変化量である。

時刻ｔ_ｆにおけるＴ’_{ｆ（Ｌ＋１）}は、あらかじめ計算された時刻ｔ_Ｌ（Ｌ＝１，２，．．．，ｎ）に対して周期Ｐの整数（ｎ_ｆ）倍を加算しｔ_ｆ＝ｔ_Ｌ＋ｎ_ｆＰが成り立つ時刻ｔ_Ｌのデータから得ることができる。

もしｔ_ｆ＝ｔ_Ｌ＋ｎ_ｆＰが成り立つ時刻ｔ_Ｌが存在せずｔ_ｆ＝ｔ_Ｌ＋Δｔ_Ｌ＋ｎ_ｆＰとなる場合には、｜Δｔ_Ｌ｜が最も小さくなる時刻ｔ_Ｌを含む１つまたは複数の時刻のデータを用いて、外挿あるいは内挿等の方法により計算することができる。

実際に人工衛星３０が飛行する軌道を制御する際には、時刻ｔ_ｃにおける状態量の目標値として与えられるのは位置ｐ_ｔあるいは位置ｐ_ｔと速度ｖ_ｔであり、一方、時刻ｔ_ｆにおいて状態量に加えることができるのは速度の変化量Δｖ_ｆである。

実施の形態３によれば、格納部１１に記録されたデータと、未来の時刻において予測された衛星の状態量を用いて、その時刻における衛星の状態量を簡便に制御することができる。

なお、図９の場合は時刻ｔ_ｆにおける状態量の変化、つまり、状態量（ｐ、ｖ、ｄ）の変化を計算したが、状態量のうちｐ、ｖを用いて以下のように計算してもよい。

時刻ｔ_ｃにおける状態量の目標値として位置ｐ_ｔが与えられたとき、時刻ｔ_ｆにおいて状態量に加える速度の変化量Δｖ_ｆは次の式１３により計算できる。

ただし、Ｔ’_{ｆｃ＿ｐ／ｖ}はＴ’_ｆｃの部分行列であり、次の式１４により表される。

時刻ｔ_ｃにおける状態量の目標値として位置ｐ_ｔと速度ｖ_ｔが与えられたとき、２回以上の時刻で状態量に速度の変化量を加える必要がある。時刻ｔ_ｆｊ（ｊ＝１，２，．．．，Ｎ）において状態量に加える速度の変化量Δｖ_ｆｊは次の式１５により計算できる。

ただし、行列の右肩の記号＃は一般逆行列を意味する。Ｔ’_{ｆｊｃ＿ｐｖ／ｖ}はＴ’_ｆｊｃの部分行列であり、次の式１６により表される。

時刻ｔ_ｆｊにおけるＴ’_ｆｊｃについても、時刻ｔ_ｆにおけるＴ’_ｆｃと同様の方法で計算することができる。

以上のように、計算部１２は、基準軌道情報１１Ａと、状態量の予測値ｘ_ｃとを用いて、中間変化量Δｘ_ｆを計算する。中間変化量Δｘ_ｆは、観測量ｚ_ｚの存在する時刻ｔ_０よりも未来の時刻ｔ_ｃにおける人工衛星３０の予測値ｘ_ｃを異なる予測値ｘ_ｔに変更するための状態量の変化量であり、未来の時刻ｔ_ｃと観測量ｚ_ｚの存在する時刻ｔ_０との間の時刻ｔ_ｆにおける状態量の変化量である。

１０ 軌道計算装置、１１ 格納部、１２ 計算部、１３ 受信部、２０ 基準軌道情報生成装置、２１ 格納部、２２ 生成部、２３ 送信部、３０ 人工衛星、４０ 基準軌道、８１ プロセッサ、８２ メモリ、８３ 通信装置。

Claims (6)

1. 人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する格納部と、
   前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いて、前記人工衛星の状態量の推定値を計算する計算部と
   を備える軌道計算装置。
2. 前記計算部は、
   前記基準軌道情報を用いて、時刻ｔ_０における、状態量の変化量Δｘ’_０を計算し、時刻ｔ_０における基準軌道状態量ｘ’_０に、計算した変化量Δｘ’_０を加算することにより、時刻ｔ_０の前記推定値である推定値ｘ_０を計算する請求項１に記載の軌道計算装置。
3. 前記計算部は、
   前記格納部に格納された前記基準軌道情報と、計算された前記状態量の前記推定値とを用いて、前記推定値の計算に用いられた前記観測量の存在する時刻よりも未来の時刻における前記人工衛星の状態量であり、新たな観測量を用いることなく計算される前記人工衛星の状態量である予測値を計算する請求項１または請求項２に記載の軌道計算装置。
4. 前記計算部は、
   前記基準軌道情報と、前記状態量の前記予測値とを用いて、前記観測量の存在する前記時刻よりも前記未来の時刻における前記人工衛星の前記予測値を異なる予測値に変更するための状態量の変化量であって、前記未来の時刻と前記観測量の存在する前記時刻との間の時刻における状態量の変化量である中間変化量を計算する請求項３に記載の軌道計算装置。
5. 前記軌道計算装置は、
   前記基準軌道情報を生成する基準軌道情報生成装置から前記基準軌道情報を受信する受信部を備え、
   前記格納部は、
   前記受信部が受信した前記基準軌道情報を格納する請求項１から請求項４のいずれか一項に記載の軌道計算装置。
6. コンピュータに、
   人工衛星が定められた周期で同一の位置を通過する軌道である基準軌道における前記人工衛星の位置と、速度と、力学パラメータとを表す状態量である基準軌道状態量と、前記基準軌道状態量の微小変化の時間的な遷移を表す状態遷移行列である基準軌道状態遷移行列とを含む基準軌道情報を格納する処理、
   前記基準軌道情報と、前記人工衛星に対して観測された物理量である観測量とを用いて、前記人工衛星の状態量の推定値を計算する処理
   を実行させるための軌道計算プログラム。

Images