(57)【要約】【課題】摂動や外乱の影響を考慮してできるだけ効率的な準天頂衛星の軌道面内保持制御計画を立案する。【解決手段】本方法は、中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出し、摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られるeベクトルの変化に基づき、eベクトルの制御量を算出し、軌道長半径目標制御量及びeベクトルの制御量と、制御実施可能期間における3つの制御実施緯度と、各制御実施緯度における軌道接線方向の制御量との関係を表す3元連立方程式を解いて軌道接線方向の制御量の初期値を算出し、外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに軌道長半径の差分及びeベクトルの差分が所定レベル以下になるように、軌道接線方向の制御量の初期値を補正する。【選択図】図７
2009227019.tif 000002
【特許請求の範囲】
【請求項１】
準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成プログラムであって、
コンピュータに、
中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出するステップと、
今回制御実施時における離心率ベクトルを摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、前記離心率ベクトルの制御量を算出するステップと、
前記軌道長半径目標制御量及び前記離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各前記制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、前記軌道接線方向の制御量の初期値を算出するステップと、
外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに前記軌道長半径の差分及び前記離心率ベクトルの差分を所定レベル以下になるように、前記軌道接線方向の制御量の初期値を補正する補正ステップと、
を実行させるための準天頂衛星の制御計画生成プログラム。
【請求項２】
前記補正ステップが、
アンローディング外乱の影響を考慮して次回制御実施時まで軌道伝播することにより実現中心経度を算出し、当該実現中心経度と前記目標中心経度との差が第１の所定レベル以下となるように、前記軌道接線方向の制御量の初期値を補正するステップと、
前記準天頂衛星のスラスタ取り付け方向のずれの影響を考慮して次回制御時まで軌道伝播することにより実現軌道長半径及び実現離心率ベクトルを算出し、前記実現軌道長半径から得られる前記軌道長半径についての差分及び前記実現離心率ベクトルから得られる前記離心率ベクトルの差分が第２の所定レベル以下となるように、補正された前記軌道接線方向の制御量の初期値をさらに補正するステップと、
を含む請求項１記載の準天頂衛星の制御計画生成プログラム。
【請求項３】
準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成方法であって、
中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出するステップと、
今回制御実施時における離心率ベクトルを摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、前記離心率ベクトルの制御量を算出するステップと、
前記軌道長半径目標制御量及び前記離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各前記制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、前記軌道接線方向の制御量の初期値を算出するステップと、
外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに前記軌道長半径の差分及び前記離心率ベクトルの差分を所定レベル以下になるように、前記軌道接線方向の制御量の初期値を補正する補正ステップと、
を含み、コンピュータにより実行される準天頂衛星の制御計画生成方法。
【請求項４】
準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成装置であって、
中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出する手段と、
今回制御実施時における離心率ベクトルを摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、前記離心率ベクトルの制御量を算出する手段と、
前記軌道長半径目標制御量及び前記離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各前記制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、前記軌道接線方向の制御量の初期値を算出する手段と、
外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに前記軌道長半径の差分及び前記離心率ベクトルの差分を所定レベル以下になるように、前記軌道接線方向の制御量の初期値を補正する補正手段と、
を有する準天頂衛星の制御計画生成装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、人工衛星の軌道制御技術に関し、より詳しくは人工衛星の追跡管制における軌道保持制御技術に関する。
【背景技術】
【０００２】
まず図１Ａ乃至図５を用いて、人工衛星の軌道制御の基本的事項について説明する。図１Ａには、地球１０００と、人工衛星の軌道１００２とが示されている。ここで、地球１０００の中心を通り且つ軌道１００２上の２点を結ぶ線分であってその長さが最大となる線分の端点のうち地球１０００に近い方の点をペリジ点と呼び、遠い方の点をアポジ点と呼ぶ。地球１０００の中心とペリジ点との距離はペリジ半径と呼ばれ、地球１０００の中心とアポジ点との距離はアポジ半径と呼ばれる。また、軌道長半径ａは上記線分の長さの半分である。すなわち、上記線分の長さは軌道長半径ａの２倍となっている。さらに、離心率をｅとすると、アポジ半径はａ×（１＋ｅ）で算出され、ペリジ半径はａ×（１−ｅ）で算出される。なお、地球１０００の中心からペリジ点に向かう長さｅのベクトルをｅベクトル（離心率ベクトル）と呼ぶ。
【０００３】
また、図１Ｂに３次元的に示すように、人工衛星が地球１０００の赤道面を下から上に突き抜ける点を昇交点と呼び、軌道１００２の面と赤道面とが成す角度を軌道傾斜角ｉと呼ぶ。また、地球１０００の春分点を基準として、地球１０００の中心と昇交点とを結ぶ直線の角度を昇交点赤経Ωと呼ぶ。また、地球１０００の中心と昇交点とを結ぶ直線とｅベクトルとが成す角度を近地点引数ωと呼ぶ。さらにｅベクトルと、人工衛星１００４と地球１０００の中心を結ぶ直線とが成す角度を真近点離角ｆと呼ぶ。
【０００４】
人工衛星の軌道には、地球の非球状効果、太陽や月の引力など、さまざまな外乱（摂動）が作用するので、目的とする軌道状態を保持するためには摂動を打ち消す軌道保持制御を定期的に実施する必要がある。軌道保持制御は、軌道面、すなわち軌道傾斜角ｉと昇交点赤経Ω を補正する軌道面制御と、軌道の形、すなわち軌道長半径ａと離心率ｅと近地点引数ωを補正する軌道面内制御に分けられる。通常、離心率ｅと近地点引数ωの代わりにｅベクトルの表現を用いる。なお、ｅベクトル＝（ξ，η）＝（ｅ・cosω，ｅ・sinω）である。
【０００５】
軌道制御とは、人工衛星に搭載された推進系によって、所定の方向に速度ベクトルを付加することにより、軌道パラメータを変化させることである。また、どのような速度ベクトルを、いつ又はどこで付加すべきかを決定するのが軌道制御計画である。
【０００６】
ここで軌道面内制御に必要な軌道長半径ａの制御は、軌道１００２の接線（tangential）方向に速度ベクトルを付与する、すなわち加速又は減速することで可能であり、円軌道に近い場合には制御実施位置は特に制限されない。一方、ｅベクトルの制御は、制御実施位置と、付加する速度ベクトルの方向とを考慮して行わなければならない。図２に示すように、人工衛星が真近点離角ｆに位置するとすると、（１）方向に速度ベクトルΔＶを設定するとｅベクトルは（１ａ）方向に変化し、（２）方向に速度ベクトルΔＶを設定するとｅベクトルは（２ａ）方向に変化し、（３）方向に速度ベクトルΔＶを設定するとｅベクトルは（３ａ）方向に変化し、（４）方向に速度ベクトルΔＶを設定するとｅベクトルは（４ａ）方向に変化する。すなわち、付加する速度ベクトルΔＶ に対して、時計回りに９０°回転させた方向にｅベクトルは変化する。
【０００７】
準天頂軌道は１９９０年台に考案された軌道で以下の特徴を持つ。
・静止軌道と同じように軌道周期は地球自転と同期
・軌道傾斜角が日本の緯度と同じくらい
静止軌道では日本付近（北緯３５度）からの仰角はせいぜい５０°であるが、準天頂軌道は日本の天頂付近を通過するので、都市部のビルの谷間や山間部のように、周囲の視界が開けていない地域でも可視性が良いことを狙ったものである。当然、２４時間、四六時中天頂付近にあるわけではないので、複数衛星を配置することにより、次から次へと天頂付近に衛星が滞在するような複数衛星のコンステレーションを組まなければならない。
【０００８】
現在設計されている準天頂衛星の傾斜角は４０乃至５０°である。この軌道傾斜角のため、静止衛星と同じ地球自転同期軌道であっても、地上軌跡は南北と東西に変動し、図３に示すように、８の字を描く。図３（ａ）は、ｅ＝０すなわち円軌道の場合を示しており、図３（ｂ）は、ｅ＝０．０９９の場合を示しており、図３（ｃ）は、ｅ＝０．２８５の場合を示す。
【０００９】
準天頂衛星の軌道保持制御は、中心経度（図３に示すような８の字の交差点）と８の字の形を保持することを目的として実施される。中心経度は、軌道長半径が地球の非球状効果で増大することにより変動する。中心経度と軌道長半径が２次関数の関係になることを利用し、軌道長半径を適切な周期で減少させてやることにより、中心経度を所定範囲に保持することができる。すなわち、図４において、縦軸を軌道長半径ａとし、横軸を中心経度λとする。そして、制御前の軌道長半径ａbから、ΔＶだけ軌道長半径を減少させて軌道長半径の制御目標値ａtまで軌道長半径を変化させるものとする。その後準天頂衛星が制御無しで運用されると、ΔＴm後には、２次曲線１１０１を経て星印で示す位置まで軌道長半径及び中心経度の組み合わせは変化する。このような変化を前提として、中心経度の保持範囲内に中心経度λが収まるようにΔＶを決定しなければならない。この考え方は静止衛星の直下点経度保持と同じである。
【００１０】
８の字の形は軌道傾斜角とｅベクトルに依存するが、現在想定されている準天頂衛星では軌道傾斜角の変動はミッションに重大な支障を及ぼさないため、ｅベクトルの保持が問題になる。準天頂軌道のｅベクトルは、地球の非球状効果、太陽／月引力により極めて複雑な変動を示す。また、変動の大きさも静止軌道に比べると一桁以上大きい。この変動を打ち消すようなｅベクトルの補正を与えることにより、ｅベクトルを保持する。
【００１１】
ｅベクトルの制御は、制御実施場所と付加する速度ベクトルの方向によって、ｅベクトルが図３のように変化することに基づいて実施される。ここで、tangential方向ΔＶの方がradial方向ΔＶに比べてｅベクトルの変化は２倍、すなわち制御効率が良いことがポイントとなる。
【００１２】
従来、静止衛星等で行われている軌道面内保持制御は、軌道長半径ａを制御するためのtangential方向ΔＶでｅベクトルを同時に制御する方式である。軌道長半径ａの制御に必要なΔＶaとｅベクトル制御に必要なΔＶeの大きさは通常異なるので、両者の辻褄を合わせるため、２つのΔＶに分けて制御を実施する方式（２パートマヌーバ方式）がある。軌道上の制御実施位置はｅベクトルを補正すべき方向で決まる。２パートの場合は、軌道上の位相が１８０°離れた２地点で制御を実施することになる。
【００１３】
ところで準天頂軌道は、静止軌道と同じ周期（すなわち２４時間）で、日本上空通過時に天頂に位置するよう軌道傾斜角を４０°乃至５０°として設計された軌道である。離心率ｅは０．１程度で、アポジ点が天頂に位置する。準天頂衛星システムは、このような軌道の昇交点赤経Ωを等間隔に配した複数の人工衛星で構成され、３機構成の場合を図５に示す。図５は、３機構成の準天頂軌道を北極上空から見た図である。人工衛星ＳＡＴ１の昇交点赤経Ω1、人工衛星ＳＡＴ２の昇交点赤経Ω2及び人工衛星ＳＡＴ３の昇交点赤経Ω3は１２０°間隔で配置される。なお、矢印１０１４は地球１０００の自転方向である。
【００１４】
そして、人工衛星ＳＡＴ１の軌道１０１１のアポジ点通過時を中心とする８時間（真近点離角ｆについては１２０°から２４０°）であって太線で示された部分、人工衛星ＳＡＴ２の軌道１０１２のアポジ点通過時を中心とする８時間であって太線で示された部分、及び人工衛星ＳＡＴ３の軌道１０１３のアポジ点通過時を中心とする８時間であって太線で示された部分がそれぞれサービス時間帯（＝ミッション運用時間帯）である。
【００１５】
軌道保持制御は衛星のミッション（通信、測位、観測等）への影響が大きいので、軌道制御実施中はミッションを制約せざるを得ない場合がある。従って、準天頂軌道のサービス時間帯（φ：３０°〜１５０°）では、軌道制御実施を避けるべきである。
【００１６】
現状想定されている準天頂衛星では、軌道制御実施範囲は南半球に限定される。軌道制御実施前後の姿勢制御時間等も考慮すると、実際に軌道制御実施可能な範囲はφについて１８７°〜３５３°となる。
【００１７】
また、軌道面内保持制御の実施頻度は年２回程度が予定されているが、この間に実施されるアンローディング（アンローディングとは衛星姿勢をコントロールするホイール角運動量をリセットするために行うスラスタ噴射であり、軌道に影響を及ぼす場合が多い）による外乱ΔＶが無視できない大きさとなる。さらに、スラスタ取り付け方向のずれにより軌道面内制御時にtangential方向以外のΔＶが発生する。これらの軌道長半径、ｅベクトルへの影響を考慮して制御計画を作成する必要がある。
【００１８】
準天頂軌道の軌道長半径とｅベクトルを同時に制御するためには２パートマヌーバが必要なのだが、軌道制御実施可能範囲が１８０°に満たないため、従来の面内保持方式では対応できない。
【００１９】
準天頂衛星の軌道制御に関しては、特開２００５−９６７０２号公報には、軌道制御実施位置が限定された場合の軌道長半径ａとｅベクトルを補正する方式が開示されている。具体的には、radial方向ΔＶでｅベクトルを、tangential方向ΔＶでａを補正する制御を基本としている。しかしながら、上で述べたように、radial方向ΔＶでｅベクトルを補正するので、制御量が増加して効率が低くなってしまう。
【特許文献１】特開２００５−９６７０２号公報
【発明の開示】
【発明が解決しようとする課題】
【００２０】
従って、本発明の目的は、制御実施可能範囲が限定されている状況において、摂動や外乱の影響を考慮してできるだけ効率的な軌道面内保持制御計画を立案するための技術を提供することである。
【課題を解決するための手段】
【００２１】
本発明は、準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成方法であり、中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出するステップと、今回制御実施時における離心率ベクトル（ｅベクトル）を摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、離心率ベクトルの制御量を算出するステップと、軌道長半径目標制御量及び離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、軌道接線方向の制御量の初期値を算出するステップと、外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに軌道長半径の差分及び離心率ベクトルの差分を所定レベル以下になるように、軌道接線方向の制御量の初期値を補正する補正ステップとを含む。
【００２２】
このような処理を行うことによって、軌道接線方向に適切に制御を行うことができるようになる。すなわち、効率良く軌道面内保持制御を行うことができる。
【００２３】
上で述べた補正ステップが、アンローディング外乱の影響を考慮して次回制御実施時まで軌道伝播することにより実現中心経度を算出し、当該実現中心経度と目標中心経度との差が第１の所定レベル以下となるように、軌道接線方向の制御量の初期値を補正するステップと、準天頂衛星のスラスタ取り付け方向のずれの影響を考慮して次回制御時まで軌道伝播することにより実現軌道長半径及び実現離心率ベクトルを算出し、実現軌道長半径から得られる軌道長半径についての差分及び実現離心率ベクトルから得られる離心率ベクトルの差分が第２の所定レベル以下となるように、補正された軌道接線方向の制御量の初期値をさらに補正するステップとを含むようにしてもよい。このように２つのステップに分けることにより、適切に軌道接線方向の制御量を補正することができるようになる。
【００２４】
なお、上記方法は、コンピュータと当該コンピュータによって実行されるプログラムとの組み合わせにて実行される場合があり、当該プログラムは、例えばフレキシブルディスク、ＣＤ−ＲＯＭ、光磁気ディスク、半導体メモリ、ハードディスク等の記憶媒体又は記憶装置に格納される。また、ネットワークなどを介してデジタル信号として配信される場合もある。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。
【発明の効果】
【００２５】
本発明によれば、制御実施可能範囲が限定されている状況において、摂動や外乱の影響を考慮してできるだけ効率的な軌道面内保持制御計画を立案できるようになる。
【発明を実施するための最良の形態】
【００２６】
本発明の一実施の形態に係る、準天頂衛星の軌道面内保持制御計画装置の機能ブロック図を図６に示す。軌道面内保持制御計画装置１００は、周知の軌道決定装置２００が算出する初期軌道要素のデータを格納する初期軌道要素格納部２０１と、周知の制御管制系システム２０２が算出するアンローディングΔＶのデータを格納するアンローディングΔＶデータ格納部２０３とを参照するようになっている。また、軌道面内保持制御計画装置１００は、ユーザから所定のパラメータの値の入力を受け付けるデータ入力部１０２と、データ入力部１０２によって入力されたデータを格納する入力データ格納部１０３と、各種パラメータの値を格納するパラメータ格納部１０５と、初期軌道要素格納部２０１と入力データ格納部１０３とパラメータ格納部１０５とに格納されているデータを用いて軌道長半径目標値及び軌道長半径制御量を算出する軌道長半径算出処理部１０１と、軌道長半径算出処理部１０１によって算出されたデータを格納する軌道長半径データ格納部１０４と、アンローディングΔＶデータ格納部２０３などのデータを用いて周知の軌道伝播処理を実施する軌道伝播処理部１０６と、軌道伝播処理部１０６の処理結果を格納する軌道伝播結果データ格納部１０７と、初期軌道要素格納部２０１と入力データ格納部１０３と軌道伝播結果データ格納部１０７とに格納されたデータを用いてｅベクトルの制御目標値及び制御量を算出するｅベクトル算出処理部１０８と、ｅベクトル算出処理部１０８の処理結果を格納するｅベクトルデータ格納部１０９と、入力データ格納部１０３とｅベクトルデータ格納部１０９と軌道長半径データ格納部１０４とに格納されているデータを用いて３パートマヌーバの計画値を算出する３パートマヌーバ計画値算出部１１０と、３パートマヌーバ計画値算出部１１０による処理結果を格納する計画値格納部１１１と、中心経度ずれによる制御目標値等の修正などを実施する第１補正処理部１１２と、スラスタ取り付け方向のずれによるピッチ・ヨー方向の速度ベクトルを考慮して制御目標値等の修正を行う第２補正処理部１１３と、計画値格納部１１１に格納される最終的な制御計画値などを出力する出力部１１４とを有する。
【００２７】
軌道伝播処理部１０６は、以下で述べる複数のケースの各々に応じた軌道伝播処理を実施する。従って、各々のケースに必要なデータを、図６に示した各データ格納部から読み出して軌道伝播処理を実施する。また、第１補正処理部１１２は、軌道長半径データ格納部１０４と軌道伝播結果データ格納部１０７とパラメータ格納部１０５と入力データ格納部１０３とに格納されたデータを用いて処理を行い、処理結果を軌道長半径データ格納部１０４に格納すると共に、必要な場合には３パートマヌーバ計画値算出部１１０に処理を実施するように命ずる。また、第２補正処理部１１３は、軌道長半径データ格納部１０４とｅベクトルデータ格納部１０９と軌道伝播結果データ格納部１０７とに格納されたデータを用いて処理を行い、処理結果を軌道長半径データ格納部１０４とｅベクトルデータ格納部１０９とに格納すると共に、必要な場合には３パートマヌーバ計画値算出部１１０に処理を実施するように命ずる。
【００２８】
次に、図７乃至図３０を用いて図６に示した軌道面内保持制御計画装置１００の処理内容について説明する。まず、軌道長半径算出処理部１０１は、初期軌道要素格納部２０１と入力データ格納部１０３とパラメータ格納部１０５とに格納されているデータを用いて、軌道長半径制御目標値及び軌道長半径目標制御量を算出し、軌道長半径データ格納部１０４に格納する（ステップＳ１）。
【００２９】
静止軌道と同様に、地球重力ポテンシャルの非球状成分により軌道長半径は永年摂動を受ける。この永年摂動によって東経１３５°付近では、軌道長半径は増加する。これにより中心経度が変化するが、軌道長半径と中心経度がほぼ２次式の関係であることから、静止衛星の東西方向保持と同様の考え方に基づき、一定期間毎に減速制御を行うことにより中心経度を保持する。以下では、次回制御までの期間ΔＴmと目標中心経度λtから、ΔＴm後に中心経度がλtとなるように制御目標軌道長半径ａt及び軌道長半径目標制御量Δａvを算出する。
【００３０】
本実施の形態では、初期軌道要素格納部２０１には、制御実施日の初期軌道長半径ａbが格納されている。また、ユーザは、データ入力部１０２に対して、次回制御までの予定期間ΔＴm、ΔＴm後の目標中心経度λt、及び摂動モデルから算出される軌道長半径変化率ａｄｏｔを入力し、これらのデータは入力データ格納部１０３に格納される。さらに、パラメータ格納部１０５には、ミッション軌道における平均軌道長半径ａ0、及びミッション軌道における平均運動ｎ0が格納されている。
【００３１】
ここで、任意時刻における経度λと軌道長半径ａとの関係は、以下のとおりである。
λ＝−（３ｎ0／４ａ0・ａｄｏｔ）Δａ2＋λ0        （１）
ここで、Ｋa＝３ｎ0／４ａ0・ａｄｏｔ  （２）
とする。また、Δａは平均軌道長半径である基準値ａ0からのずれで、Δａ＝０の時、λ＝λ0である。
【００３２】
制御実施直後の状態（軌道長半径制御量，中心経度）は、
（Δａt，λb）  （３）
であり、ΔＴm後の状態は、
（Δａ(ΔＴm)，λt）  （４）
とする。但し、Δａ(ΔＴm)＝Δａt＋ａｄｏｔ×ΔＴm    （５）
とする。
【００３３】
このとき、（３）式及び（４）式を（１）式に入力すると、以下のようになる。
λb＝−Ｋa・Δａt2＋λ0    （６）
λt＝−Ｋa・（Δａt＋ａｄｏｔ×ΔＴm）2＋λ0    （７）
（６）式及び（７）式をΔａtについて解くと以下のようになる。
Δａt＝｛（λb−λt）／Ｋa−（ａｄｏｔ×ΔＴm）2／２ａｄｏｔ・ΔＴm      （８）
ここで（λb−λt）では、λ0はキャンセルされるので必要ない。
【００３４】
よって、（２）式及び（８）式並びに以下で示す（９）式及び（１０）式で、軌道長半径の制御目標値ａtと、軌道長半径の目標制御量Δａvとが算出される。
ａt＝ａ0＋Δａt    （９）
Δａv＝ａt−ａb    （１０）
【００３５】
次に、ｅベクトル算出処理部１０８等は、ｅベクトル制御目標値ｅt（＝（ｅtx，ｅty））及びｅベクトル目標制御量Δｅ_v（＝（Δｅvx，Δｅvy）。なお、「_」は、ｅの下に付くアンダーバーを表し、ベクトルを意味する。以下同じ。）算出処理を実施する（ステップＳ３）。この処理については、図８Ａ及び図８Ｂを用いて説明する。
【００３６】
まず、軌道伝播処理部１０６は、制御実施日の初期軌道要素であるｅベクトルｅ_bを次回制御予定時期（すなわちΔＴm後）まで伝播する処理を実施する（ステップＳ２１）。ｅベクトルは、地球形状偏平の効果、太陽／月引力の摂動で変動する。変動の大きさはΩに依存し、且つ太陽／月引力の１８．６年周期の変化にも影響を受ける。従って、ミッション期間中のｅベクトルについては、その時々の摂動量を計算しなければならない。本ステップでは、制御実施日から次回制御予定期間ΔＴm後まで、軌道伝播処理部１０６によりアンローディングΔＶデータ格納部２０３に格納されているデータを用いて周知の特別摂動で軌道伝播し、摂動によるｅベクトル変動量Δｅ_p（近地点引数の変動量Δωp）を算出し、軌道伝播結果データ格納部１０７に格納する。
【００３７】
図に示すと図８Ｂのようになる。すなわち、制御実施日のｅベクトルｅ_bが摂動によるｅベクトル変動量Δｅ_pだけ動いてしまい、ｅベクトルｅ_bpとなってしまう。この際、ｅ_bからｅ_bpまで、Δωpだけ回転が生ずる。
【００３８】
次に、ｅベクトル算出処理部１０８は、ｅベクトルモードがユーザによって指示されているか判断する（ステップＳ２３）。入力データ格納部１０３に格納されているＩＦＥフラグの値が、「０」であればｅベクトルモードであり、「１」であればωモードであると判断する。ｅベクトルモードは、Δｅ_pの中心が指定されたｅベクトル保持中心ｅc（入力データ格納部１０３に格納されているデータ）になるようにｅベクトルの制御目標を決定するモードである。一方、ωモードは、Δωpの中心が、指定ω保持中心ωc（入力データ格納部１０３に格納されているデータ）になるようにｅベクトルの制御目標を決定するモードである。準天頂衛星軌道のΩが１８０°付近の場合にはｅベクトルの摂動はωの回転が主で、ｅベクトルの大きさの変化が少ないため、ωモードでも十分な場合がある。
【００３９】
図８Ｂに示すように、ｅベクトルモードの場合には、ｅベクトルの保持中心ｅcとΔｅ_p／２とが一致するように、制御目標のｅベクトルｅ_t0を決定する。同様にωモードの場合には、摂動による回転の１／２であるΔωp／２とω保持中心ωcとが一致するように、制御目標のｅベクトルｅt1を決定する。
【００４０】
そして、ｅベクトルモードであれば、ｅベクトル算出処理部１０８は、ｅベクトルの変化分Δｅ_p（＝ｅ_bp−ｅ_b）を基に、ｅベクトル制御目標値ｅt（＝（ｅtx，ｅty））及びｅベクトル目標制御量Δｅ_v（＝（Δｅvx，Δｅvy））を算出し、ｅベクトルデータ格納部１０９に格納する（ステップＳ２５）。
【００４１】
ステップＳ２１ではｅ_bpが算出されるので、ｅベクトルの変化分Δｅ_pは上で述べたように定義される。従って、図８Ｂから分かるように、以下のようにｅベクトル制御目標値ｅt（＝（ｅtx，ｅty））及びｅベクトル目標制御量Δｅ_v（＝（Δｅvx，Δｅvy））は算出される。
ｅ_t＝ｅ_c−Δｅ_p／２    （１１）
Δｅ_v＝ｅ_t−ｅ_b        （１２）
その後、図７の処理に戻る。
【００４２】
一方、ωモードであれば、ｅベクトル算出処理部１０８は、回転角Δωpに基づきｅベクトル制御目標値ｅt（＝（ｅtx，ｅty））及びｅベクトル目標制御量Δｅ_v（＝（Δｅvx，Δｅvy））を算出し、ｅベクトルデータ格納部１０９に格納する（ステップＳ２７）。
【００４３】
初期軌道要素である初期ｅベクトルｅ_b及びΔＴm後のｅベクトルｅ_bpとから、Δωpは以下のように表される。
Δωp＝ａｒｇ（ｅ_bp）−ａｒｇ（ｅ_b）    （１３）
なお、ａｒｇ（Ｐ_）は、ベクトルＰ_を引数としてＰ_の角度を出力する関数であり、以下のように算出される。ＡＴＡＮはアークタンジェントである。
ａｒｇ（Ｐ_）＝ａｒｇ（Ｘ，Ｙ）＝ＡＴＡＮ（−Ｘ，−Ｙ）＋π
ωモードでは、ｅベクトルの大きさは変化させないので、以下の式が成り立つ。
ωt＝ωc−Δωp／２    （１４）
ｅt＝ｅb                          （１５）
ｅ_t＝（ｅt・ｃｏｓωt，ｅt・ｓｉｎωt）      （１６）
ωモードでも、ｅベクトル目標制御量については（１２）式が適用される。
【００４４】
図７の説明に戻って、３パートマヌーバ計画値算出部１１０は、パラメータ格納部１０５と軌道長半径データ格納部１０４と入力データ格納部１０３とｅベクトルデータ格納部１０９とに格納されているデータを用いて、以下で述べる三元連立方程式を解き、３パートマヌーバの初期制御量（１パート目の制御量Ｖ１，２パート目の制御量Ｖ２，３パート目の制御量Ｖ３）を算出し、計画値格納部１１１に格納する（ステップＳ５）。
【００４５】
本実施の形態では、図９に示すように、準天頂衛星のサービス期間は実施制御実施不可範囲として指定されており、それ以外の範囲において緯度φ1、φ2及びφ3の３箇所で制御を軌道の接線方向（tangential方向）に実施する。すなわち、φ1でΔＶ1だけ制御を実施し、φ2でΔＶ2だけ制御を実施し、φ3でΔＶ3だけ制御を実施し、制御前ｅベクトルを３段階に渡って変化させて、制御目標ｅベクトルになるようにする。なお、増速方向を＋方向とする。
【００４６】
ここで、最初の制御実施時刻における緯度φ1は、入力データ格納部１０３に格納されており、最初の制御実施時刻から最後の制御実施時刻までの時間Δτf又は最後の制御実施時刻における緯度φ3も入力データ格納部１０３に格納されている。Δτfが指定されている場合には、衛星平均運動よりφ1からΔτf後の緯度φ3を周知の方法を用いて算出する。また、２回目の制御実施時刻における緯度φ2は、（φ1＋φ3）／２で算出する。なお、［φ1，φ3］＞１８０°の場合には、φ3＝φ1＋１８０°とする。
【００４７】
また、ミッション軌道における平均運動ｎ0及びミッション軌道における平均速度の大きさＶ0は、パラメータ格納部１０５に格納されている。さらに、軌道長半径目標制御量Δａvは、軌道長半径データ格納部１０４に格納されており、ｅベクトル目標制御量Δｅ_v＝（Δｅvx，Δｅvy）は、ｅベクトルデータ格納部１０９に格納されている。
【００４８】
以上のデータを用いて、軌道長半径方向の増速量Ｖa、ｅベクトルのｘ方向の増速量Ｖex、ｅベクトルのｙ方向の増速量Ｖeyを以下のように算出する。
Ｖa＝（ｎ0／２）Δａv      （１７）
Ｖex＝（Ｖ0／２）Δｅvx    （１８）
Ｖey＝（Ｖ0／２）Δｅvy    （１９）
【００４９】
また、緯度φiで増速制御量Ｖiを与えた場合、ｅベクトルのＸ及びＹ成分に寄与する増速量は、Ｖi・cosφi及びＶi・sinφiで与えられる。
【００５０】
従って、（１７）乃至（１９）式から、以下の式が成り立つ。
Ｖ1＋Ｖ2＋Ｖ3＝Ｖa        （２０）
Ｖ1・cosφ1＋Ｖ2・cosφ2＋Ｖ3・cosφ3＝Ｖex      （２１）
Ｖ2・sinφ1＋Ｖ2・sinφ2＋Ｖ3・sinφ3＝Ｖey      （２２）
【００５１】
３パートマヌーバ計画値算出部１１０は、（２０）乃至（２２）式をＶ1、Ｖ2及びＶ3の三元連立方程式として解き、結果をφ1、φ2及びφ3と共に、計画値格納部１１１に格納する。
【００５２】
次に、第１補正処理部１１２等は、中心経度ずれに基づく補正処理を実施する（ステップＳ７）。この中心経度ずれに基づく補正処理については、図１０及び図１１Ａを用いて説明する。
【００５３】
ステップＳ１では軌道長半径と中心経度との関係について２次形式を仮定して軌道長半径制御目標値及び軌道長半径目標制御量を算出していたが、アンローディングΔＶ等の効果により、上記関係は厳密には２次形式ではないので、実際にΔＴm後まで特別摂動で軌道伝播すると、目標値λtからずれてしまう。また、面内制御時のスラスタ噴射で、進行方向以外のピッチ・ヨー方向（ピッチ方向（紙面に垂直な方向）、ヨー方向（進行方向に垂直な方向）については図１０に示す方向）の速度成分ΔＶが発生する。本ステップは、前者についての補正を行う処理である。
【００５４】
まず、軌道伝播処理部１０６は、アンローディングΔＶデータ格納部２０３に格納されているアンローディング外乱ΔＶ及び計画値格納部１１１に格納されている３パートマヌーバのデータ等を考慮して、次回制御時（ΔＴm後）まで周知の軌道伝播処理を実施することにより、伝播後の中心経度λrを算出し、軌道伝播結果データ格納部１０７に格納する（図１１Ａ：ステップＳ３１）。この処理は周知であるからこれ以上述べない。
【００５５】
次に、第１補正処理部１１２は、軌道伝播結果データ格納部１０７に格納されている伝播後の中心経度λrと、入力データ格納部１０３に格納されている目標経度λtとの差δλt（＝λr−λt）を算出し、メインメモリなどの記憶装置に格納する（ステップＳ３３）。
【００５６】
そして、第１補正処理部１１２は、δλtが所定の閾値未満であるか否かで収束したか判断する（ステップＳ３５）。伝播後の中心経度λrと目標経度λtとの差δλtが十分小さくなっている場合には、適切な制御量並びに軌道長半径目標制御値及び軌道長半径目標制御量が算出されていることになるので、元の処理に戻る。通常ステップＳ３３を１回実施しただけでは収束したと判断されることはないが、収束したと判断された場合には、中心経度についてはステップＳ１乃至Ｓ５が適切に実施されたと判断されたことになる。
【００５７】
一方、収束したと判断されなかった場合、第１補正処理部１１２は、δλtを補正するための軌道長半径目標値の変化分δａｔを算出し、当該δａｔを用いて軌道長半径目標値ａt’及び軌道長半径目標制御量Δａv’を算出し、軌道長半径データ格納部１０４に格納する（ステップＳ３７）。
【００５８】
具体的には、中心経度λと軌道長半径ａとの関係は（１）式に示したとおりである。
λ＝−Ｋa・Δａ2＋λ0      （１）'
ここで、Ｋa＝３ｎ0／４ａ0・ａｄｏｔである。
【００５９】
ΔＴm後における軌道長半径は、以下のように表される。
Δａm＝ａt＋ａｄｏｔ・ΔＴm−ａ0    （２３）
制御直後の軌道長半径にδａの変動があった場合のΔＴm後の中心経度変動は以下のように表される。
δλ＝−２Ｋa・Δａm・δａ        （２４）
従って、中心経度の実現値と目標値との差分δλt（＝λr−λt）を補正するための軌道長半径目標値の変化分は、以下のとおりである。
δａt＝−（１／２Ｋa・Δａm）×（−δλt）
＝δλt／（２Ｋa・Δａm）        （２５）
【００６０】
従って、軌道長半径制御目標値の修正値ａt'と軌道長半径目標制御量の修正値Δａv’は以下のように算出される。
ａt'＝ａt＋δａt        （２６）
Δａv'＝Δａv＋δａt        （２７）
【００６１】
さらに、第１補正処理部１１２は、３パートマヌーバ計画値算出部１１０に処理を指示する。
【００６２】
この指示に応じて、３パートマヌーバ計画値算出部１１０は、ステップＳ３７で算出され且つ軌道長半径データ格納部１０４に格納されている軌道長半径目標制御量Δａv’及びｅベクトルデータ格納部１０９に格納されている従前のｅベクトル目標制御量Δｅ_vを用いて三元連立方程式（２０）乃至（２２）式を解き、３パートマヌーバの修正制御量Ｖ1'、Ｖ2'、Ｖ3'を算出して、計画値格納部１１１に格納する（ステップＳ３９）。そしてステップＳ３１に戻る。
【００６３】
以上のような処理を収束、すなわちδλtが十分小さくなるまで繰り返す。
【００６４】
図７の説明に戻って、第２補正処理部１１３等は、スラスタの影響を考慮した補正処理を実施する（ステップＳ９）。この補正処理については、図１１Ｂを用いて説明する。
【００６５】
本ステップは、上で述べたように、面内制御時のスラスタ噴射で、進行方向以外のピッチ及びヨー方向の速度成分ΔＶが発生するため、制御直後の軌道がａtｅtからずれてしまうのを補正する処理である。従って、目標値の修正値ａt’及びｅt’を算出する。
【００６６】
まず、軌道伝播処理部１０６は、パラメータ格納部１０５に格納されているスラスタ取り付け方向のずれによるピッチ及びヨー方向の速度成分ΔＶ及び計画値格納部１１１に格納されている３パートマヌーバのデータ等を考慮して、次回制御時（ΔＴm後）まで周知の軌道伝播処理を実施することにより、伝播後の軌道長半径実現値ａr及びｅベクトルの実現値ｅ_rを算出し、軌道伝播結果データ格納部１０７に格納する（図１１Ｂ：ステップＳ４１）。この処理は周知であるからこれ以上述べない。
【００６７】
次に、第２補正処理部１１３は、軌道伝播結果データ格納部１０７に格納されている伝播後の軌道長半径実現値ａrと、軌道長半径データ格納部１０４に格納されている軌道長半径制御目標値ａtとの差δａtと、軌道伝播結果データ格納部１０７に格納されている伝播後のｅベクトル制御実現値ｅ_rと、ｅベクトルデータ格納部１０９に格納されているｅベクトル制御目標値ｅ_tとの差δｅ_tを算出し、メインメモリなどの記憶装置に格納する（ステップＳ４３）。
【００６８】
そして、第２補正処理部１１３は、δａt及びδｅ_tが所定の閾値未満であるか否かで収束したか判断する（ステップＳ４５）。上記δａt及びδｅ_tが十分小さくなっている場合には、適切な制御量、軌道長半径制御目標値並びに軌道長半径目標制御量、ｅベクトル制御目標値及びｅベクトル目標制御量が算出されていることになるので、元の処理に戻る。通常ステップＳ４１を１回実施しただけでは収束したと判断されることはないが、収束したと判断された場合には、スラスタの噴射の影響については、ステップＳ１乃至Ｓ５が適切に処理されたと判断されたことになる。
【００６９】
一方、収束したと判断されなかった場合、第２補正処理部１１３は、δａtを用いて軌道長半径制御目標値ａt’及び軌道長半径目標制御量Δａv’を算出し、軌道長半径データ格納部１０４に格納すると共に、δ_ｅtを用いてｅベクトル制御目標値ｅ_t’及びｅベクトル目標制御量Δｅ_v’を算出し、ｅベクトルデータ格納部１０９に格納する（ステップＳ４７）。以下のように算出する。
ａt’＝ａt−δａt    （２８）  
Δａv’＝Δａt−δａt      （２９）
ｅ_t’＝ｅ_t−δｅ_t        （３０）
Δ_ｅv’＝Δ_ｅv−δｅ_t      （３１）
【００７０】
さらに、第２補正処理部１１３は、３パートマヌーバ計画値算出部１１０に処理を指示する。
【００７１】
この指示に応じて、３パートマヌーバ計画値算出部１１０は、ステップＳ４７で算出され且つ軌道長半径データ格納部１０４に格納されているた軌道長半径目標制御量Δａv’及びステップＳ４７で算出され且つｅベクトルデータ格納部１０９に格納されているｅベクトル目標制御量Δｅ_v’を用いて三元連立方程式（２０）乃至（２２）式を解き、３パートマヌーバの修正制御量Ｖ1'、Ｖ2'、Ｖ3'を算出して、計画値格納部１１１に格納する（ステップＳ４９）。そしてステップＳ４１に戻る。
【００７２】
以上のような処理をδａt及びδｅ_tが十分小さくなるまで繰り返す。
【００７３】
図７の説明に戻って、出力部１１４は、計画値格納部１１１に格納されているＶ1’、Ｖ2’及びＶ3’を、φ1、φ2及びφ3と共に表示装置や印刷装置などの出力装置に出力する（ステップＳ１１）。その後、実際に準天頂衛星の制御が行われる。
【００７４】
以上のような処理を実施することによって、非サービス期間において、軌道の接線方向に実施すべき制御量等を得ることができるようになる。特に、近似を用いて三元連立方程式の解に問題を落とし込んだ上で、軌道伝播によってアンローディングΔＶやスラスタ噴射のずれなどを考慮して制御量等に微調整を行うことによって、容易に制御量を算出することができるようになっている。また、軌道の接線方向に制御を実施するので、使用燃料が少なくて済むという効果もある。
【００７５】
例えば、図１２乃至図１７に、西暦２０１０年から２０２５年までの軌道保持について、制御実施頻度を年２回とした場合で、保持開始時点（２０１０年）におけるΩ＝９０°の場合の制御例を示す。なお、ＮＯＭ９０−Ｅはｅベクトルモードの場合、ＮＯＭ９０−Ｗはωモードの場合を示す。図１２乃至図１７において、縦軸は制御量を表し、横軸は西暦を表す。図１２には、ｅベクトルモードにおいて、軌道長半径ａの制御量Ｖa、ｅベクトルのｘ成分ｅxの制御量Ｖex及びｙ成分ｅyの制御量Ｖeyとが各年においてプロットされている。図１３には、ωモードにおいて、軌道長半径ａの制御量Ｖa、ｅベクトルのｘ成分ｅxの制御量Ｖex及びｙ成分ｅyの制御量Ｖeyとが、各年においてプロットされている。図１４には、ｅベクトルモードにおいて、３パートマヌーバの各パート制御量Ｖ1、Ｖ2及びＶ3と全制御量Ｖ123とが、各年においてプロットされている。図１５には、ωモードにおいて、３パートマヌーバの各パートの制御量Ｖ1、Ｖ2及びＶ3と全制御量Ｖ123とが各年においてプロットされている。図１６には、ｅベクトルモードにおいて、２０１０年から各年の１年の制御量Ｖyearと、その累積値である累積制御量Ｖaccumとが、各年においてプロットされている。図１７には、ωモードにおいて、２０１０年から各年の１年の制御量Ｖyearと、その累積値である累積制御量Ｖaccumとが、各年においてプロットされている。
【００７６】
このように上で述べた状況においては、必要な制御量は、単調に変化している。また、ωモードの方が累積制御量Ｖaccumがやや少なく済んでいることも分かる。
【００７７】
また、図１８乃至図２３に、西暦２０１０年から２０２５年までの軌道保持について、制御実施頻度を年２回とした場合で、保持開始時点（２０１０年）におけるΩ＝２１０°の場合の制御例を示す。なお、ＮＯＭ２１０−Ｅはｅベクトルモードの場合、ＮＯＭ２１０−Ｗはωモードの場合を示す。図１８乃至図２３において、縦軸は制御量を表し、横軸は西暦を表す。図１８には、ｅベクトルモードにおいて、軌道長半径ａの制御量Ｖa、ｅベクトルのｘ成分ｅxの制御量Ｖex及びｙ成分ｅyの制御量Ｖeyとが各年においてプロットされている。図１９には、ωモードにおいて、軌道長半径ａの制御量Ｖa、ｅベクトルのｘ成分ｅxの制御量Ｖex及びｙ成分ｅyの制御量Ｖeyとが、各年においてプロットされている。図２０には、ｅベクトルモードにおいて、３パートマヌーバの各パート制御量Ｖ1、Ｖ2及びＶ3と全制御量Ｖ123とが、各年においてプロットされている。図２１には、ωモードにおいて、３パートマヌーバの各パートの制御量Ｖ1、Ｖ2及びＶ3と全制御量Ｖ123とが各年においてプロットされている。図２２には、ｅベクトルモードにおいて、２０１０年から各年の１年の制御量Ｖyearと、その累積値である累積制御量Ｖaccumとが、各年においてプロットされている。図２３には、ωモードにおいて、２０１０年から各年の１年の制御量Ｖyearと、その累積値である累積制御量Ｖaccumとが、各年においてプロットされている。
【００７８】
このように上で述べた状況においては、必要な制御量は、大きく変化している部分もあるが、累積制御量ＶaccumについてはΩ＝９０°の場合に比して大幅に削減されていることが分かる。また、ωモードの方が累積制御量Ｖaccumが少なく済んでいることも分かる。
【００７９】
さらに、図２４乃至図２９に、西暦２０１０年から２０２５年までの軌道保持について、制御実施頻度を年２回とした場合で、保持開始時点（２０１０年）におけるΩ＝３３０°の場合の制御例を示す。なお、ＮＯＭ３３０−Ｅはｅベクトルモードの場合、ＮＯＭ３３０−Ｗはωモードの場合を示す。図２４乃至図２９において、縦軸は制御量を表し、横軸は西暦を表す。図２４には、ｅベクトルモードにおいて、軌道長半径ａの制御量Ｖa、ｅベクトルのｘ成分ｅxの制御量Ｖex及びｙ成分ｅyの制御量Ｖeyとが各年においてプロットされている。図２５には、ωモードにおいて、軌道長半径ａの制御量Ｖa、ｅベクトルのｘ成分ｅxの制御量Ｖex及びｙ成分ｅyの制御量Ｖeyとが、各年においてプロットされている。図２６には、ｅベクトルモードにおいて、３パートマヌーバの各パート制御量Ｖ1、Ｖ2及びＶ3と全制御量Ｖ123とが、各年においてプロットされている。図２７には、ωモードにおいて、３パートマヌーバの各パートの制御量Ｖ1、Ｖ2及びＶ3と全制御量Ｖ123とが各年においてプロットされている。図２８には、ｅベクトルモードにおいて、２０１０年から各年の１年の制御量Ｖyearと、その累積値である累積制御量Ｖaccumとが、各年においてプロットされている。図２９には、ωモードにおいて、２０１０年から各年の１年の制御量Ｖyearと、その累積値である累積制御量Ｖaccumとが、各年においてプロットされている。
【００８０】
このように上で述べた状況においては、必要な制御量は、ほぼ単調に変化する。また、累積制御量ＶaccumについてはΩ＝９０°の場合とほぼ同様であることも分かる。但し、ωモードの方が累積制御量Ｖaccumが少なく済んでいることも分かる。
【００８１】
さらに、図３０に他の手法との差についても示す。図３０では、２０１０年の年間保持制御量ΔＶ（ｍ／ｓｅｃ）を、本方式、制御実施位置に制約のない２パートマヌーバ方式（２パート）、及び特開２００５−９６７０２号公報記載の手法（radialΔＶ）のそれぞれについて算出した場合の数値を列挙する。なお、本方式と２パートマヌーバ方式の場合についてはｅベクトルモードとωモードの両方について算出している。２パートマヌーバ方式は、制御実施位置の制約がないので、最も効率的に制御を行うことができる方式であるが、本方式は、ｅベクトルモードでΩ＝２１０°の場合を除き、２パートマヌーバ方式とほぼ同じ制御量となっている。すなわち、本方式は、人工衛星の燃料効率を向上させている。特開２００５−９６７０２号公報と比べれば効率の良さは顕著である。
【００８２】
以上本発明の一実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、軌道面内保持制御計画装置１００は、ネットワークに接続されており、同様に当該ネットワークに接続されているクライアント端末からの指示に応じて上で述べたような処理を実施するようにしても良い。
【００８３】
また、図６に示したような機能ブロック図は一例であって、必ずしもプログラムモジュール構成とは一致しない場合がある。また、処理結果さえ同一であれば、処理ステップを同時に実行したり、順番を入れ替えたりしても良い場合がある。
【００８４】
なお、軌道面内保持制御計画装置１００は、コンピュータ装置であって、このコンピュータ装置においては、図３１に示すように、メモリ２５０１（記憶部）とＣＰＵ２５０３（処理部）とハードディスク・ドライブ（ＨＤＤ）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。ＯＳ及びＷｅｂブラウザを含むアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。必要に応じてＣＰＵ２５０３は、表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、必要な動作を行わせる。また、処理途中のデータについては、メモリ２５０１に格納され、必要があればＨＤＤ２５０５に格納される。このようなコンピュータは、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及び必要なアプリケーション・プログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【００８５】
（付記１）
準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成プログラムであって、
コンピュータに、
中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出するステップと、
今回制御実施時における離心率ベクトルを摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、前記離心率ベクトルの制御量を算出するステップと、
前記軌道長半径目標制御量及び前記離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各前記制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、前記軌道接線方向の制御量の初期値を算出するステップと、
外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに前記軌道長半径の差分及び前記離心率ベクトルの差分を所定レベル以下になるように、前記軌道接線方向の制御量の初期値を補正する補正ステップと、
を実行させるための準天頂衛星の制御計画生成プログラム。
【００８６】
（付記２）
前記補正ステップが、
アンローディング外乱の影響を考慮して次回制御実施時まで軌道伝播することにより実現中心経度を算出し、当該実現中心経度と前記目標中心経度との差が第１の所定レベル以下となるように、前記軌道接線方向の制御量の初期値を補正するステップと、
前記準天頂衛星のスラスタ取り付け方向のずれの影響を考慮して次回制御時まで軌道伝播することにより実現軌道長半径及び実現離心率ベクトルを算出し、前記実現軌道長半径から得られる前記軌道長半径についての差分及び前記実現離心率ベクトルから得られる前記離心率ベクトルの差分が第２の所定レベル以下となるように、補正された前記軌道接線方向の制御量の初期値をさらに補正するステップと、
を含む請求項１記載の準天頂衛星の制御計画生成プログラム。
【００８７】
（付記３）
準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成方法であって、
中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出するステップと、
今回制御実施時における離心率ベクトルを摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、前記離心率ベクトルの制御量を算出するステップと、
前記軌道長半径目標制御量及び前記離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各前記制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、前記軌道接線方向の制御量の初期値を算出するステップと、
外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに前記軌道長半径の差分及び前記離心率ベクトルの差分を所定レベル以下になるように、前記軌道接線方向の制御量の初期値を補正する補正ステップと、
を含み、コンピュータにより実行される準天頂衛星の制御計画生成方法。
【００８８】
（付記４）
準天頂衛星の軌道を所定の範囲に保持するための制御計画を生成する制御計画生成装置であって、
中心経度が次回制御実施時において指定の目標中心経度になるように、今回制御実施時点における軌道長半径に対する、軌道長半径目標制御量を近似的に算出する手段と、
今回制御実施時における離心率ベクトルを摂動の影響を考慮した形で次回制御実施時まで軌道伝播することにより得られる前記離心率ベクトルの変化に基づき、前記離心率ベクトルの制御量を算出する手段と、
前記軌道長半径目標制御量及び前記離心率ベクトルの制御量と、制御実施可能期間における３つの制御実施緯度と、各前記制御実施緯度における軌道接線方向の制御量との関係を表す３元連立方程式を解くことにより、前記軌道接線方向の制御量の初期値を算出する手段と、
外乱要素の影響を考慮して次回制御実施時まで軌道伝播することにより算出される中心経度ずれ並びに前記軌道長半径の差分及び前記離心率ベクトルの差分を所定レベル以下になるように、前記軌道接線方向の制御量の初期値を補正する補正手段と、
を有する準天頂衛星の制御計画生成装置。
【図面の簡単な説明】
【００８９】
【図１Ａ】人工衛星の軌道制御の前提を説明するための図である。
【図１Ｂ】人工衛星の軌道制御の前提を説明するための図である。
【図２】人工衛星の軌道制御の前提を説明するための図である。
【図３】人工衛星の軌道制御の前提を説明するための図である。
【図４】人工衛星の軌道制御の前提を説明するための図である。
【図５】人工衛星の軌道制御の前提を説明するための図である。
【図６】本発明の実施の形態に係る軌道面内保持制御計画装置の機能ブロック図である。
【図７】本発明の実施の形態におけるメインフローを示す図である。
【図８Ａ】ｅベクトル制御目標値及びｅベクトル目標制御量算出処理の処理フローを示す図である。
【図８Ｂ】ｅベクトルモード及びωモードについて説明するための図である。
【図９】ｅベクトルの変化について説明するための図である。
【図１０】ピッチ及びヨー方向を説明するための図である。
【図１１Ａ】中心経度ずれに基づく補正処理の処理フローを示す図である。
【図１１Ｂ】スラスタの影響を考慮した補正処理の処理フローを示す図である。
【図１２】Ω＝９０°でｅベクトルモードの場合に、ａ、ｅx及びｅyの制御に必要な制御量の推移を表す図である。
【図１３】Ω＝９０°でωモードの場合に、ａ、ｅx及びｅyの制御に必要な制御量の推移を表す図である。
【図１４】Ω＝９０°でｅベクトルモードの場合に、３パートマヌーバの各パートにおける制御量及び全制御量の推移を表す図である。
【図１５】Ω＝９０°でωモードの場合に、３パートマヌーバの各パートにおける制御量及び全制御量の推移を表す図である。
【図１６】Ω＝９０°でｅベクトルモードの場合に、年間制御量と累積制御量の推移を表す図である。
【図１７】Ω＝９０°でωモードの場合に、年間制御量と累積制御量の推移を表す図である。
【図１８】Ω＝２１０°でｅベクトルモードの場合に、ａ、ｅx及びｅyの制御に必要な制御量の推移を表す図である。
【図１９】Ω＝２１０°でωモードの場合に、ａ、ｅx及びｅyの制御に必要な制御量の推移を表す図である。
【図２０】Ω＝２１０°でｅベクトルモードの場合に、３パートマヌーバの各パートにおける制御量及び全制御量の推移を表す図である。
【図２１】Ω＝２１０°でωモードの場合に、３パートマヌーバの各パートにおける制御量及び全制御量の推移を表す図である。
【図２２】Ω＝２１０°でｅベクトルモードの場合に、年間制御量と累積制御量の推移を表す図である。
【図２３】Ω＝２１０°でωモードの場合に、年間制御量と累積制御量の推移を表す図である。
【図２４】Ω＝３３０°でｅベクトルモードの場合に、ａ、ｅx及びｅyの制御に必要な制御量の推移を表す図である。
【図２５】Ω＝３３０°でωモードの場合に、ａ、ｅx及びｅyの制御に必要な制御量の推移を表す図である。
【図２６】Ω＝３３０°でｅベクトルモードの場合に、３パートマヌーバの各パートにおける制御量及び全制御量の推移を表す図である。
【図２７】Ω＝３３０°でωモードの場合に、３パートマヌーバの各パートにおける制御量及び全制御量の推移を表す図である。
【図２８】Ω＝３３０°でｅベクトルモードの場合に、年間制御量と累積制御量の推移を表す図である。
【図２９】Ω＝３３０°でωモードの場合に、年間制御量と累積制御量の推移を表す図である。
【図３０】本発明の実施の形態の効果を説明するための図である。
【図３１】コンピュータの機能ブロック図である。
【符号の説明】
【００９０】
１００    軌道面内保持制御計画装置
１０１  軌道長半径算出処理部    １０２  データ入力部
１０３  入力データ格納部    １０４  軌道長半径データ格納部
１０５  パラメータ格納部    １０６  軌道伝播処理部
１０７  軌道伝播結果データ格納部    １０８  ｅベクトル算出処理部
１０９  ｅベクトルデータ格納部    １１０  ３パートマヌーバ計画値算出部
１１１  計画値格納部    １１２  第１補正処理部
１１３  第２補正処理部    １１４  出力部
２００  軌道決定装置    ２０１  初期軌道要素格納部
２０２  制御管制系システム    ２０３  アンローディングΔＶデータ格納部
【図１Ａ】
2009227019.tif 000003
【図１Ｂ】
2009227019.tif 000004
【図２】
2009227019.tif 000005
【図３】
2009227019.tif 000006
【図４】
2009227019.tif 000007
【図５】
2009227019.tif 000008
【図６】
2009227019.tif 000009
【図７】
2009227019.tif 000010
【図８Ａ】
2009227019.tif 000011
【図８Ｂ】
2009227019.tif 000012
【図９】
2009227019.tif 000013
【図１０】
2009227019.tif 000014
【図１１Ａ】
2009227019.tif 000015
【図１１Ｂ】
2009227019.tif 000016
【図１２】
2009227019.tif 000017
【図１３】
2009227019.tif 000018
【図１４】
2009227019.tif 000019
【図１５】
2009227019.tif 000020
【図１６】
2009227019.tif 000021
【図１７】
2009227019.tif 000022
【図１８】
2009227019.tif 000023
【図１９】
2009227019.tif 000024
【図２０】
2009227019.tif 000025
【図２１】
2009227019.tif 000026
【図２２】
2009227019.tif 000027
【図２３】
2009227019.tif 000028
【図２４】
2009227019.tif 000029
【図２５】
2009227019.tif 000030
【図２６】
2009227019.tif 000031
【図２７】
2009227019.tif 000032
【図２８】
2009227019.tif 000033
【図２９】
2009227019.tif 000034
【図３０】
2009227019.tif 000035
【図３１】
2009227019.tif 000036