特許 5977040 軌道位置推定方法、軌道位置推定装置及び軌道位置推定プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5977040

(24)【登録日】2016年7月29日

(45)【発行日】2016年8月24日

(54)【発明の名称】軌道位置推定方法、軌道位置推定装置及び軌道位置推定プログラム

(51)【国際特許分類】

   G01S 19/14 20100101AFI20160817BHJP

   G01S 19/42 20100101ALI20160817BHJP

   G01C 21/28 20060101ALI20160817BHJP

   B64G 1/10 20060101ALI20160817BHJP

【ＦＩ】

   G01S19/14

   G01S19/42

   G01C21/28

   B64G1/10

【請求項の数】7

【全頁数】14

(21)【出願番号】特願2012-32479(P2012-32479)

(22)【出願日】2012年2月17日

(65)【公開番号】特開2013-167606(P2013-167606A)

(43)【公開日】2013年8月29日

【審査請求日】2015年1月19日

(73)【特許権者】

【識別番号】301072650

【氏名又は名称】ＮＥＣスペーステクノロジー株式会社

(74)【代理人】

【識別番号】100109313

【弁理士】

【氏名又は名称】机 昌彦

(74)【代理人】

【識別番号】100124154

【弁理士】

【氏名又は名称】下坂 直樹

(72)【発明者】

【氏名】棚町 健彦

【審査官】 三田村 陽平

(56)【参考文献】

【文献】 特表２００２−５４２４７０（ＪＰ，Ａ）

【文献】 特開２００５−１０６７２３（ＪＰ，Ａ）

【文献】 特表２００８−５１３７７５（ＪＰ，Ａ）

【文献】 特開２００８−２４９６３９（ＪＰ，Ａ）

【文献】 特開２０１０−２６６４６８（ＪＰ，Ａ）

【文献】 特表２００２−５３０６５４（ＪＰ，Ａ）

【文献】 特開２００２−０３９７８５（ＪＰ，Ａ）

【文献】 特開２００３−０４３１２９（ＪＰ，Ａ）

【文献】 米国特許第０６０８９５０７（ＵＳ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｓ １９／００−１９／５５

Ｇ０１Ｃ ２１／００−２１／３６

Ｂ６４Ｇ １／００−９９／００

(57)【特許請求の範囲】

【請求項1】

移動体の軌道位置を表す軌道要素を推定する推定方法であって、
前記移動体の複数の変数を含む軌道要素の観測値を取得して、前記観測値から前記複数の変数のうちの第１の変数について非線形であり、時間伝搬される前の第１の変数に基づいて求められる非線形項を除去するステップと、
前記非線形項が除去された観測値を時間伝搬させるステップと、
前記時間伝搬された観測値に時間伝搬された第１の変数に基づいて求められた非線形項を加算して、前記移動体の軌道要素の推定値とするステップを含むことを特徴とする推定方法。

【請求項2】

前記観測値を時間伝搬させるステップにおいて、カルマンフィルタの推定方法を利用することを特徴とする請求項１に記載の推定方法。

【請求項3】

前記観測値が取得できないときに、前記観測値の取得を停止し、
前記移動体の軌道要素の推定に使用される数値を与えて、前記非線形項が除去された観測値を時間伝搬させることを特徴とする請求項１又は２に記載の推定方法。

【請求項4】

移動体の軌道位置を表す軌道要素を推定する推定装置であって、
前記移動体の複数の変数を含む軌道要素の観測値を取得して、前記観測値から前記複数の変数のうちの第１の変数について非線形であり、時間伝搬される前の第１の変数に基づいて求められる非線形項を除去する平均化部と、
前記非線形項が除去された観測値を時間伝搬させる時間伝搬部と、
前記時間伝搬された観測値に時間伝搬された第１の変数に基づいて求められる非線形項を加算して、前記移動体の軌道要素の推定値とする加算部とを含むことを特徴とする推定装置。

【請求項5】

前記時間伝搬部は、カルマンフィルタの推定方法を利用して、前記非線形項が除去され
た観測値を時間伝搬させることを特徴とする請求項４に記載の推定装置。

【請求項6】

前記位置情報が取得できないときに、前記移動体の軌道要素の推定に使用される数値を受信して前記時間伝搬部に送信する状態変数更新部をさらに備え、
前記平均化部は前記観測値の取得を停止し、
前記時間伝搬部は前記数値に基づいて、前記非線形項が除去された観測値を時間伝搬させることを特徴とする請求項４又は５に記載の推定装置。

【請求項7】

移動体の軌道位置を表す軌道要素を推定する軌道位置推定プログラムであって、
前記移動体の複数の変数を含む軌道要素の観測値を取得して、前記観測値から前記複数の変数のうちの第１の変数について非線形であり、時間伝搬される前の第１の変数に基づいて求められる非線形項を除去する平均化処理と、
前記非線形項が除去された観測値を時間伝搬させる時間伝搬処理と、
前記時間伝搬された観測値に時間伝搬された第１の変数に基づいて求められる非線形項を加算して、前記移動体の軌道要素の推定値とする加算処理をコンピュータに行わせることを特徴とする軌道位置推定プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、移動体の軌道位置を推定する方法、推定装置及び推定プログラムに関し、特に、宇宙機に搭載されて宇宙機の時間伝搬された軌道位置を、オンボード軌道モデルのオンボードパラメータを介して、推定する方法、推定装置及び推定プログラムに関する。

【背景技術】

【0002】

地上からの距離が大きい高度で飛行する宇宙機においては、地表付近を航行する飛行体に比較して、ノイズや信号の減衰などにより地上との通信環境が貧弱である。

【0003】

一方、宇宙機は地表に対して大きな速度で移動しており、地表面の所定の領域を観測するなどの用途に対しては、地表からの高度が大きいことも考慮して地表面に対する宇宙機の位置を高い精度で制御する必要がある。

【0004】

地上からの通信による軌道位置の制御方法の他に、宇宙機に搭載した装置が自機の位置情報を収集し、この位置情報に基づいて軌道位置を推定する方法がある。

【0005】

宇宙機のオンボード軌道の決定には、宇宙機の軌道位置のセンサ情報として、ＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）が用いられる。

【0006】

ＧＰＳ受信機の出力には、ランダム或いは短周期のノイズがある。また、電波障害などの外部からの要因により、ＧＰＳ受信機における処理に大きな誤差が生じることがある。これらは宇宙機のオンボード軌道の決定値の擾乱源の一つである。

【0007】

地球を指向する宇宙機においては指向制御に宇宙機の軌道位置情報が必要であるが、ＧＰＳ受信機の処理における誤差はまた宇宙機の軌道位置の誤差として表れるため、宇宙機の指向制御の精度劣化の要因の一つになっている。

【0008】

このように、地球観測衛星等の宇宙機には高精度の指向安定度が必要であるが、このような宇宙機においては、指向安定度の劣化を防止するために、ＧＰＳ受信機（ＧＰＳＲ：ＧＰＳ Ｒｅｃｅｉｖｅｒ）の出力は直接には使用されない。このような宇宙機に対しては、次のような処理が行われて制御に使用される。すなわち、ＧＰＳＲの出力データは地上に降ろされ平滑化処理を受ける。宇宙機は所定の期間内では軌道位置を予測するオンボード軌道モデルを搭載している。地上では、平滑化処理を受けたデータを元に宇宙機の軌道を精密に決定し、オンボード軌道モデルのパラメータを算出する。算出されたパラメータ値は地上から宇宙機にコマンド送信され、宇宙機はこのパラメータ値に基づいて軌道位置を決定する。

【0009】

この軌道決定方法においては、オンボード軌道決定値にランダム或いは短周期ノイズの影響がないことから、宇宙機の指向安定度の劣化を防止できる。また、ＧＰＳ衛星システムや地上からの電波干渉等の、外部からの要因による短期的なＧＰＳＲ出力異常の影響も回避できる。

【0010】

上記の方法はこれらの利点を享受できるが、その一方で、宇宙機の運用処理が煩雑であるという問題がある。すなわち、数日毎にオンボード軌道モデルのパラメータ値を地上で算出して宇宙機へコマンド送信し、モデルのパラメータ値を更新する作業が必要である。これらが、宇宙機の運用のコストを引き上げる要因の一つになっていた。

【0011】

一般的に、衛星の軌道運動は大局的には楕円運動である。衛星は、この楕円運動の軌道上を短周期で変動している。以下、前者の大局的な軌道運動を表す軌道要素を平均要素と呼び、後者の変動運動まで考慮した瞬時の（それぞれの時刻の）軌道位置を表す軌道要素を接触軌道要素と呼ぶ。

【0012】

衛星の軌道位置の伝搬をオンボードで算出する場合、前者の大局的な楕円運動である平均要素は、時刻を引数にして解析的に求められる。一方、後者の接触軌道要素については、求める時刻までのそれぞれの時刻における軌道の変動量を逐次積算する必要がある。

【0013】

この変動量は摂動的に算出されるが、衛星軌道位置を高精度に決定するためには、短期摂動項を考慮した接触軌道要素を用いて、すなわち、摂動計算における時間刻みを小さくして、積算される切片を求める必要がある。しかし、長時間の軌道伝搬に対しては、逐次積分演算により演算誤差（高次の摂動項の影響及び数値的な誤差）が蓄積されるという問題がある。また、摂動計算と積分計算に要する計算量が増大するため、オンボード計算への負荷が大きくなる。

【0014】

上記の接触軌道要素の問題を解決した、オンボード軌道モデルのパラメータの算出式の一例を以下に示す。

【0015】

【数1】

なお、上記の算出式において、各変数と各定数は以下の量を意味する。

【0016】

【数2】

オンボード軌道パラメータＰ_１乃至Ｐ_１５を正確に求めることができれば、比較的容易な演算で正確に衛星軌道位置を伝搬することができる。

【0017】

特許文献１は、衛星軌道位置を予想する自律的な軌道伝搬システム及び方法を開示する。位置を予想する処理において、衛星の位置に関して以前に受信した情報に基づいて、衛星の補正加速度を発生させ衛星の位置を決定する。

【0018】

特許文献２は、ＧＰＳ衛星信号を観測値として用いてカルマンフィルタの理論に基づく位置算出処理を行う位置算出方法及び位置算出装置を開示する。位置算出処理において、最新の観測値と最新の予測値との差で表される最新観測予測差値を算出し、この最新観測予測差値をもちいてカルマンフィルタの補正処理に使用する観測誤差行列を設定する。

【先行技術文献】

【特許文献】

【0019】

【特許文献1】特表２０１１−５０３５５４号公報

【特許文献2】特開２０１０−２６６４６８号公報

【発明の概要】

【発明が解決しようとする課題】

【0020】

上記に示されたオンボード軌道モデルに算出式における、オンボード軌道パラメータＰ_１乃至Ｐ_１５は、天体運動や太陽活動の影響で時間的に変化する。このため、これらのパラメータの値を定期的に更新する必要がある。たとえば、低極軌道の地球観測衛星の場合には、３日に一度程度の頻度で地上で軌道決定を行った結果を基にオンボード軌道パラメータを算出し、このパラメータの値を衛星にコマンドで送信している。

【0021】

従って、関連技術におけるオンボード軌道モデルには、軌道位置を高精度に維持するための運用コストが大きいという問題があった。

【0022】

特許文献１の開示する軌道伝搬システム及び方法は、軌道の大局的な運動の軌道を蓄積されたデータに基づいて補正するのみであり、接触軌道要素には対応していない。このため、接触軌道要素の積算による軌道位置の偏位を算出できず、軌道の精度を維持できない。

【0023】

特許文献２の開示する位置算出方法は、観測値と予測値との差にカルマンフィルタの理論を適用するものであり、この差について線形性を前提にしている。上記のオンボード軌道モデルに示すように、宇宙機の軌道位置の算出には三角関数などの非線形な項の寄与を考慮する必要がある。特許文献２の開示する位置算出方法は非線形項の寄与を考慮しないため、接触軌道要素の寄与が蓄積した場合の軌道の位置を精度よく評価することができない。

【0024】

本発明は上述した点に鑑みてなされたもので、宇宙機の指向安定度の劣化を防ぎながら運用の煩雑性を解決するために、オンボード軌道モデルのパラメータをオンボードで自動的に推定しながらオンボード軌道モデルで軌道伝搬を行う、オンボード軌道モデルのオンボードパラメータ推定方法、推定装置及び推定プログラムを提供することを目的とする。

【課題を解決するための手段】

【0025】

上記の目的を達成するため、本発明のオンボードパラメータの推定方法は、移動体の軌道位置を表す軌道要素を推定する推定方法であって、前記移動体の複数の変数を含む軌道要素の観測値を取得して、前記観測値から前記複数の変数のうちの第１の変数について非線形であり、時間伝搬される前の第１の変数に基づいて求められる非線形項を除去するステップと、前記非線形項が除去された観測値を時間伝搬させるステップと、前記時間伝搬された観測値に時間伝搬された第１の変数に基づいて求められた非線形項を加算して、前記移動体の軌道要素の推定値とするステップを含むことを特徴とする。

【0026】

また、本発明のオンボードパラメータの推定装置は、移動体の軌道位置を表す軌道要素を推定する推定装置であって、前記移動体の複数の変数を含む軌道要素の観測値を取得して、前記観測値から前記複数の変数のうちの第１の変数について非線形であり、時間伝搬される前の第１の変数に基づいて求められる非線形項を除去する平均化部と、前記非線形項が除去された観測値を時間伝搬させる時間伝搬部と、前記時間伝搬された観測値に時間伝搬された第１の変数に基づいて求められる非線形項を加算して、前記移動体の軌道要素の推定値とする加算部とを含むことを特徴とする。

【0027】

さらに、本発明のオンボードパラメータの推定プログラムは、移動体の軌道位置を表す軌道要素を推定する軌道位置推定プログラムであって、前記移動体の複数の変数を含む軌道要素の観測値を取得して、前記観測値から前記複数の変数のうちの第１の変数について非線形であり、時間伝搬される前の第１の変数に基づいて求められる非線形項を除去する平均化処理と、前記非線形項が除去された観測値を時間伝搬させる時間伝搬処理と、前記時間伝搬された観測値に時間伝搬された第１の変数に基づいて求められる非線形項を加算して、前記移動体の軌道要素の推定値とする加算処理をコンピュータに行わせることを特徴とする。

【発明の効果】

【0028】

本発明によれば、定期的な地上からのコマンド送信により設定されていたオンボード軌道モデルのパラメータがオンボードで自動的に推定されるので、地上からのコマンド送信が不要になる。

【図面の簡単な説明】

【0029】

【図1】本発明の実施形態による軌道位置推定システムの構成と、軌道推定処理におけるデータの流れの一例を示す。

【図2】本発明の実施形態による軌道位置推定方法の処理の流れを示すフローチャートである。

【図3】本発明の実施形態による軌道位置推定システム構成と、ＧＰＳＲのデータ出力に問題が検出された場合の、軌道推定処理におけるデータの流れの一例を示す。

【発明を実施するための形態】

【0030】

本発明の実施形態を図面に基づいて説明する。ただし本発明は以下に示す実施形態に限定されない。

【0031】

宇宙機の軌道位置の伝搬の算出には、一般にカルマンフィルタを用いた推定方法が用いられる。カルマンフィルタは、観測値を元に誤差を考慮して導出したシステムの状態変数に対してその時間的変化を推定する線形モデルの一つとして広く用いられている。

【0032】

関連技術におけるオンボード軌道推定モデルのパラメータ算出式は、上記の算出式に示されるように三角関数を含む非線形な方程式であり、カルマンフィルタに基づく推定方法をそのまま用いることはできない。

【0033】

関連技術においては、オンボード軌道モデルの特徴として平均緯度引数αの平均的な値であるα（ｂａｒ）を引数とした三角関数を含む項が存在する。この項は地球の重力ポテンシャルに依存し、地球周回の定数倍で振動する短期摂動項と呼ばれる。

【0034】

本発明の実施形態に係る推定方法におけるオンボードパラメータの算出式を以下に説明する。本算出式は、状態方程式と観測方程式を含む。

【0035】

状態方程式を以下に示す。

【0036】

【数3】

上記に示した状態方程式は三角関数を含まず、線形項から構成される。この状態方程式に、カルマンフィルタによる推定方法が適用される。この推定処理により、パラメータ値が仮推定される。

【0037】

本実施形態における算出式に含まれる状態方程式は、関連技術におけるオンボード軌道モデル式から、数日間程度の時間では軌道伝搬の精度に与える影響が小さいと考えられる長周期の摂動項を除外したものである。すなわち、（θ_Ｇ−Ｐ_７）・ｔを引数に含む三角関数に比例する項を、以下に述べる時間伝搬処理の期間中では値の変動が小さいとして時間依存性を無視し定数項とみなす。長周期の摂動項は値が小さいとして無視するか、或いは状態変数に繰り込むものとする。

【0038】

さらに、このオンボード軌道モデル式から短期摂動項を除外して状態変数に対応させ、状態方程式に与える。本発明の実施形態に係る状態方程式は、接触要素の軌道６要素である、平均緯度引数α、軌道長半径ａ、昇交点赤経Ω、離心率ベクトル（ｅ_ｘ、ｅ_ｙ）、軌道傾斜角ｉに対して、平均要素に相当する部分、すなわちそれぞれの符号の上部にバーを付けた要素、α（ｂａｒ）、ａ（ｂａｒ）、Ω（ｂａｒ）、ｅ（ｂａｒ）_ｘ、ｅ（ｂａｒ）_ｙ、ｉ（ｂａｒ）を推定する方程式である。

【0039】

本実施形態における状態方程式において、各状態変数は時刻ｔの２次式で表される。それぞれの時間に関する微分を符号の上部に点（ｄｏｔ）をつけて表記する。さらに、時間に関する２階微分は符号の上部に点を２つ（ｄｏｕｂｌｅｄｏｔ）つけて表記する。

【0040】

本実施形態に係る状態方程式においては、上記に示されるように離心率ベクトル（ｅ（ｂａｒ）_ｘ、ｅ（ｂａｒ）_ｙ）と軌道傾斜角ｉ（ｂａｒ）は時刻ｔを含まず時間微分がゼロである。また、軌道長半径ａ（ｂａｒ）と昇交点赤経Ω（ｂａｒ）は時刻ｔの一次式であり２階微分がゼロになる。たとえば、軌道長半径ａ（ｂａｒ）とその１階微分ａ（ｂａｒ−ｄｏｔ）を２つの成分とする２行１列の行列を考えれば、その行列の時間微分はａ（ｂａｒ−ｄｏｔ）とゼロを成分とする行列である。すなわち、ａ（ｂａｒ）とａ（ｂａｒ−ｄｏｔ）を２つの成分とする行列を変数とする線形の微分方程式として表現される。

【0041】

平均緯度引数α（ｂａｒ）は時刻ｔの２次式であるので、３階微分がゼロになる。上記と同様にして、α（ｂａｒ）、その１階微分α（ｂａｒ−ｄｏｔ）、および２階微分α（ｂａｒ−ｄｏｕｂｌｅｄｏｔ）を３つの成分とする３行１列の行列を変数とする線形の微分方程式が得られる。

【0042】

これらの方程式をカルマンフィルタを用いた推定方法に従って数値的に解き、所定の期間における時間伝搬を得る。

【0043】

次に、本発明の実施形態に係る推定方法におけるオンボードパラメータの観測方程式を以下に示す。

【0044】

【数4】

上記の観測方程式においては、観測値から短期摂動項（αを引数とする三角関数を含む項）を差し引いた値が観測値として置き換えられる。これは状態空間中の状態変数（ｂａｒを付した要素）に対応付けられる。上記の状態方程式と観測方程式によりカルマンフィルタが構成される。すなわち、まず状態方程式により短期摂動項を取り除いた平均要素である、α（ｂａｒ）、ａ（ｂａｒ）、Ω（ｂａｒ）、ｅ（ｂａｒ）_ｘ、ｅ（ｂａｒ）_ｙ、ｉ（ｂａｒ）が推定される。次に、推定されたα（ｂａｒ）を用いて三角関数を含む短期摂動項が算出され、観測方程式においてこれを元に状態空間中の状態変数が観測値の空間中の観測値に対応付けられ、接触要素の軌道６要素である、平均緯度引数α、軌道長半径ａ、昇交点赤経Ω、離心率ベクトル（ｅ_ｘ、ｅ_ｙ）、軌道傾斜角ｉが算出される。

【0045】

短期摂動項は、α（ｂａｒ）を引数とする三角関数を含むが、α（ｂａｒ）の時間的変動が小さいことと、短期摂動項の値がそれぞれの観測量の平均要素の値に対して小さいことから、上記のように状態方程式において短期摂動項を無視すること、及び観測方程式において観測値から短期摂動項を差し引くという仮定はよい近似を与える。これにより、カルマンフィルタによる線形的な推定処理が可能である。

【0046】

本発明の実施形態にかかるオンボードパラメータ推定方法においては、ＧＰＳ受信機（ＧＰＳＲ）から出力される接触軌道６要素の観測値をもとに、上述の状態方程式と観測方程式により時間伝搬（プロパゲーション）と観測更新（アップデート）が行われる。

【0047】

ＧＰＳＲの出力にランダム或いは短周期ノイズが含まれる場合、カルマンフィルタの推定方法を適用する際に平均化によりノイズの影響が除去され、これにより指向安定度の劣化が抑制される。

【0048】

また、ＧＰＳＲの出力に短期的（数秒乃至数分間）に異常が発生した場合は、カルマンフィルタの観測更新を停止し時間伝搬のみを実施する。これにより、ＧＰＳＲの短期的出力異常に影響を受けないというシステムのロバスト性が維持される。

【0049】

さらに、ＧＰＳＲが長期間異常（故障等により）の場合は、カルマンフィルタの時間伝搬のみを実施し、従来のオンボード軌道モデルと同様に定期的に地上からオンボード軌道パラメータＰを更新する。これは、関連技術におけるオンボード軌道モデルと同等の軌道伝搬である。すなわち、ＧＰＳＲ故障時のバックアップの軌道モデルとしての利用が可能である。

【0050】

次に、上記の状態方程式と観測方程式における観測値の処理について、図１及び２を参照して具体的に説明する。

【0051】

本発明の実施形態にかかるオンボードパラメータ推定処理におけるデータの流れの一例を、概念図として図１に示す。

【0052】

図２は、本発明の実施形態にかかる軌道推定方法の処理の流れを示すフローチャートである。

【0053】

本発明の実施形態にかかるオンボード軌道モデルが稼働する、宇宙機に搭載される軌道位置推定システム１は、ＧＰＳ受信機（ＧＰＳＲ２０１）、オンボード計算機（搭載装置１０１）を含む。
搭載装置１０１は、図示しない記憶手段及び制御手段を備え、搭載ソフトウェアが記憶手段に格納される。制御手段は、搭載ソフトウェアが出力する指示に従って、各部の動作を制御する。

【0054】

ＧＰＳＲ２０１の出力は搭載装置１０１に入力され、入力された観測値は前処理部１０２において前処理される（ステップＳ１００１）。すなわち、ＧＰＳＲ２０１から出力される軌道データが、カルテジアン軌道要素の場合、または、離心率ｅ、近地点引数ω、真近点離角ｆを含む形式で表現されるケプラリアン軌道要素の場合は、これらのデータは離心率ベクトル（ｅ_ｘ、ｅ_ｙ）、平均緯度引数αで表現される形式のケプラリアン要素に変換される。

【0055】

次に、平均化部１０３において、前処理部１０２で得られた観測値から短期摂動項を差し引く（ステップＳ１００２）。短期摂動項は、一つ前の処理サイクルにおいて推定して得られた平均緯度引数α（ｂａｒ）の値を用いて求められる。短期摂動項を差し引いた平均 化した観測値が、カルマンフィルタ部１０４に入力される。

【0056】

カルマンフィルタ部１０４においては、平均緯度引数α（ｂａｒ）について非線形な短期摂動項を除いた状態変数の推定が行われる。

【0057】

カルマンフィルタ部１０４は観測更新部１０５と時間伝搬部１０６を含む。

【0058】

観測更新部１０５には、平均化部１０３から観測値が、短期摂動項を除いた平均化した値として入力される。観測方程式に従って、この観測値は観測値の空間から状態空間に対応付けられて、状態変数の値が更新される（ステップＳ１００３）。

【0059】

時間伝搬部１０６において、カルマンフィルタによる推定方法を利用して、状態方程式に従って、状態変数に対する時間伝搬が行われる（ステップＳ１００４）。

【0060】

時間伝搬部１０６において時間伝搬された平均緯度引数α（ｂａｒ）の推定値は、一次保持部１０７に一時的に保持されているが、カルマンフィルタ部１０４で算出された平均緯度引数α（ｂａｒ）の値により、一時保持部１０７に保持された値が更新される（ステップＳ１００５）。

【0061】

更新された平均緯度引数α（ｂａｒ）の値はまた平均化部１０３と加算部１０８に送信され、それぞれにおいて短期摂動項が算出される。平均化部１０３においては、上記の通り前処理部１０２で得られた観測値から短期摂動項が差し引かれて平均化した観測値が算出される。観測更新部１０５において、状態方程式により時間伝搬された状態変数が観測方程式により状態空間から観測値の空間に対応付けられ、時間伝搬され平均化された観測値が得られる。加算部１０８は、この時間伝搬され平均化された観測値のそれぞれに短期摂動項を加算して（ステップＳ１００６）、時間伝搬を行った観測値を接触軌道要素推定値として出力部１０９に最終的に出力する（ステップＳ１００７）。

【0062】

図３は、何らかの要因によって航法ステータスが示すＧＰＳＲ２０１データの信憑性が低いと検出された場合の、データ処理の流れの一例を概略的に示す。このとき、平均化部１０３と観測更新部１０５との間のデータ経路を一時的に遮断し、観測値の更新を停止して時間伝搬のみを行う。

【0063】

また、ＧＰＳＲ２０１が故障するなどして長時間観測値が入力されない場合は、カルマンフィルタ部１０４の中で時間伝搬部１０６の時間伝搬の処理のみを継続する。

【0064】

このとき、関連技術と同様にして、図３に示すように地上局３０１で算出したオンボード軌道パラメータＰ_１乃至Ｐ_９及びＰ_１５を定期的に衛星にコマンド送信してもよい。搭載装置１０１はオンボード軌道パラメータを受信し、状態変数更新部１１０において時間伝搬の状態変数、α（ｂａｒ）、ａ（ｂａｒ）、Ω（ｂａｒ）、ｅ（ｂａｒ）_ｘ、ｅ（ｂａｒ）_ｙ、ｉ（ｂａｒ）、α（ｂａｒ−ｄｏｔ）、α（ｂａｒ−ｄｏｕｂｌｅｄｏｔ）に対して本実施形態における状態方程式に従って変換し、時間伝搬部１０６における時間伝搬の中で使われる状態変数の更新を行う。

【0065】

観測値の時間伝搬をカルマンフィルタの推定方法を用いて推定する際に用いられる関係式を具体的に説明する。

【0066】

平均緯度引数αを算出するカルマンフィルタの推定式を以下に示す。

【0067】

【数5】

上記の推定式において、平均緯度引数αについての状態変数Ｘ_αを、α（ｂａｒ）、その１階微分α（ｂａｒ−ｄｏｔ）、および２階微分α（ｂａｒ−ｄｏｕｂｌｅｄｏｔ）を３つの成分とする３行１列の行列とする。

【0068】

状態方程式において、Ｆ_αはシステムの時間遷移に関する行列、Ｇ_ｗαは時間遷移に関するノイズである。この線形の微分方程式を解くことにより、状態変数α（ｂａｒ）の時間伝搬が得られる。

【0069】

また、観測方程式において、Ｈ_αは状態変数を観測値の空間に写像する行列であり、ｖ_αは観測に関するノイズである。すなわち、この観測方程式を解くことにより、観測できない状態変数Ｘ_αと観測量Ｙ_αとの関係を求め、状態方程式の解である状態変数Ｘ_αが観測値Ｙ_αと対応づけられる。

【0070】

軌道長半径ａを算出するカルマンフィルタの推定式を以下に示す。

【0071】

【数6】

上記の推定式において、軌道長半径ａについての状態変数Ｘ_ａを、ａ（ｂａｒ）、その１階微分ａ（ｂａｒ−ｄｏｔ）を２つの成分とする２行１列の行列とする。

【0072】

状態方程式において、Ｆ_ａはシステムの時間遷移に関する行列、Ｇ_ｗａは時間遷移に関するノイズである。この線形の微分方程式を解くことにより、状態変数ａ（ｂａｒ）の時間伝搬が得られる。

【0073】

観測方程式において、Ｈ_ａは状態変数を観測値の空間に写像する行列であり、ｖ_ａは観測に関するノイズである。すなわち、この観測方程式を解くことにより、観測できない状態変数Ｘ_ａと観測量Ｙ_ａとの関係を求め、状態方程式の解である状態変数Ｘ_ａが観測値Ｙ_ａと対応づけられる。

【0074】

観測方程式において、観測値の補正には平均緯度引数α（ｂａｒ）を引数とする項Ｂ_ａが含まれるが、この項には時間伝搬する直前の平均緯度引数α（ｂａｒ）の値を代入すればよい。

【0075】

軌道傾斜角ｉを算出するカルマンフィルタの推定式を以下に示す。

【0076】

【数7】

上記の推定式において、軌道傾斜角ｉについての状態変数Ｘ_ｉを、ｉ（ｂａｒ）を１つの成分とする１行１列の行列とする。

【0077】

状態方程式において、Ｆ_ｉはシステムの時間遷移に関する行列、Ｇ_ｗｉは時間遷移に関するノイズである。この線形の微分方程式を解くことにより、状態変数ｉ（ｂａｒ）の時間伝搬が得られる。

【0078】

観測方程式において、Ｈ_ｉは状態変数を観測値の空間に写像する行列であり、ｖ_ｉは観測に関するノイズである。すなわち、この観測方程式を解くことにより、観測できない状態変数Ｘ_ｉと観測量Ｙ_ｉとの関係を求め、状態方程式の解である状態変数Ｘ_ｉが観測値Ｙ_ｉと対応づけられる。

【0079】

観測方程式において、観測値の補正には平均緯度引数α（ｂａｒ）を引数とする項Ｂ_ｉが含まれるが、この項には時間伝搬する直前の平均緯度引数α（ｂａｒ）の値を代入すればよい。

【0080】

また、平均緯度引数α（ｂａｒ）の時間伝搬を求める関係式は、軌道長半径ａ（ｂａｒ）及び軌道傾斜角ｉ（ｂａｒ）をあらわには含まないので、観測方程式に代入する平均緯度引数α（ｂａｒ）は時間伝搬した直後の値でもよい。

【0081】

平均緯度引数αは、他の軌道要素に依存せず独立して推定できる。従って、まずは平均緯度引数の平均的な値α（ｂａｒ）を推定して、他の軌道要素はα（ｂａｒ）を外部からの入力として、それぞれ独立したカルマンフィルタを構成して推定することが可能である。昇交点赤経Ωは軌道長半径ａと同様に求められ、離心率ベクトル（ｅ_ｘ、ｅ_ｙ）は軌道傾斜角ｉ（ｂａｒ）と同様に求められる。

【0082】

以上、本発明の実施形態によれば、ＧＰＳＲの出力にランダム或いは短周期ノイズが含まれていても、カルマンフィルタにより平滑化が行われるので、指向安定度の劣化を防ぐことができる。

【0083】

また、ＧＰＳＲの出力が短期的（数秒乃至数分間）に異常になった場合は、カルマンフィルタの観測更新を止めて時間伝搬のみを実施することで、ＧＰＳＲの短期的出力異常に影響を受けないロバスト性を保つことができる。

【0084】

さらに、ＧＰＳＲが長期間異常（故障等）の場合には、カルマンフィルタの時間伝搬のみを実施し、従来のオンボード軌道モデルと同様に定期的に地上からオンボード軌道パラメータＰを更新することにより、従来のオンボード軌道モデルと同等の軌道伝搬を行うことも可能であり、ＧＰＳＲ故障時のバックアップ軌道モデルとしてもそのまま活用可能であり、ＧＰＳＲ故障時のバックアップ軌道モデルを別途搭載する必要がなくなる。

【0085】

本発明の実施形態に係る軌道位置推定方法においては、状態変数から平均緯度引数について非線形な項を差し引いて線形の状態方程式を得る。この線形の状態方程式にはカルマンフィルタの手法を適用することができ、観測誤差を考慮した時間伝搬を比較的小さな計算量で高い精度で求めることができる。時間伝搬した状態変数に対して、非線形項を摂動的に追加して、観測量に対応して軌道位置を表す接触要素が求められる。

【0086】

本実施形態に係る軌道位置推定方法においては、推定処理にかかる時間を短くして摂動項の寄与を積算して軌道位置を推定するので、軌道位置を表す接触要素から長周期の摂動項の時間的変動を無視できると仮定する。この近似により、状態方程式から非線形項が排除され、線形な状態方程式が導出される。

【0087】

以上、実施形態を参照して本発明を説明したが、本発明は上記各実施形態に限定されない。本発明の構成や詳細については当業者が理解し得るさまざまな変更を加えることができる。また、本発明には上記各実施形態の構成の一部又は全部を相互に適宜組み合わせたものも含まれる。

【産業上の利用可能性】

【0088】

本発明は上記の実施形態に限定されるものではなく、軌道位置の時間伝搬をオンボードで算出する移動体に好適に適用可能である。

【符号の説明】

【0089】

１ 軌道位置推定システム
１０１ 搭載装置
１０２ 前処理部
１０３ 平均化部
１０４ カルマンフィルタ部
１０５ 観測更新部
１０６ 時間伝搬部
１０７ 一時保持部
１０８ 加算部
１０９ 出力部
１１０ 状態変数更新部
２０１ ＧＰＳＲ
３０１ 地上局

【図1】

【図2】

【図3】