軌道予測装置、及び軌道予測方法
【課題】演算にかかるコストを削減できる軌道予測装置、及び軌道予測方法を提供することを目的としている。
【解決手段】測位衛星からの信号を受信して受信データを取得する情報取得部と、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出部と、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽から前記測位衛星への光の圧力である太陽輻射圧の影響を示す値、及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ部と、前記推定パラメータ部で推定された推定パラメータを用いて前記測位衛星の位置を推定する軌道予測部と、を備える。
【解決手段】測位衛星からの信号を受信して受信データを取得する情報取得部と、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出部と、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽から前記測位衛星への光の圧力である太陽輻射圧の影響を示す値、及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ部と、前記推定パラメータ部で推定された推定パラメータを用いて前記測位衛星の位置を推定する軌道予測部と、を備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、軌道予測装置、及び軌道予測方法に関する。
【背景技術】
【０００２】
ＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）は、位置測定、測量などに用いられている。地球上空の６つの軌道を２４機以上のＧＰＳ衛星が周回し、各ＧＰＳ衛星は、アルマナック（Ａｌｍａｎａｃ）データとエフェメリス（Ｅｐｈｅｍｅｒｉｓ）データの二種類の軌道情報を含む航法メッセージを送信している。
カーナビゲーション装置などのＧＰＳ受信機は、この航法メッセージを用いてＧＰＳ衛星からＧＰＳ受信機までの信号到達時間を測定する。そしてＧＰＳ受信機は、信号到達時間に光速を掛けることでＧＰＳ衛星とＧＰＳ受信機間の距離（擬似距離）を求めている。この作業を４機以上の衛星に対して行うことで、受信機の位置を三次元上で算出している。アルマナック・データは、軌道上における全ての衛星に関する軌道情報であり、１日に１回、データが更新される。エフェメリス・データには、各衛星の正確な位置情報と信号を発射した時刻情報が含まれており、２時間毎にデータが更新される。
【０００３】
ＧＰＳ受信機は、受信した各情報を自機の記憶部に、例えば電源を切る前に受信したアルマナック・データとエフェメリス・データを保存している。これらの情報には、各々有効期限がある。ＧＰＳ受信機の電源がオン状態で、ＧＰＳ衛星からの送信情報を送信毎に受信できていない場合など、ＧＰＳ受信機は、再度、送信情報を受信する必要がある。
アルマナック・データの再受信に要する時間は、例えば１２分３０秒であり、エフェメリス・データの再受信に要する時間は、３０秒である。
どちらかのデータの有効期限が切れにより再受信している期間、ＧＰＳ受信機は、位置測定を行えない。このため、エフェメリス・データの有効期限が切れた場合、ＧＰＳ受信機は、以前受信したエフェメリス・データを用いて、現在のＧＰＳ衛星の位置を予測し、現在のエフェメリス・データが利用できない場合に現在のＧＰＳ衛星の状態を使用して位置測定を行うことが提案されている（例えば、特許文献１参照）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】米国特許第７８３９３３０号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、特許文献１に記載の従来技術では、測位を行う際に用いる推定パラメータの個数が多いため、演算にかかるコストが高いという課題があった。
【０００６】
本発明は、上記の問題点に鑑みてなされたものであって、演算にかかるコストを削減できる軌道予測装置、及び軌道予測方法を提供することを目的としている。
【課題を解決するための手段】
【０００７】
上記目的を達成するため、本発明に係る軌道予測装置は、測位衛星からの信号を受信して受信データを取得する情報取得部と、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出部と、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽から前記測位衛星への光の圧力である太陽輻射圧の影響を示す値、及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ部と、前記推定パラメータ部で推定された推定パラメータを用いて前記測位衛星の位置を推定する軌道予測部と、を備えることを特徴としている。
【０００８】
また、本発明の軌道予測装置において、前記推定パラメータ部は、前記測位衛星の半径方向と進行方向及び軌道面垂直方向との成分による座標系に変換する行列を用いて前記推定パラメータを推定するようにしてもよい。
【０００９】
また、本発明の軌道予測装置において、前記推定パラメータ部は、前記測位衛星の位置の進行方向成分に対して予め算出された値を用い、前記測位衛星の速度の軌道面垂直方向成分に対して予め算出された値を用い、前記推定パラメータを推定するようにしてもよい。
【００１０】
また、本発明の軌道予測装置において、前記推定パラメータ部は、前記測位衛星に対する太陽輻射圧のワイバイアスとスケール係数とに対して各々予め算出した値を用い、前記推定パラメータを推定するようにしてもよい。
【００１１】
上記目的を達成するため、本発明は、軌道予測装置における軌道予測方法であって、情報取得部が、測位衛星からの信号を受信して受信データを取得する情報取得過程と、衛星位置速度算出部が、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出過程と、推定パラメータ部が、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽輻射圧の影響を示す値及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ過程と、軌道予測部が、前記推定パラメータ手順で推定された前記推定パラメータを用いて前記測位衛星の位置を推定する軌道予測過程と、を含むことを特徴としている。
【発明の効果】
【００１２】
本発明によれば、軌道における測位衛星の半径方向と進行方向及び軌道面垂直方向の各単位ベクトルと、位置速度情報とを用いて、測位衛星の位置と測位衛星の速度と太陽から測位衛星への光の圧力である太陽輻射圧の影響及び地球自転の影響を含む推定パラメータを算出するようにした。このように測位衛星の半径方向と進行方向及び軌道面垂直方向の各単位ベクトルを用いて推定パラメータを推定するため、推定するパラメータを削減することができる。この結果、最小二乗法で用いる行列を従来技術より小さくできるので、推定に必要なメモリ容量を削減できる。また、最小二乗法で用いる行列を小さくできるので、演算に係る時間を短くできる。
【図面の簡単な説明】
【００１３】
【図１】第１実施形態に係る軌道予測装置１０の概略ブロック図である。
【図２】同実施形態に係る過去のＧＰＳ衛星位置と速度データを記憶させる処理のフローチャートである。
【図３】同実施形態に係る衛星位置・速度記憶部１４０に記憶される情報と、削除される情報の一例を説明する図である。
【図４】保存間隔毎の予測誤差の一例を説明する図である。
【図５】同実施形態に係る軌道予測を行う処理手順のフローチャートである。
【図６】第２実施形態に係る衛星位置ｐ０（ｉ）の軌道計算の処理手順のフローチャートである。
【発明を実施するための形態】
【００１４】
［第１実施形態］
以下、図面を用いて、本発明の実施形態について説明する。
まず、一般的なＧＰＳ衛星の軌道予測の概略原理について説明する。
【００１５】
地球上空の６つの軌道を２４機以上のＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）衛星（以下、衛星という）が周回し、各衛星は、航法メッセージを送信している。航法メッセージは、フレームを単位に送信している。１フレームの送信には、３０秒を要する。１フレームは、５組のサブフレームで構成され、サブフレーム１から５が順番に送信される。また、サブフレーム１から３は、各衛星に固有の情報を含み、毎回同じ内容で構成される。サブフレーム４と５は、全衛星が同じ内容の情報で構成され、サブフレーム毎にページ１から２５で構成されている。このため、サブフレーム４と５の全ての情報を送信するには、２５フレームが必要である。このため、ＧＰＳ受信機が航法メッセージの全情報を得るには、１２分３０秒の時間がかかる。サブフレーム１は、各衛星の状態（正常に動作しているか否か）、衛星が送信する衛星のクロック誤差を補正するための係数であるクロック補正係数などにより構成されている。サブフレーム２は、各衛星の軌道情報（エフェメリス（Ｅｐｈｅｍｅｒｉｓ）・データ）１／２で構成されている。サブフレーム３は、各衛星の軌道情報（エフェメリス・データ）２／２で構成されている。サブフレーム４は、ＧＰＳ受信機が受信する信号が、電離層により遅延する量を補正するための係数である電離層遅延補正係数、ＧＰＳ時刻とＵＴＣ（協定世界時間；Universal Time, Coordinated））との関係を表す情報であるＵＴＣ関係情報、全衛星の軌道情報（アルマナック（Ａｌｍａｎａｃ）・データ）１／２などで構成されている。サブフレーム５は、全衛星の軌道情報（アルマナック・データ）２／２で構成されている。また、各サブフレームの先頭には、ＧＰＳ時刻を示す情報が含まれている。なお、ＧＰＳ時刻とは、１週間を単位として衛星で管理されている時間であり、毎週、日曜日の０時からの経過時間で表される情報である。
【００１６】
また、エフェメリス・データは、衛星の位置を計算するために必要な軌道の６要素（昇交点赤経、軌道傾斜角、近地点引数、軌道長半径、離心率、真近点角）のデータ、各補正値、及び軌道のエポック時刻ｔｏｅ（軌道の元期（ｅｐｈｅｍｅｒｉｓ ｒｅｆｅｒｅｎｃｅ ｔｉｍｅ））などにより構成されている。エフェメリス・データは、２時間毎に更新される。また、エフェメリス・データの有効期間は、２時間±２時間である。
アルマナック・データは、各衛星の軌道情報である。また、アルマナック・データは、約１日に１度の頻度で更新されている。また、アルマナック・データの有効期間は、エポック時刻±３日間である。
【００１７】
ＧＰＳ受信機は、受信したこれらのデータを記憶部に記憶する。記憶されているアルマナック・データ及びエフェメリス・データ共に有効な状態から、ＧＰＳ受信機が測位を開始することをホット・スタートという。また、エフェメリス・データが無効でありアルマナック・データのみ有効な状態から、ＧＰＳ受信機が測位を開始することをウォーム・スタートという。アルマナック・データ及びエフェメリス・データ共に無効な状態から、ＧＰＳ受信機が測位を開始することをコールド・スタートという。
そして、ホット・スタートの場合、ＧＰＳ受信機は、電源投入後に最初の位置情報が出力されるまでの時間は、記憶されているエフェメリス・データを用いて測位できるので数秒程度である。一方、ウォーム・スタートの場合、ＧＰＳ受信機は、エフェメリス・データを再取得する必要があるので、取得に３０秒以上の時間がかかり、さらに取得したエフェメリス・データの検査にも時間がかかる。このため、ウォーム・スタート場合、電源投入後に最初の位置情報が出力されるまでの時間は、３０秒以上である。さらに、コールド・スタートの場合、全ての航法メッセージを取得する必要があるので、測位信号の捕獲と航法メッセージの取得に数分〜１２分３０秒かかる。このため、電源投入後に最初の位置情報が出力されるまでの時間は、数分〜１２分３０秒以上である。
【００１８】
軌道予測とは、以前収集されたエフェメリス・データを利用し、未来の衛星位置を予測する技術である。軌道予測装置は、この予測した衛星位置を測位に用いることで、ウォーム・スタート時にエフェメリス・データの受信完了を待たずに測位開始することができるため、ＧＰＳ初期位置算出時間（ＴＴＦＦ）を短縮することができる。ＧＰＳ受信機が衛星からの情報取得後、電源オフ状態の期間が短ければ、以前取得したデータを用いて、衛星の位置を把握でききるため、データの再取得までの時間が短くなる。
そして、軌道予測装置は、エフェメリス・データから計算される衛星位置・速度をバックアップ・メモリに保存する。次に、軌道予測装置は、軌道予測計算に必要なデータが揃い次第、軌道予測パラメータを推定し、推定した軌道予測パラメータを使用して軌道を予測する。次に、軌道予測装置は、予測した軌道をバックアップ・メモリへ保存する。そして、軌道予測装置は、エフェメリス・データが無効の時、バックアップ・メモリに保存された予測した軌道を使用して衛星位置を計算する。
【００１９】
衛星の位置は、衛星の初期位置と初期速度が与えられ、衛星にかかる加速度が計算され、積分が行われることで、計算される。衛星にかかる加速度ベクトルａｓａｔは、地球重力定数をＧＭｅ、地球中心から衛星までの位置ベクトルをｒとし、地球重力の地球非球面による加速度ベクトルをａ１、太陽重力による加速度ベクトルをａ２、月重力による加速度ベクトルをａ３、太陽輻射圧による加速度ベクトルをａ４とすると、次式（１）のように表される。なお、太陽輻射圧とは、太陽から飛んでくる光子（Ｐｈｏｔｏｎ）がＧＰＳ衛星２の表面に衝突して発生する力である。
【００２０】
【数１】
【００２１】
一般的な軌道予測プログラムでは、式（１）の加速度ａｓａｔを数値積分により時間積分することで、ＧＰＳ衛星の位置と速度を算出する。用いる数値積分法は、例えば、ＲＫＦ４５（ルンゲクッタフェールベルグ法；Runge-Kutta-Fehlberg method）及びＣｏｗｅｌｌである。
【００２２】
次に、本実施形態に係る軌道予測装置の構成について、図１を用いて説明する。図１は、本実施形態に係るＧＰＳの概略ブロック図である。
図１に示すように、ＧＰＳは、測位装置１、ＧＰＳ衛星２から構成されている。測位装置（ＧＰＳ受信機）１は、軌道予測装置１０、測位算出部２０を備えている。なお、ＧＰＳ衛星２は、１機ではなく、２４機以上の衛星から構成される。
軌道予測装置１０は、情報取得部１１０、記憶部１１５、アルマナック取得部１２０、衛星型抽出部１２５、エフェメリス取得部１３０、衛星位置速度算出部１３５、衛星位置・速度記憶部１４０、演算開始判定部１４５、パラメータ推定部１５０、軌道予測演算部１５５、予測結果記憶部１６０を備えている。
【００２３】
情報取得部１１０は、ＧＰＳ衛星２が送信する航法メッセージから、ＵＴＣ関係情報、ＧＰＳ時刻、電離層遅延補正係数を取得し、取得したＧＰＳ時刻、ＵＴＣ関係情報と電離層遅延補正係数を記憶部１１５に記憶させる。
記憶部１１５には、ＧＰＳ時刻、ＵＴＣ関係情報と電離層遅延補正係数が記憶されている。
【００２４】
アルマナック取得部１２０（情報取得部）は、ＧＰＳ衛星２が送信したアルマナック・データを取得し、取得したアルマナック・データを衛星型抽出部１２５に出力する。
衛星型抽出部１２５は、アルマナック取得部１２０が出力するアルマナック・データから衛星型番号を抽出し、抽出した衛星型番号を衛星位置・速度記憶部１４０に記憶させる。なお、衛星型番号とは、Ｂｌｏｃｋ ＩＩＡ，Ｂｌｏｃｋ ＩＩＲ，Ｂｌｏｃｋ ＩＩＲ−Ｍ，Ｂｌｏｃｋ ＩＩＦなどのＧＰＳ衛星２の型を示す番号である。
【００２５】
エフェメリス取得部１３０（情報取得部）は、ＧＰＳ衛星２が送信するエフェメリス・データを取得し、取得したエフェメリス・データを衛星位置速度算出部１３５に出力する。
【００２６】
衛星位置速度算出部１３５は、エフェメリス取得部１３０が出力したエフェメリス・データから衛星の状態（Ｈｅａｌｔｈ）を抽出する。衛星位置速度算出部１３５は、抽出した衛星の状態が健康であるか否かを判別する。衛星の状態とは、衛星の状態を表すコードで、０以外は何らかの異常が衛星にあることを示している。
衛星の状態が健康であると判別した場合、衛星位置速度算出部１３５は、取得したエフェメリス・データをエポック時刻ｔｏｅに近い偶数正時ＴＯＥ２に変換する。衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に記憶されているＧＰＳ衛星２の位置と速度情報（以下、位置速度情報という）の内、算出した時刻ＴＯＥ２−３日より以前の情報を保存領域から削除する。時刻ＴＯＥ２−３日より以前の情報を保存領域から削除する理由は、エフェメリス・データを均等に３日分記憶させるためである。後述するように、このように均等に３日分、記憶させたエフェメリス・データを用いて軌道予測装置１０が軌道予測を行った場合、予測誤差が少なくなるためである。
衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されているか否かを判別する。衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていると判別した場合、後述するように、衛星位置速度算出部１３５は、ＴＯＥ２の位置速度情報を削除後、衛星位置・速度記憶部１４０の保存領域の空きが、２個分の位置速度情報を記憶できる領域になるまで情報を削除する。衛星位置速度算出部１３５は、時刻ＴＯＥ２と時刻ＴＯＥ２＋２時間の各時刻の各ＧＰＳ衛星２の各位置と各速度を算出し、算出した２つの時刻の位置と速度を位置速度情報として衛星位置・速度記憶部１４０に記憶させる。
衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていないと判別した場合、後述するように、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０の保存領域の空きが、３個分の位置速度情報を記憶できる領域になるまで情報を削除する。衛星位置速度算出部１３５は、時刻ＴＯＥ２、時刻ＴＯＥ−２時間、及び時刻ＴＯＥ２＋２時間の各時刻のＧＰＳ衛星２の各位置と各速度を算出し、算出した２つの時刻の位置と速度を位置速度情報として衛星位置・速度記憶部１４０に記憶させる。
【００２７】
衛星位置・速度記憶部１４０には、衛星型番号、各ＧＰＳ衛星２における各時刻の位置速度情報が記憶されている。
演算開始判定部１４５は、衛星位置・速度記憶部１４０から過去の位置速度情報、衛星型番号を読み出し、記憶部１１５からＧＰＳ時刻、ＵＴＣ関係情報及び電離層遅延補正係数を読み出す。演算開始判定部１４５は、読み出した情報に基づき、後述する演算開始条件を満たしているか否かを判定する。演算開始条件を満たしていると判別した場合、演算開始判定部１４５は、パラメータ推定部１５０に読み出した情報と判定結果を出力する。
【００２８】
パラメータ推定部１５０は、演算開始判定部１４５が出力する読み出された情報と判定結果を用いて、推定パラメータｓを推定する。パラメータ推定部１５０は、推定した推定パラメータを軌道予測演算部１５５に出力する。なお、推定パラメータｓは、後述するように１１個の成分を有するベクトルである。
軌道予測演算部（軌道予測部）１５５は、パラメータ推定部１５０が出力する推定パラメータｓを使用し、後述するように、例えば将来３日間２時間ごとの時刻ｔe（ｋ）（ｋ＝０，・・・，３６）での衛星位置ｐe（ｉ）（ｉ＝０，・・・，１１０）を軌道計算する。軌道予測演算部１５５は、算出した予測軌道値を予測結果記憶部１６０に記憶させる。
予測結果記憶部１６０には、予測軌道値が記憶されている。
測位算出部２０は、予測結果記憶部１６０に記憶されている予測軌道値を読み出し、読み出した予測軌道値を外挿して衛星位置を算出し、公知の手法で位置測位を行う。
【００２９】
なお、軌道予測装置１０は、航法メッセージを受信し、受信した航法メッセージを復調して、復調した航法メッセージから各データを抽出する。情報取得部１１０とアルマナック取得部１２０及びエフェメリス取得部１３０は、この抽出された各データを取得している。例えば、情報取得部１１０が航法メッセージを復調し、復調した航法メッセージから各データを抽出し、抽出した各情報を衛星型抽出部１２５及び衛星位置速度算出部１３５に出力し、記憶部１１５に書き込むようにしてもよい。また、記憶部１１５、衛星位置・速度記憶部１４０、及び予測結果記憶部１６０に分ける例を説明したが、これらの記憶部をまとめて１つの記憶部に記憶させるようにしてもよい。
【００３０】
軌道予測装置１０が行う処理について説明する。
まず、軌道予測装置１０は、航法データの取得を開始する。
次に、アルマナック取得部１２０は、ＧＰＳ衛星２が送信するアルマナック・データを取得し、取得したアルマナック・データを衛星型抽出部１２５に出力する。
次に、衛星型抽出部１２５は、アルマナック取得部１２０が出力するアルマナック・データから衛星型番号を抽出し、抽出した衛星型番号を衛星位置・速度記憶部１４０に記憶させる。
次に、情報取得部１１０は、ＧＰＳ衛星２が送信する航法メッセージのサブフレーム４のデータから、電離層遅延補正係数、ＧＰＳ時刻、及びＵＴＣ関係情報を取得し、取得した電離層遅延補正係数、ＧＰＳ時刻、及びＵＴＣ関係情報を記憶部１１５に記憶させる。
【００３１】
次に、過去のＧＰＳ衛星の位置と速度データを記憶させる処理について、図２を用いて説明する。図２は、本実施形態に係る過去のＧＰＳ衛星の位置と速度データを記憶させる処理のフローチャートである。
【００３２】
（ステップＳ１）エフェメリス取得部１３０は、ＧＰＳ衛星２からエフェメリス・データが送信されるのを待つ。ＧＰＳ衛星２からエフェメリス・データが送信された場合、エフェメリス取得部１３０は、ＧＰＳ衛星２が送信するエフェメリス・データを取得し、取得したエフェメリス・データを衛星位置速度算出部１３５に出力する。ステップＳ１終了後、ステップＳ２に進む。
【００３３】
（ステップＳ２）衛星位置速度算出部１３５は、エフェメリス取得部１３０が出力するエフェメリス・データから各衛星の状態を抽出する。衛星位置速度算出部１３５は、抽出した各衛星の状態が異常を示す値ではないか否か且つエポック時刻ｔｏｅが更新されたデータであるか否かを判別する。
各衛星の状態が異常を示す値ではない、すなわち健康であり且つエポック時刻ｔｏｅが更新されたデータであると判別した場合（ステップＳ２；Ｙｅｓ）、ステップＳ３に進む。各衛星の状態が異常を示す値である、または、エポック時刻ｔｏｅが更新されたデータではないと判別した場合（ステップＳ２；Ｎｏ）、ステップＳ１に戻る。
【００３４】
（ステップＳ３）各衛星の状態が異常を示す値ではない且つエポック時刻ｔｏｅが更新されたデータであると判別した場合、衛星位置速度算出部１３５は、取得したエフェメリス・データをエポック時刻ｔｏｅに近い偶数正時ＴＯＥ２に変換する。なお、取得したエフェメリス・データをエポック時刻ｔｏｅに近い偶数正時ＴＯＥ２に変換する理由は、ＧＰＳ衛星２から送信されたエポック時刻ｔｏｅが、２時間間隔ではなく、まれにずれることがあるためである。
時刻ＴＯＥ２を時刻ｔ（ｋ）（ｋ＝０，・・・，ｎ−１）と表現する。以降の処理上、位置はｒ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での位置のｘ成分がｒ（３ｋ）、位置のｙ成分がｒ（３ｋ＋１）、位置のz成分がｒ（３ｋ＋２）である。速度はｖ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での速度ｘ成分がｖ（３ｋ）、速度ｙ成分がｖ（３ｋ＋１）、速度ｚ成分がｖ（３ｋ＋２）である。
衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に記憶されている位置速度情報の内、算出した時刻ＴＯＥ２−３日より以前の情報を保存領域から削除する。ステップＳ３終了後、ステップＳ４に進む。
【００３５】
（ステップＳ４）衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されているか否かを判別する。
衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていると判別した場合（ステップＳ４；Ｙｅｓ）、ステップＳ５に進む。衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていないと判別した場合（ステップＳ４；Ｎｏ）、ステップＳ７に進む。なお、衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていない状態は、例えば、ＧＰＳ衛星２が測位装置１から見て遮蔽されていて、エフェメリス・データを受信できない場合、測位装置１の電源がオフ状態などである。
【００３６】
（ステップＳ５）衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていると判別した場合、衛星位置速度算出部１３５は、ＴＯＥ２の位置速度情報を削除後、衛星位置・速度記憶部１４０の保存領域の空きが、２個分の位置速度情報を記憶できる領域になるまで位置速度情報を削除する。衛星位置速度算出部１３５は、経過時間の古いものから、２ｎ時間（ｎは１以上の自然数）毎の隣り合う情報の内、新しい時刻の情報を削除する。なお、経過時間とは、衛星位置・速度記憶部１４０に１つも記憶されていない状態から位置速度情報の記憶を開始してからの時間の経過であり、図３における縦軸に示した時間である。
具体的には、２時間毎の隣り合う情報の内、新しい時刻の情報を削除する。２時間毎の隣り合う情報を比較して、情報を削除しても、空き領域を確保できない場合は、４時間毎に隣り合う情報の内、新しい方の情報を削除する。空き領域を確保できない場合、衛星位置速度算出部１３５は、６時間毎に隣り合う情報について比較、削除を行う。空き領域を確保できない場合、衛星位置速度算出部１３５は、８時間毎に隣り合う情報について比較、削除を行う。ステップＳ５終了後、ステップＳ６に進む。
【００３７】
（ステップＳ６）次に、衛星位置速度算出部１３５は、時刻ＴＯＥ２と時刻ＴＯＥ２＋２時間の各時刻のＧＰＳ衛星２の各位置と各速度を算出する。次に、衛星位置速度算出部１３５は、算出した２つの時刻の各位置速度情報を衛星位置・速度記憶部１４０に記憶させる。なお、記憶されるＧＰＳ衛星２の位置と速度はＥＣＥＦ（地球中心・地球固定）座標系である。ステップＳ６終了後、ステップＳ１に戻る。
【００３８】
（ステップＳ７）衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていないと判別した場合、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０の保存領域の空きが、３個分の位置速度情報を記憶できる領域になるまで位置速度情報を削除する。衛星位置速度算出部１３５は、経過時間の古いものから、２ｎ（ｎは１以上の自然数）時間毎の隣り合う情報の内、新しい時刻の情報を削除する。
具体的には、２時間毎の隣り合う情報を比較して、情報を削除しても、空き領域を確保できない場合は、４時間毎に隣り合う情報の内、新しい方の情報を削除する。空き領域を確保できない場合、衛星位置速度算出部１３５は、６時間毎に隣り合う情報について比較、削除を行う。空き領域を確保できない場合、衛星位置速度算出部１３５は、８時間毎に隣り合う情報について比較、削除を行う。ステップＳ７終了後、ステップＳ８に進む。
【００３９】
（ステップＳ８）次に、衛星位置速度算出部１３５は、時刻ＴＯＥ２、時刻ＴＯＥ２−２時間、及び時刻ＴＯＥ２＋２時間の各時刻のＧＰＳ衛星２の各位置と各速度を算出する。次に、衛星位置速度算出部１３５は、算出した３つの時刻の各位置速度情報を衛星位置・速度記憶部１４０に記憶させる。ステップＳ８終了後、ステップＳ１に戻る。
【００４０】
次に、衛星位置・速度記憶部１４０に記憶される位置速度情報の順番と、削除される順番の具体的な例について、図３を用いて説明する。図３は、本実施形態に係る衛星位置・速度記憶部１４０に記憶される位置速度情報と、削除される位置速度情報の一例を説明する図である。図３において、横軸は保存回数、縦軸は経過時間である。なお、経過時間において、２時間毎のデータであり、「２」は４時間を表している。また、保存回数とは、エフェメリス・データから、ＧＰＳ衛星２の位置と速度を計算し、衛星位置・速度記憶部１４０へ位置速度情報を保存した回数である。経過時間は、保存開始の時刻から各保存された情報の時刻までの経過時間である。
【００４１】
まず、保存回数１では、まだ時刻ＴＯＥ２の位置速度情報が衛星位置・速度記憶部１４０に記憶されていない。このため、図３に示すように、衛星位置速度算出部１３５は、経過時間０、１（２時間後）、及び２（４時間後）の位置速度情報を衛星位置・速度記憶部１４０に記憶する。この場合、衛星位置速度算出部１３５は、ステップＳ８で説明したように、時刻ＴＯＥ２（ｇ２）と、時刻ＴＯＥ２−２時間（ｇ１）、及び時刻ＴＯＥ２＋２時間（ｇ３）の３つの位置速度情報を衛星位置・速度記憶部１４０に記憶する。
保存回数２において、ステップＳ６で説明したように、衛星位置速度算出部１３５は、２つの位置速度情報が衛星位置・速度記憶部１４０に記憶する。この場合、衛星位置速度算出部１３５は、この時点で最新の位置速度情報である位置速度情報ｇ３を削除した後、時刻ＴＯＥ２の位置速度情報を位置速度情報ｇ３’として衛星位置・速度記憶部１４０に記憶させる。位置速度情報ｇ３を削除する理由は、保存回数２における時刻ＴＯＥ２時間は、保存回数１における時刻ＴＯＥ２＋２時間であるためである。そして、衛星位置速度算出部１３５は、時刻ＴＯＥ２＋２時間の位置速度情報を位置速度情報ｇ４として記憶させる。
同様に、保存回数３において、衛星位置速度算出部１３５は、位置速度情報ｇ４を削除した後、時刻ＴＯＥ２の位置速度情報を位置速度情報ｇ４’として記憶させる。そして、時刻ＴＯＥ２＋２時間の位置速度情報を、位置速度情報ｇ５として記憶する。
【００４２】
経過時間４から１０の期間、記憶されている情報がないのは、エフェメリス・データを軌道予測装置１０が取得できなかったことを表している。
保存回数４において、保存回数４における時刻ＴＯＥ２から直近の偶数正時、すなわち経過時間１０（２０時間）にＧＰＳ衛星２の位置と速度情報が衛星位置・速度記憶部１４０に記憶されていないため、ステップＳ８のように、衛星位置速度算出部１３５は、再び３つの時刻の位置速度情報を衛星位置・速度記憶部１４０に記憶させる。このような動作を繰り返していくと、保存回数６で、衛星位置・速度記憶部１４０に保存される位置速度情報が１０個になる。
【００４３】
保存回数６に達した後の衛星位置・速度記憶部１４０に記憶されている位置速度情報を、便宜的にｇ１１〜ｇ２０とする。
保存回数７において、ステップＳ５で説明したように、衛星位置速度算出部１３５は、２つの位置速度情報分の記憶領域を確保する。このため、まず、保存回数６までと同様に、衛星位置速度算出部１３５は、保存回数６で保存された位置速度情報ｇ２０を削除する。次に、衛星位置速度算出部１３５は、位置速度情報を２つ書き込める領域が確保できるまで、衛星位置・速度記憶部１４０に保存されている位置速度情報を削除する。
具体的には、衛星位置速度算出部１３５は、経過時間が一番古い２時間間隔のペアの位置速度情報ｇ１１とｇ１２を比較する。位置速度情報ｇ１１は、保存回数１で記憶された位置速度情報ｇ１であるので、保存回数１における時刻ＴＯＥ２−２時間の位置速度情報である。位置速度情報ｇ１２は、保存回数１で記憶された位置速度情報ｇ１であり、保存回数１における時刻ＴＯＥ２の位置速度情報である。位置速度情報ｇ１２の方が新しい情報であるので、衛星位置速度算出部１３５は、位置速度情報ｇ１２を削除する。このようにして、衛星位置速度算出部１３５は、２つの空き領域を確保する。
次に、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に、時刻ＴＯＥ２時間の位置速度情報を新たな位置速度情報ｇ２０として記憶させ、時刻ＴＯＥ２＋２時間の位置速度情報をｇ２１として記憶させる。
【００４４】
なお、本実施形態では、ステップＳ５及び保存回数７において、２つの位置速度情報を削除した後に、時刻ＴＯＥの位置速度情報と時刻ＴＯＥ＋２時間の位置速度情報を衛星位置・速度記憶部１４０に記憶させる例を説明した。しかしながら、１つの位置速度情報を削除した後に、時刻ＴＯＥ＋２時間の位置速度情報を記憶させ、時刻ＴＯＥの位置速度情報と上書きするようにしてもよい。具体的には、保存回数７において衛星位置速度算出部１３５は、位置速度情報ｇ１１とｇ１２を比較する。衛星位置速度算出部１３５は、比較結果に基づき、時刻の新しい方の位置速度情報ｇ１２を削除する。
次に、衛星位置速度算出部１４は、時刻ＴＯＥ２時間の位置速度情報を、位置速度情報ｇ２０に上書きし、時刻ＴＯＥ２＋２時間の位置速度情報をｇ２１として記憶させるようにしてもよい。
【００４５】
このように、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０の位置速度情報に対して削除と追加を繰り返す。
保存回数１３において、衛星位置・速度記憶部１４０に記憶されている位置速度情報を、便宜上ｇ２１〜ｇ３０とする。
保存回数１４において、時刻ＴＯＥ２の位置速度情報が衛星位置・速度記憶部１４０に記憶されていないため、３つの位置速度情報を削除する必要がある。この場合、衛星位置速度算出部１３５は、まず、経過時間が一番古い２時間間隔のペアの位置速度情報ｇ２８とｇ２９を比較し、比較結果に基づき、位置速度情報ｇ２９を削除する。次に、衛星位置速度算出部１３５は、経過時間が一番古い４時間間隔のペアの位置速度情報ｇ２４とｇ２５を比較し、比較結果に基づき、位置速度情報ｇ２５を削除する。次に、衛星位置速度算出部１３５は、経過時間が二番古い４時間間隔のペアの位置速度情報ｇ２６とｇ２７を比較し、比較結果に基づき、位置速度情報ｇ２７を削除する。
以上の処理により、保存回数１３で保存されていた位置速度情報から３個の位置速度情報の削除を完了する。
また、保存回数１５回において、衛星位置速度算出部１３５は、経過時間０の情報が、ステップＳ３で説明した時刻ＴＯＥ２−３日の位置速度情報であるため削除する。
【００４６】
次に、このように、衛星位置・速度記憶部１４０に記憶されている位置速度情報を、間引きながら保存する理由について説明する。図４は、保存間隔毎の予測誤差の一例を説明する図である。図４（ａ）と図４（ｂ）は、位置速度情報を４時間毎、均等に３日間分取得した例である。図４（ｃ）と図４（ｄ）は、位置速度情報を３日間分の内、前後にのみ取取得し、３日に渡って中抜けした例である。図４（ｅ）と図４（ｆ）は、位置速度情報を３日の内、２日間に渡って中抜けした例である。
図４（ａ）、図４（ｃ）及び図４（ｅ）において、横軸は予測誤差の最大値を表し、縦軸はその割合を表している。図４（ｂ）、図（ｄ）及び図（ｆ）において、横軸は時間を表し、白丸は位置速度情報を取得できなかった「中抜け」状態を表し、黒丸は位置速度情報を取得できた状態を表している。
４時間間隔の場合、このように取得された位置速度情報を用いて、ＧＰＳ衛星２の位置を予測すると、図４（ａ）に示すように、予測誤差の最大値は７［ｍ］以下である。中抜けが３日に渡る場合、このように取得された位置速度情報を用いて、ＧＰＳ衛星２の位置を予測すると、図４（ｃ）に示すように、予測誤差の最大値は２０［ｍ］以上もあり、その割合もランダムである。中抜けが２日に渡る場合、このように取得された位置速度情報を用いて、ＧＰＳ衛星２の位置を予測すると、図４（ｅ）に示すように、予測誤差の最大値は２０［ｍ］以上が５０％もある。
このように、予測誤差を抑えるには、中抜けの期間が長い状態は好ましくなく、位置速度情報が均等に保存されていることが好ましい。
【００４７】
このため、本実施形態では、単に１０個の位置速度情報を時系列に記憶させているのではなく、限られた記憶容量内で３日分の位置速度情報を時間間隔が均等に記憶するように、保存した位置速度情報を削除して、新たな位置速度情報を記録するようにした。この結果、予測誤差を抑えることができるので、精度良く軌道予測を行うことができる。
【００４８】
次に、軌道予測を行う処理手順について、図５を用いて説明する。図５は、本実施形態に係る軌道予測を行う処理手順のフローチャートである。
【００４９】
（ステップＳ１０１）軌道予測装置１０は、ステップＳ１〜Ｓ８の過去の位置速度情報を記憶させる処理を行う。
ここで、エポック時刻ｔｏｅに近い偶数正時ＴＯＥ２を時刻ｔ（ｋ）（ｋ＝０，・・・，ｎ−１）と表現する。以降の処理上、位置はｒ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での位置のｘ成分がｒ（３ｋ）、位置のｙ成分がｒ（３ｋ＋１）、位置のz成分がｒ（３ｋ＋２）である。速度はｖ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での速度ｘ成分がｖ（３ｋ）、速度ｙ成分がｖ（３ｋ＋１）、速度ｚ成分がｖ（３ｋ＋２）である。
ステップＳ１０１終了後、ステップＳ１０２に進む。
【００５０】
（ステップＳ１０２）演算開始判定部１４５は、衛星位置・速度記憶部１４０から過去の位置速度情報、衛星型番号を読み出し、記憶部１１５から電離層遅延補正係数、ＧＰＳ時刻、及びＵＴＣ関係情報を読み出す。演算開始判定部１４５は、読み出した情報に基づき、演算開始条件を満たしているか否かを判定する。演算開始条件とは、以下の５つである。
１．過去の位置速度情報が９個以上記憶されていること
２．記憶部１１５に記憶されている衛星型番号が有効で、衛星型がＢＬＯＣＫ ＩＩＡ、ＢＬＯＣＫ ＩＩＲ、ＢＬＯＣＫ ＩＩＲ−Ｍであること
３．軌道予測が未完了であること
４．ＲＴＣ（実時間のクロック）の補正が完了していること
５．アルマナック・データが有効であること
演算開始条件を満たしていないと判別した場合、ステップＳ１０１に戻る。
演算開始条件を満たしていると判別した場合、演算開始判定部１４５は、パラメータ推定部１５０に読み出した情報と判定結果を出力する。
ステップＳ１０２終了後、ステップＳ１０３に進む。
【００５１】
（ステップＳ１０３）パラメータ推定部１５０は、演算開始判定部１４５が出力する読み出された情報と判定結果を用いて、以下の１１個の成分を有するベクトルである推定パラメータｓを推定する。
成分１ ｓ（０）；時刻ｔ（０）での衛星位置のｘ成分（初期値ｒ（０））
成分２ ｓ（１）；時刻ｔ（０）での衛星位置のｙ成分（初期値ｒ（１））
成分３ ｓ（２）；時刻ｔ（０）での衛星位置のｚ成分（初期値ｒ（２））
成分４ ｓ（３）；時刻ｔ（０）での衛星速度のｘ成分（初期値ｖ（０））
成分５ ｓ（４）；時刻ｔ（０）での衛星速度のｙ成分（初期値ｖ（１））
成分６ ｓ（５）；時刻ｔ（０）での衛星速度のｚ成分（初期値ｖ（２））
成分７ ｓ（６）；太陽輻射圧パラメータのワイバイアス（Ｙ−ｂｉａｓ）（初期値０）
成分８ ｓ（７）；太陽輻射圧パラメータのスケールファクタ（ｓｃａｌｅ Ｆａｃｔｏｒ）（初期値１）
成分９ ｓ（８）；地球自転パラメータの極運動ｘ方向（初期値０）
成分１０ ｓ（９）；地球自転パラメータの極運動ｙ方向（初期値０）
成分１１ ｓ（１０）；地球自転パラメータの地球自転速度変化（初期値０）
パラメータ推定部１５０は、推定した推定パラメータを軌道予測演算部１５５に出力する。ステップＳ１０３終了後、ステップＳ１０４に進む。なお、ワイバイアスとは、軌道におけるＧＰＳ衛星２の太陽パネル（ソーラーパネル）の梁方向（Ｂｅａｍ Ｄｉｒｅｃｔｉｏｎｓ）であるｙ軸に対するバイアス成分であり、スケールファクタとは、スケール係数である。太陽輻射圧、ワイバイアス及びスケースファクタについては、例えば、非特許文献１（” New Empirically Derived Solar Radiation Pressure Model for Global Positioning System Satellites”, Y.Bar-Sever and D.Kuang, IPN Progress Report 42-159, Nov 15, 2004）参照。
【００５２】
（ステップＳ１０４）軌道予測演算部１５５は、パラメータ推定部１５０が出力する推定パラメータｓを使用し、後述するように、例えば将来３日間２時間ごとの時刻ｔe（ｋ）（ｋ＝０，・・・，３６）での衛星位置ｐe（ｉ）（ｉ＝０，・・・，１１０）を軌道計算する。ステップＳ１０４終了後、ステップＳ１０５に進む。
【００５３】
（ステップＳ１０５）軌道予測演算部１５５は、算出した予測軌道値を予測結果記憶部１６０に記憶させる。ステップＳ１０５終了後、ステップＳ１０６に進む。
【００５４】
（ステップＳ１０６）測位算出部２０は、予測結果記憶部１６０から予測軌道値を読み出し、読み出した予測軌道値を外挿し衛星位置を算出し、公知の手法で測位を行う。
【００５５】
次に、ステップＳ１０３〜Ｓ１０４で行う衛星位置ｐ０（ｉ）の軌道計算の方法について、図６を用いて説明する。図６は、本実施形態に係る衛星位置ｐ０（ｉ）の軌道計算の処理手順のフローチャートである。
（ステップＳ２０１）まず、パラメータ推定部１５０は、演算開始判定部１４５が出力する読み出された情報と判定結果を用いて、推定パラメータｓで軌道計算する（ステップＳ１０３）。
時刻ｔ（０）での軌道における衛星位置において、半径方向、進行方向、軌道面垂直方向の各単位ベクトルｄ１、ｄ２、ｄ３は、次式（２）〜（４）のように表される。
【００５６】
【数２】
【００５７】
【数３】
【００５８】
【数４】
【００５９】
なお、式（２）において、ｒ（０）、ｒ（１）、ｒ（２）はＧＰＳ衛星２の位置を表し、√（（ｒ（０））２＋（ｒ（１））２＋（ｒ（２））２）は地球からＧＰＳ衛星２との距離を表している。
【００６０】
（ステップＳ２０２）パラメータ推定部１５０は、最小二乗法による推定を行うために、推定パラメータの値をΔｓずらす。パラメータ推定部１５０は、ｊ＝０に対して推定パラメータのｓ（０），ｓ（１），ｓ（２）をｓ（０）＋ｄ１（０）Δｓ（０），ｓ（１）＋ｄ１（１）Δｓ（０），ｓ（２）＋ｄ１（２）Δｓ（０）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。
【００６１】
Ａ（ｉ，０）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（０） ・・・（５）
【００６２】
式（５）において、ｐ０（ｉ）は、Δｓ（０）ずらす前のＧＰＳ衛星２の位置であり、ｐ（ｉ）は、Δｓ（０）ずらした後のＧＰＳ衛星２の位置である。このため、式（５）において、行列Ａは、各パラメータをΔｓ（０）ずらした時に、ＧＰＳ衛星２の位置がどれだけずれたかの変化量である偏微分値の行列を表している。また、ステップＳ２０２の計算は、ＧＰＳ衛星２の位置を半径方向ｄ１にΔｓ（０）ずらした衛星位置の計算を行っている。また、位置をずらす量Δｓ（０）は、予め定められている量を用いる。そして、Δｓ（０）は、例えば、実験により算出する。この場合、パラメータ推定部１５０での桁落ちを考慮し、ｓの値に応じて、Δｓｊを設定するようにしてもよい。
【００６３】
（ステップＳ２０３）軌道予測演算部１５５は、ｊ＝１に対して推定パラメータのｓ（０），ｓ（１），ｓ（２）をｓ（０）＋ｄ３（０）Δｓ（１），ｓ（１）＋ｄ３（１）Δｓ（１），ｓ（２）＋ｄ３（２）Δｓ（１）に代え、時刻ｔ（ｋ）でのＧＰＳ衛星２の位置ｐ（ｉ）を軌道計算する。また、ステップＳ２０３の計算は、ＧＰＳ衛星２の位置を垂直方向ｄ３にΔｓ（１）ずらしたＧＰＳ衛星２の位置の計算を行っている。
【００６４】
Ａ（ｉ，１）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（１） ・・・（６）
【００６５】
（ステップＳ２０４）パラメータ推定部１５０は、ｊ＝２に対して推定パラメータのｓ（３），ｓ（４），ｓ（５）をｓ（３）＋ｄ１（０）Δｓ（２），ｓ（４）＋ｄ１（１）Δｓ（２），ｓ（５）＋ｄ１（２）Δｓ（２）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。また、ステップＳ２０４の計算は、ＧＰＳ衛星２の位置を半径方向ｄ１にΔｓ（２）ずらした位置の計算を行っている。
【００６６】
Ａ（ｉ，２）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（２） ・・・（７）
【００６７】
（ステップＳ２０５）パラメータ推定部１５０は、ｊ＝３に対して推定パラメータのｓ（３），ｓ（４），ｓ（５）をｓ（３）＋ｄ２（０）Δｓ（３），ｓ（４）＋ｄ２（１）Δｓ（３），ｓ（５）＋ｄ２（２）Δｓ（３）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。また、ステップＳ２０５の計算は、ＧＰＳ衛星２の位置を進行方向ｄ２にΔｓ（３）ずらした位置の計算を行っている。
【００６８】
Ａ（ｉ，３）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（３） ・・・（８）
【００６９】
すなわち、本実施形態では、ＧＰＳ衛星２の位置の半径方向と垂直方向、ＧＰＳ衛星２の速度の半径方向と進行方向の計４つのパラメータを推定する。なお、ステップＳ２０２とＳ２０３において、半径方向と垂直方向にのみずらし、進行方向にずらした場合の位置の計算を行わない理由は、進行方向にずらしても、将来のＧＰＳ衛星２の位置への影響が少ないためである。すなわち、地球からの重力の位置への変化量の影響が少ないため、進行方向は推定しない。同様に、ステップＳ２０４とＳ２０５において、半径方向と進行方向にのみずらし、垂直方向にずらした場合の位置の計算を行わない理由は、垂直方向にずらしてもＧＰＳ衛星２の速度への影響が少ないためである。すなわち、地球からの重力の加速度への変化量の影響が少ないため、進行方向は推定しない。
【００７０】
（ステップＳ２０６）パラメータ推定部１５０は、ｊ＝４，・・・，８（成分７〜１１の太陽輻射圧と地球自転との５つパラメータ）に対して推定パラメータの１つの成分ｓ（ｊ＋２）をｓ（ｊ＋２）＋Δｓ（ｊ）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。
【００７１】
Ａ（ｉ，ｊ）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（ｊ） ・・・（９）
【００７２】
ただし、式（５）〜式（７）において、ｉ＝０，・・・，３ｎ−１であり、式（９）において、ｊ＝４，・・・，８である。
【００７３】
（ステップＳ２０７）パラメータ推定部１５０は、次式（１０）を用いて、最小二乗法により推定パラメータｓを算出する。
【００７４】
ｓ＝ｓ＋Ｂ（Ａ×ＡＴ）−１×ＡＴ（ｒ−ｐ０） ・・・（１０）
【００７５】
式（１０）において、Ｔは転置行列を表し、ｒはエフェメリス・データから求めたＧＰＳ衛星２の位置を表し、ｐ０は、推定パラメータから軌道計算したＧＰＳ衛星２の位置を表し、「−１」は逆行列を表している。推定したＧＰＳ衛星２のずれは、推定パラメータのずれによるものである。このため、パラメータ推定部１５０は、推定パラメータのずれ量を式（５）〜（９）で求めて更新する。
また、式（１０）において、Ｔは転置行列を表し、行列Ｂは、次式（１１）である。なお、行列Ａは、ＧＰＳ衛星２の位置の半径方向と垂直方向の成分、ＧＰＳ衛星２の速度の半径方向と進行方向の成分である。行列Ｂは、この行列Ａの成分を、元のｘｙｚ方向の成分に座標系を変換する役割をしている。
【００７６】
【数５】
【００７７】
（ステップＳ２０８）パラメータ推定部１５０は、推定パラメータｓが収束したか否かを判別する。例えば、ステップＳ２０１〜Ｓ２０８を繰り返すことで、Δｓ変化させても、予め定められている範囲内のずれしか発生しなくなった場合、収束したとする。
推定パラメータｓが収束したと判別した場合（ステップＳ２０８；Ｙｅｓ）、ステップＳ２０９に進む。推定パラメータの１つの成分ｓが収束していないと判別した場合（ステップＳ２０８；Ｎｏ）、ステップＳ２０１に戻る。
（ステップＳ２０９）軌道予測演算部１５５は、この推定パラメータの１つの成分sを使用し、例えば将来３日間２時間ごとの時刻ｔe（ｋ）（ｋ＝０，・・・，３６）での衛星位置ｐe（ｉ）（ｉ＝０，・・・，１１０）を軌道計算する。
【００７８】
次に、パラメータ推定部１５０および軌道予測演算部１５５が行う軌道計算の例について説明する。
（手順１）推定パラメータｓの最初の６個のパラメータ衛星位置速度をＥＣＥＦ座標系
からＥＣＩ（地心慣性；Earth Centered Inertial）座標系に変換する。
（手順２）ＧＰＳ時刻をＵＴＣ時刻に変換する。
（手順３）ＧＰＳ衛星２に加わる力（地球からの重力、月からの重力、太陽からの重力、太陽輻射圧）を計算する。
（手順４）ＧＰＳ衛星２に加わる力を加速度に変換し数値積分し、将来のＧＰＳ衛星２の位置と速度を算出する。用いる数値積分法は、例えば、ＲＫＦ４５及びＣｏｗｅｌｌである。座標系、時刻系、衛星に加わる太陽輻射圧以外の力、数値積分については、例えば、非特許文献２（‘’Satellite Orbits, Models, Methods and Application’’, O.Montenbruck and E.Gill, Springer,2005）の２章「Ｉｎｔｒｏｄｕｃｔｏｒｙ Ａｓｔｒｏｄｙｎａｍｉｃｓ」、３章「Ｆｏｒｃｅ Ｍｏｄｅｌ」、４章「Ｎｕｍｅｒｉｃａｌ Ｉｎｔｅｇｒａｔｉｏｎ」、５章「Ｔｉｍｅ ａｎｄ Ｒｅｆｅｒｅｎｃｅ Ｓｙｓｔｅｍｓ」、及び８章の「Ｏｒｂｉｔ Ｄｅｔｅｒｍｉｎａｔｉｏｎ ａｎｄ Ｐａｒａｍｅｔｅｒ Ｅｓｔｉｍａｔｉｏｎ」参照。
（手順５）算出したＧＰＳ衛星２の位置と速度をＥＣＩ座標系からＥＣＥＦ座標系に変換する。
【００７９】
このように算出された位置速度情報、及び実測により取得された位置速度情報、ＵＴＣ関係情報、ＧＰＳ時刻、電離層遅延補正係数を用いて、測位算出部２０は、公知の手法で位置測定を行う。
【００８０】
なお、従来行われている軌道予測処理は、例えば以下である。
（手順１０１）本実施形態と同様に、従来の軌道予測装置は、各種の情報をＧＰＳ衛星２から受信する。
（手順１０２）従来の軌道予測速装置は、本実施形態と同様に１１個の成分を有するベクトルである推定パラメータｓを推定する。
（手順１０３）従来の軌道予測装置は、最小二乗法による予測推定を行うために、推定パラメータの値をΔｓずらす。従来の軌道予測装置は、推定パラメータｓの１つの成分ｓ（ｊ）をｓ（ｊ）＋Δｓ（ｊ）に代え、時刻ｔ（ｋ）（ｋ＝０，・・・，ｎ−１）での衛星位置ｐ（ｉ）（ｉ＝０，・・・，３ｎ−１）を軌道計算する。なお、ｉは、ＧＰＳ衛星２の位置の何番目がずれるのかを表している。ｊは、０から１０であり、どのパラメータをずらしたかを表している。
【００８１】
Ａ（ｉ，ｊ）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（ｊ） ・・・（１２）
【００８２】
式（１２）において、ｐ０（ｉ）は、Δｓ（ｊ）ずらす前のＧＰＳ衛星２の位置であり、ｐ（ｉ）は、Δｓ（ｊ）ずらした後のＧＰＳ衛星２の位置である。
（手順１０４）従来の軌道予測装置は、式（１２）で求めた誤差量を用いて、次式（１３）を用いて、推定パラメータｓを算出する。
【００８３】
ｓ＝ｓ＋（Ａ×ＡＴ）−１×ＡＴ（ｒ−ｐ０） ・・・（１３）
【００８４】
式（１３）において、Ｔは転置行列を表し、ｒはエフェメリス・データから求めたＧＰＳ衛星２の位置を表し、ｐ０は、推定パラメータから軌道計算したＧＰＳ衛星２の位置を表し、「−１」は逆行列を表している。推定したＧＰＳ衛星２の位置に対するずれは、推定パラメータのずれによるものであるため、推定パラメータのずれ量を式（１２）で求めて更新する。
従来の軌道予測装置は、推定パラメータｓが収束するまで、手順１０２から１０４を繰り返す。従来の軌道予測装置は、このように算出された推定パラメータｓを用いて将来のＧＰＳ衛星２の位置ｐｅを軌道計算していた。
また、従来技術では、位置と速度について、衛星位置の半径方向と進行方向及び垂直方向、衛星速度の半径方向と進行方向及び垂直方向の計６つのパラメータについて推測していた。
【００８５】
一方、本実施形態では、上述したように、推定パラメータｓの算出において、半径方向と垂直方向に対してΔｓずらして衛星の位置の推定し、半径方向と進行方向に対してΔｓずらして衛星の速度の推定を行うようにした。このように、本実施形態では、１１個の成分の内、９つのパラメータを推定するようにしたので、従来技術の１１個のパラメータを推定する場合と比較して推測にかかる演算コストを下げることができる。すなわち、パラメータ数を減らしたので、最小二乗法を行う行列を小さくできるため、演算に要する記憶部の容量を小さくでき、さらに演算に係る時間を短くできる。このように推定するパラメータを削減しても、従来技術と同等に、精度良く軌道予測を行うことができる。
【００８６】
以上のように、本実施形態の軌道予測装置１０は、エフェメリス・データから算出した位置速度情報を、時間間隔が均等になるように衛星位置・速度記憶部１４０に記憶させる。また、軌道予測装置１０は、衛星位置・速度記憶部１４０に記憶されている位置速度情報が、予め定めた個数以上の場合、時間間隔が均等になるように衛星位置・速度記憶部１４０から位置速度情報を削除する。このように均等に記憶された位置速度情報を用いて、軌道予測を行うので、従来技術と比較して、予測誤差を小さくできる効果がある。
また、推定パラメータの内、ＧＰＳ衛星２の位置の進行方向、ＧＰＳ衛星２の速度の垂直方向の計２つのパラメータについて推定しないようにした。このため、最小二乗法で用いる行列を従来技術より小さくできるので、推定に必要なメモリ容量を削減できる。また、最小二乗法で用いる行列を小さくできるので、演算に係る時間を短くできる効果もある。このため、軌道予測装置１０のコストを下げることができ、さらに演算時間を短くできるため、起動後に従来技術より短時間で測位を開始できる。この結果、利用者にとって利便性が増すので、使用頻度を増す効果もある。
【００８７】
また、本実施形態の軌道予測装置１０は、衛星位置の半径方向と進行方向及び軌道面垂直方向の各単位ベクトルｄ１とｄ２及びｄ３を用いて衛星の位置と速度との各偏微分値である行列Ａを算出し、算出した行列Ａを元のｘｙｚ方向の成分に座標系を変換する行列Ｂを用いて推定パラメータを算出するようにした。このため、本実施形態の軌道予測装置１０は、推定に使用するパラメータを従来技術より削減でき、未知数が少なくなるので、推定に必要な観測方程式も削減できる。この結果、本実施形態の軌道予測装置１０は、従来技術と比較して少ないエフェメリス・データで軌道予測を行えるので、軌道予測に要する時間を削減することができる。この結果、起動後に従来技術より短時間で位置測位を開始できるので、利用者にとって利便性が増すので、使用頻度を増す効果もある。
例えば、一般的なＧＰＳ受信装置では、電源オフ状態から電源オン状態になった後、エフェメリス・データの取得を開始する。そして、エフェメリス・データ取得後、パラメータの推定及び軌道推定を開始する。パラメータの推定にかかる時間が長い場合、推定途中で測位装置の電源をオフ状態にされてしまうと、それまでの演算結果が無効なままである。このため、電源オン状態になった後、再度、エフェメリス・データの取得からやり直す必要が生じる場合があった。
一方、本実施形態の軌道予測装置１０は、軌道推定にかかる時間が短いため、利用者が電源をオフ状態にする前に推定を完了し、推定結果をＧＰＳ受信部内の記憶部に記憶させておくことが可能になる。そして、このように記憶された推定結果を、電源オフ状態からオン状態になった場合に用いて推定を行えるので、本実施形態の軌道予測装置１０は、電源オン状態後に位置測位開始までの時間を短縮できる。
また、従来技術と同数のエフェメリス・データを用いて軌道予測を行った場合、推定に用いるパラメータの数を削減したため、観測方程式の数に対して推定パラメータの数が減ったため、パラメータ間の共分散が減少し，精度の良いパラメータ値を推定できる。この結果、本実施形態の軌道予測装置１０は、推定精度を改善することができる。なお、推定しないパラメータに関しては、実測により、軌道予測に十分な精度のデータであることを確認している。
【００８８】
［第２実施形態］
第１実施形態では、推定パラメータの内、９つのパラメータを推定する例を説明した。本実施形態では、第１実施形態で用いた９個の成分のうち、さらに、太陽輻射圧パラメータのワイバイアスとスケールファクタの２つを推定せず、７個の成分のベクトルによる推定パラメータｓを推定する例を説明する。
なお、太陽輻射圧パラメータのワイバイアスとスケールファクタは、ＧＰＳ衛星２毎に固有の固定値を、予めＧＰＳ衛星２毎に算出して軌道予測装置１０内のいずれかの記憶部か機能部に、予め記憶させておく。例えば、所定の期間のデータを算出して、平均値を記憶させるようにしてもよい。
【００８９】
次に、過去のＧＰＳ衛星２の位置と速度データを記憶させる処理について、図２、図５及び図６を用いて説明する。
ＧＰＳ衛星２の位置と速度情報の記憶と削除は、第１実施形態と同様に行う（図２、ステップＳ１〜Ｓ８）。
第１実施形態と異なるのは、衛星位置ｐ０（ｉ）の軌道計算の方法である。
（ステップＳ２０１’）軌道予測演算部１５５は、推定パラメータｓを使用し、時刻ｔ（ｋ）での衛星位置ｐ０（ｉ）（ｉ、・・・、３ｎ−１）を軌道計算する。
推定パラメータｓは、以下の９個の成分のベクトルである。
成分１ ｓ（０）；時刻ｔ（０）での衛星位置のｘ成分（初期値ｒ（０））
成分２ ｓ（１）；時刻ｔ（０）での衛星位置のｙ成分（初期値ｒ（１））
成分３ ｓ（２）；時刻ｔ（０）での衛星位置のｚ成分（初期値ｒ（２））
成分４ ｓ（３）；時刻ｔ（０）での衛星速度のｘ成分（初期値ｖ（０））
成分５ ｓ（４）；時刻ｔ（０）での衛星速度のｙ成分（初期値ｖ（１））
成分６ ｓ（５）；時刻ｔ（０）での衛星速度のｚ成分（初期値ｖ（２））
成分７ ｓ（６）；地球自転パラメータの極運動ｘ方向（初期値０）
成分８ ｓ（７）；地球自転パラメータの極運動ｙ方向（初期値０）
成分９ ｓ（８）；地球自転パラメータの地球自転速度変化（初期値０）
【００９０】
時刻ｔ（０）での衛星位置において、半径方向、進行方向、軌道面垂直方向の各単位ベクトルｄ１、ｄ２、ｄ３は、第１実施形態と同様に式（２）〜（４）のように表される。
【００９１】
ステップＳ２０２〜ステップＳ２０５は、第１実施形態と同様に行う。
（ステップＳ２０６’）軌道予測演算部１５５は、ｊ＝４，５，６に対して推定パラメータの１つの成分ｓ（ｊ＋２）をｓ（ｊ＋２）＋Δｓ（ｊ）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を式（９）のように軌道計算する。ただし、式（９）おいて、ｊ＝４，５，６である。
【００９２】
（ステップＳ２０７’）軌道予測演算部１５５は、次式（１４）を用いて、推定パラメータｓを算出する。
【００９３】
ｓ＝ｓ＋Ｂ（Ａ×ＡＴ）−１×ＡＴ（ｒ−ｐ０） ・・・（１４）
【００９４】
式（１４）は、右辺の演算を実行して、演算結果を新たなｓとして更新することを表している。なお、式（１４）において、Ｔは転置行列を表し、行列Ｂは、次式（１５）である。また、行列Ｂは、行列Ａの成分を、元のｘｙｚ方向の成分に座標系を変換する役割をしている。
【００９５】
【数６】
【００９６】
（ステップＳ２０８）軌道予測演算部１５５は、推定パラメータｓが収束したか否かを判別する。
推定パラメータｓが収束したと判別した場合（ステップＳ２０８；Ｙｅｓ）、ステップＳ２０９に進む。推定パラメータの１つの成分ｓが収束していないと判別した場合（ステップＳ２０８；Ｎｏ）、ステップＳ２０１’に戻る。
（ステップＳ２０９）軌道予測演算部１５５は、この推定パラメータの１つの成分sを使用し、将来の衛星位置ｐｅを軌道計算する。
軌道予測演算部１５５は、第１実施形態と同様に軌道計算を行う。
【００９７】
以上のように、第２実施形態によれば、第１実施形態と比較して、推定パラメータの成分の内、２つについて固定値を用いたため、パラメータを推定する際に最小二乗法で用いる行列が小さくなる。このように最小二乗法で用いる行列を小さくできた結果、本実施形態の軌道予測装置１０では、演算に用いる必要メモリの容量が小さくて済み、さらに演算時間を短くできる効果がある。演算に必要なメモリ容量を少なくすることができたので、本実施形態の軌道予測装置１０を用いた測位装置１では、コストダウンを実現できる。また、演算時間を短くできることにより、測位装置１の電源をオフ状態からオン状態にしたような場合、従来技術より早く起動後に位置測定を開始できるので、本実施形態の軌道予測装置１０では、利用者にとって利便性が増すため、利用者の使用頻度が向上する効果もある。
また、第２実施形態によれば、第１実施形態と比較して、さらに推定パラメータにおける未知数が少なくなるので、本実施形態の軌道予測装置１０では、必要な観測方程式の数を減らせる効果もある。この結果、少数のエフェメリス・データを用いて、軌道予測を行うことが可能になるので、従来技術と比較して、本実施形態の軌道予測装置１０では、利用者にとって利便性が増すため、利用者の使用頻度が向上する効果もある。
【００９８】
なお、本実施形態では、３日分の位置速度情報を衛星位置・速度記憶部１４０に記憶させる例を説明したが、本実施形態の軌道予測装置１０を用いる測位装置１の使用目的に応じて、記憶させる期間を３日より多くしても３日より少なくしてもよい。また、エフェメリス・データを記憶するタイミングとして２時間間隔の例を説明したが、測位装置１の用途に応じて、エフェメリス取得部１３０は、ＧＰＳ衛星２から送信されたエフェメリス・データを２時間毎に受信し、受信した後に、４時間間隔に記憶させるようにして間引くようにしてもいい。また、衛星位置・速度記憶部１４０に記憶させる位置速度情報が１０個の例を説明したが、測位装置１の用途に応じて、１０個以上であってもよい。
【符号の説明】
【００９９】
１・・・測位装置、２・・・ＧＰＳ衛星、１０・・・軌道予測装置、２０・・・測位算出部、１１０・・・情報取得部、１１５・・・記憶部、１２０・・・アルマナック取得部、１２５・・・衛星型抽出部、１３０・・・エフェメリス取得部、１３５・・・衛星位置速度算出部、１４０・・・衛星位置・速度記憶部、１４５・・・演算開始判定部、１５０・・・パラメータ推定部、１５５・・・軌道予測演算部、１６０・・・予測結果記憶部
【技術分野】
【０００１】
本発明は、軌道予測装置、及び軌道予測方法に関する。
【背景技術】
【０００２】
ＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）は、位置測定、測量などに用いられている。地球上空の６つの軌道を２４機以上のＧＰＳ衛星が周回し、各ＧＰＳ衛星は、アルマナック（Ａｌｍａｎａｃ）データとエフェメリス（Ｅｐｈｅｍｅｒｉｓ）データの二種類の軌道情報を含む航法メッセージを送信している。
カーナビゲーション装置などのＧＰＳ受信機は、この航法メッセージを用いてＧＰＳ衛星からＧＰＳ受信機までの信号到達時間を測定する。そしてＧＰＳ受信機は、信号到達時間に光速を掛けることでＧＰＳ衛星とＧＰＳ受信機間の距離（擬似距離）を求めている。この作業を４機以上の衛星に対して行うことで、受信機の位置を三次元上で算出している。アルマナック・データは、軌道上における全ての衛星に関する軌道情報であり、１日に１回、データが更新される。エフェメリス・データには、各衛星の正確な位置情報と信号を発射した時刻情報が含まれており、２時間毎にデータが更新される。
【０００３】
ＧＰＳ受信機は、受信した各情報を自機の記憶部に、例えば電源を切る前に受信したアルマナック・データとエフェメリス・データを保存している。これらの情報には、各々有効期限がある。ＧＰＳ受信機の電源がオン状態で、ＧＰＳ衛星からの送信情報を送信毎に受信できていない場合など、ＧＰＳ受信機は、再度、送信情報を受信する必要がある。
アルマナック・データの再受信に要する時間は、例えば１２分３０秒であり、エフェメリス・データの再受信に要する時間は、３０秒である。
どちらかのデータの有効期限が切れにより再受信している期間、ＧＰＳ受信機は、位置測定を行えない。このため、エフェメリス・データの有効期限が切れた場合、ＧＰＳ受信機は、以前受信したエフェメリス・データを用いて、現在のＧＰＳ衛星の位置を予測し、現在のエフェメリス・データが利用できない場合に現在のＧＰＳ衛星の状態を使用して位置測定を行うことが提案されている（例えば、特許文献１参照）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】米国特許第７８３９３３０号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、特許文献１に記載の従来技術では、測位を行う際に用いる推定パラメータの個数が多いため、演算にかかるコストが高いという課題があった。
【０００６】
本発明は、上記の問題点に鑑みてなされたものであって、演算にかかるコストを削減できる軌道予測装置、及び軌道予測方法を提供することを目的としている。
【課題を解決するための手段】
【０００７】
上記目的を達成するため、本発明に係る軌道予測装置は、測位衛星からの信号を受信して受信データを取得する情報取得部と、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出部と、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽から前記測位衛星への光の圧力である太陽輻射圧の影響を示す値、及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ部と、前記推定パラメータ部で推定された推定パラメータを用いて前記測位衛星の位置を推定する軌道予測部と、を備えることを特徴としている。
【０００８】
また、本発明の軌道予測装置において、前記推定パラメータ部は、前記測位衛星の半径方向と進行方向及び軌道面垂直方向との成分による座標系に変換する行列を用いて前記推定パラメータを推定するようにしてもよい。
【０００９】
また、本発明の軌道予測装置において、前記推定パラメータ部は、前記測位衛星の位置の進行方向成分に対して予め算出された値を用い、前記測位衛星の速度の軌道面垂直方向成分に対して予め算出された値を用い、前記推定パラメータを推定するようにしてもよい。
【００１０】
また、本発明の軌道予測装置において、前記推定パラメータ部は、前記測位衛星に対する太陽輻射圧のワイバイアスとスケール係数とに対して各々予め算出した値を用い、前記推定パラメータを推定するようにしてもよい。
【００１１】
上記目的を達成するため、本発明は、軌道予測装置における軌道予測方法であって、情報取得部が、測位衛星からの信号を受信して受信データを取得する情報取得過程と、衛星位置速度算出部が、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出過程と、推定パラメータ部が、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽輻射圧の影響を示す値及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ過程と、軌道予測部が、前記推定パラメータ手順で推定された前記推定パラメータを用いて前記測位衛星の位置を推定する軌道予測過程と、を含むことを特徴としている。
【発明の効果】
【００１２】
本発明によれば、軌道における測位衛星の半径方向と進行方向及び軌道面垂直方向の各単位ベクトルと、位置速度情報とを用いて、測位衛星の位置と測位衛星の速度と太陽から測位衛星への光の圧力である太陽輻射圧の影響及び地球自転の影響を含む推定パラメータを算出するようにした。このように測位衛星の半径方向と進行方向及び軌道面垂直方向の各単位ベクトルを用いて推定パラメータを推定するため、推定するパラメータを削減することができる。この結果、最小二乗法で用いる行列を従来技術より小さくできるので、推定に必要なメモリ容量を削減できる。また、最小二乗法で用いる行列を小さくできるので、演算に係る時間を短くできる。
【図面の簡単な説明】
【００１３】
【図１】第１実施形態に係る軌道予測装置１０の概略ブロック図である。
【図２】同実施形態に係る過去のＧＰＳ衛星位置と速度データを記憶させる処理のフローチャートである。
【図３】同実施形態に係る衛星位置・速度記憶部１４０に記憶される情報と、削除される情報の一例を説明する図である。
【図４】保存間隔毎の予測誤差の一例を説明する図である。
【図５】同実施形態に係る軌道予測を行う処理手順のフローチャートである。
【図６】第２実施形態に係る衛星位置ｐ０（ｉ）の軌道計算の処理手順のフローチャートである。
【発明を実施するための形態】
【００１４】
［第１実施形態］
以下、図面を用いて、本発明の実施形態について説明する。
まず、一般的なＧＰＳ衛星の軌道予測の概略原理について説明する。
【００１５】
地球上空の６つの軌道を２４機以上のＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）衛星（以下、衛星という）が周回し、各衛星は、航法メッセージを送信している。航法メッセージは、フレームを単位に送信している。１フレームの送信には、３０秒を要する。１フレームは、５組のサブフレームで構成され、サブフレーム１から５が順番に送信される。また、サブフレーム１から３は、各衛星に固有の情報を含み、毎回同じ内容で構成される。サブフレーム４と５は、全衛星が同じ内容の情報で構成され、サブフレーム毎にページ１から２５で構成されている。このため、サブフレーム４と５の全ての情報を送信するには、２５フレームが必要である。このため、ＧＰＳ受信機が航法メッセージの全情報を得るには、１２分３０秒の時間がかかる。サブフレーム１は、各衛星の状態（正常に動作しているか否か）、衛星が送信する衛星のクロック誤差を補正するための係数であるクロック補正係数などにより構成されている。サブフレーム２は、各衛星の軌道情報（エフェメリス（Ｅｐｈｅｍｅｒｉｓ）・データ）１／２で構成されている。サブフレーム３は、各衛星の軌道情報（エフェメリス・データ）２／２で構成されている。サブフレーム４は、ＧＰＳ受信機が受信する信号が、電離層により遅延する量を補正するための係数である電離層遅延補正係数、ＧＰＳ時刻とＵＴＣ（協定世界時間；Universal Time, Coordinated））との関係を表す情報であるＵＴＣ関係情報、全衛星の軌道情報（アルマナック（Ａｌｍａｎａｃ）・データ）１／２などで構成されている。サブフレーム５は、全衛星の軌道情報（アルマナック・データ）２／２で構成されている。また、各サブフレームの先頭には、ＧＰＳ時刻を示す情報が含まれている。なお、ＧＰＳ時刻とは、１週間を単位として衛星で管理されている時間であり、毎週、日曜日の０時からの経過時間で表される情報である。
【００１６】
また、エフェメリス・データは、衛星の位置を計算するために必要な軌道の６要素（昇交点赤経、軌道傾斜角、近地点引数、軌道長半径、離心率、真近点角）のデータ、各補正値、及び軌道のエポック時刻ｔｏｅ（軌道の元期（ｅｐｈｅｍｅｒｉｓ ｒｅｆｅｒｅｎｃｅ ｔｉｍｅ））などにより構成されている。エフェメリス・データは、２時間毎に更新される。また、エフェメリス・データの有効期間は、２時間±２時間である。
アルマナック・データは、各衛星の軌道情報である。また、アルマナック・データは、約１日に１度の頻度で更新されている。また、アルマナック・データの有効期間は、エポック時刻±３日間である。
【００１７】
ＧＰＳ受信機は、受信したこれらのデータを記憶部に記憶する。記憶されているアルマナック・データ及びエフェメリス・データ共に有効な状態から、ＧＰＳ受信機が測位を開始することをホット・スタートという。また、エフェメリス・データが無効でありアルマナック・データのみ有効な状態から、ＧＰＳ受信機が測位を開始することをウォーム・スタートという。アルマナック・データ及びエフェメリス・データ共に無効な状態から、ＧＰＳ受信機が測位を開始することをコールド・スタートという。
そして、ホット・スタートの場合、ＧＰＳ受信機は、電源投入後に最初の位置情報が出力されるまでの時間は、記憶されているエフェメリス・データを用いて測位できるので数秒程度である。一方、ウォーム・スタートの場合、ＧＰＳ受信機は、エフェメリス・データを再取得する必要があるので、取得に３０秒以上の時間がかかり、さらに取得したエフェメリス・データの検査にも時間がかかる。このため、ウォーム・スタート場合、電源投入後に最初の位置情報が出力されるまでの時間は、３０秒以上である。さらに、コールド・スタートの場合、全ての航法メッセージを取得する必要があるので、測位信号の捕獲と航法メッセージの取得に数分〜１２分３０秒かかる。このため、電源投入後に最初の位置情報が出力されるまでの時間は、数分〜１２分３０秒以上である。
【００１８】
軌道予測とは、以前収集されたエフェメリス・データを利用し、未来の衛星位置を予測する技術である。軌道予測装置は、この予測した衛星位置を測位に用いることで、ウォーム・スタート時にエフェメリス・データの受信完了を待たずに測位開始することができるため、ＧＰＳ初期位置算出時間（ＴＴＦＦ）を短縮することができる。ＧＰＳ受信機が衛星からの情報取得後、電源オフ状態の期間が短ければ、以前取得したデータを用いて、衛星の位置を把握でききるため、データの再取得までの時間が短くなる。
そして、軌道予測装置は、エフェメリス・データから計算される衛星位置・速度をバックアップ・メモリに保存する。次に、軌道予測装置は、軌道予測計算に必要なデータが揃い次第、軌道予測パラメータを推定し、推定した軌道予測パラメータを使用して軌道を予測する。次に、軌道予測装置は、予測した軌道をバックアップ・メモリへ保存する。そして、軌道予測装置は、エフェメリス・データが無効の時、バックアップ・メモリに保存された予測した軌道を使用して衛星位置を計算する。
【００１９】
衛星の位置は、衛星の初期位置と初期速度が与えられ、衛星にかかる加速度が計算され、積分が行われることで、計算される。衛星にかかる加速度ベクトルａｓａｔは、地球重力定数をＧＭｅ、地球中心から衛星までの位置ベクトルをｒとし、地球重力の地球非球面による加速度ベクトルをａ１、太陽重力による加速度ベクトルをａ２、月重力による加速度ベクトルをａ３、太陽輻射圧による加速度ベクトルをａ４とすると、次式（１）のように表される。なお、太陽輻射圧とは、太陽から飛んでくる光子（Ｐｈｏｔｏｎ）がＧＰＳ衛星２の表面に衝突して発生する力である。
【００２０】
【数１】
【００２１】
一般的な軌道予測プログラムでは、式（１）の加速度ａｓａｔを数値積分により時間積分することで、ＧＰＳ衛星の位置と速度を算出する。用いる数値積分法は、例えば、ＲＫＦ４５（ルンゲクッタフェールベルグ法；Runge-Kutta-Fehlberg method）及びＣｏｗｅｌｌである。
【００２２】
次に、本実施形態に係る軌道予測装置の構成について、図１を用いて説明する。図１は、本実施形態に係るＧＰＳの概略ブロック図である。
図１に示すように、ＧＰＳは、測位装置１、ＧＰＳ衛星２から構成されている。測位装置（ＧＰＳ受信機）１は、軌道予測装置１０、測位算出部２０を備えている。なお、ＧＰＳ衛星２は、１機ではなく、２４機以上の衛星から構成される。
軌道予測装置１０は、情報取得部１１０、記憶部１１５、アルマナック取得部１２０、衛星型抽出部１２５、エフェメリス取得部１３０、衛星位置速度算出部１３５、衛星位置・速度記憶部１４０、演算開始判定部１４５、パラメータ推定部１５０、軌道予測演算部１５５、予測結果記憶部１６０を備えている。
【００２３】
情報取得部１１０は、ＧＰＳ衛星２が送信する航法メッセージから、ＵＴＣ関係情報、ＧＰＳ時刻、電離層遅延補正係数を取得し、取得したＧＰＳ時刻、ＵＴＣ関係情報と電離層遅延補正係数を記憶部１１５に記憶させる。
記憶部１１５には、ＧＰＳ時刻、ＵＴＣ関係情報と電離層遅延補正係数が記憶されている。
【００２４】
アルマナック取得部１２０（情報取得部）は、ＧＰＳ衛星２が送信したアルマナック・データを取得し、取得したアルマナック・データを衛星型抽出部１２５に出力する。
衛星型抽出部１２５は、アルマナック取得部１２０が出力するアルマナック・データから衛星型番号を抽出し、抽出した衛星型番号を衛星位置・速度記憶部１４０に記憶させる。なお、衛星型番号とは、Ｂｌｏｃｋ ＩＩＡ，Ｂｌｏｃｋ ＩＩＲ，Ｂｌｏｃｋ ＩＩＲ−Ｍ，Ｂｌｏｃｋ ＩＩＦなどのＧＰＳ衛星２の型を示す番号である。
【００２５】
エフェメリス取得部１３０（情報取得部）は、ＧＰＳ衛星２が送信するエフェメリス・データを取得し、取得したエフェメリス・データを衛星位置速度算出部１３５に出力する。
【００２６】
衛星位置速度算出部１３５は、エフェメリス取得部１３０が出力したエフェメリス・データから衛星の状態（Ｈｅａｌｔｈ）を抽出する。衛星位置速度算出部１３５は、抽出した衛星の状態が健康であるか否かを判別する。衛星の状態とは、衛星の状態を表すコードで、０以外は何らかの異常が衛星にあることを示している。
衛星の状態が健康であると判別した場合、衛星位置速度算出部１３５は、取得したエフェメリス・データをエポック時刻ｔｏｅに近い偶数正時ＴＯＥ２に変換する。衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に記憶されているＧＰＳ衛星２の位置と速度情報（以下、位置速度情報という）の内、算出した時刻ＴＯＥ２−３日より以前の情報を保存領域から削除する。時刻ＴＯＥ２−３日より以前の情報を保存領域から削除する理由は、エフェメリス・データを均等に３日分記憶させるためである。後述するように、このように均等に３日分、記憶させたエフェメリス・データを用いて軌道予測装置１０が軌道予測を行った場合、予測誤差が少なくなるためである。
衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されているか否かを判別する。衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていると判別した場合、後述するように、衛星位置速度算出部１３５は、ＴＯＥ２の位置速度情報を削除後、衛星位置・速度記憶部１４０の保存領域の空きが、２個分の位置速度情報を記憶できる領域になるまで情報を削除する。衛星位置速度算出部１３５は、時刻ＴＯＥ２と時刻ＴＯＥ２＋２時間の各時刻の各ＧＰＳ衛星２の各位置と各速度を算出し、算出した２つの時刻の位置と速度を位置速度情報として衛星位置・速度記憶部１４０に記憶させる。
衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていないと判別した場合、後述するように、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０の保存領域の空きが、３個分の位置速度情報を記憶できる領域になるまで情報を削除する。衛星位置速度算出部１３５は、時刻ＴＯＥ２、時刻ＴＯＥ−２時間、及び時刻ＴＯＥ２＋２時間の各時刻のＧＰＳ衛星２の各位置と各速度を算出し、算出した２つの時刻の位置と速度を位置速度情報として衛星位置・速度記憶部１４０に記憶させる。
【００２７】
衛星位置・速度記憶部１４０には、衛星型番号、各ＧＰＳ衛星２における各時刻の位置速度情報が記憶されている。
演算開始判定部１４５は、衛星位置・速度記憶部１４０から過去の位置速度情報、衛星型番号を読み出し、記憶部１１５からＧＰＳ時刻、ＵＴＣ関係情報及び電離層遅延補正係数を読み出す。演算開始判定部１４５は、読み出した情報に基づき、後述する演算開始条件を満たしているか否かを判定する。演算開始条件を満たしていると判別した場合、演算開始判定部１４５は、パラメータ推定部１５０に読み出した情報と判定結果を出力する。
【００２８】
パラメータ推定部１５０は、演算開始判定部１４５が出力する読み出された情報と判定結果を用いて、推定パラメータｓを推定する。パラメータ推定部１５０は、推定した推定パラメータを軌道予測演算部１５５に出力する。なお、推定パラメータｓは、後述するように１１個の成分を有するベクトルである。
軌道予測演算部（軌道予測部）１５５は、パラメータ推定部１５０が出力する推定パラメータｓを使用し、後述するように、例えば将来３日間２時間ごとの時刻ｔe（ｋ）（ｋ＝０，・・・，３６）での衛星位置ｐe（ｉ）（ｉ＝０，・・・，１１０）を軌道計算する。軌道予測演算部１５５は、算出した予測軌道値を予測結果記憶部１６０に記憶させる。
予測結果記憶部１６０には、予測軌道値が記憶されている。
測位算出部２０は、予測結果記憶部１６０に記憶されている予測軌道値を読み出し、読み出した予測軌道値を外挿して衛星位置を算出し、公知の手法で位置測位を行う。
【００２９】
なお、軌道予測装置１０は、航法メッセージを受信し、受信した航法メッセージを復調して、復調した航法メッセージから各データを抽出する。情報取得部１１０とアルマナック取得部１２０及びエフェメリス取得部１３０は、この抽出された各データを取得している。例えば、情報取得部１１０が航法メッセージを復調し、復調した航法メッセージから各データを抽出し、抽出した各情報を衛星型抽出部１２５及び衛星位置速度算出部１３５に出力し、記憶部１１５に書き込むようにしてもよい。また、記憶部１１５、衛星位置・速度記憶部１４０、及び予測結果記憶部１６０に分ける例を説明したが、これらの記憶部をまとめて１つの記憶部に記憶させるようにしてもよい。
【００３０】
軌道予測装置１０が行う処理について説明する。
まず、軌道予測装置１０は、航法データの取得を開始する。
次に、アルマナック取得部１２０は、ＧＰＳ衛星２が送信するアルマナック・データを取得し、取得したアルマナック・データを衛星型抽出部１２５に出力する。
次に、衛星型抽出部１２５は、アルマナック取得部１２０が出力するアルマナック・データから衛星型番号を抽出し、抽出した衛星型番号を衛星位置・速度記憶部１４０に記憶させる。
次に、情報取得部１１０は、ＧＰＳ衛星２が送信する航法メッセージのサブフレーム４のデータから、電離層遅延補正係数、ＧＰＳ時刻、及びＵＴＣ関係情報を取得し、取得した電離層遅延補正係数、ＧＰＳ時刻、及びＵＴＣ関係情報を記憶部１１５に記憶させる。
【００３１】
次に、過去のＧＰＳ衛星の位置と速度データを記憶させる処理について、図２を用いて説明する。図２は、本実施形態に係る過去のＧＰＳ衛星の位置と速度データを記憶させる処理のフローチャートである。
【００３２】
（ステップＳ１）エフェメリス取得部１３０は、ＧＰＳ衛星２からエフェメリス・データが送信されるのを待つ。ＧＰＳ衛星２からエフェメリス・データが送信された場合、エフェメリス取得部１３０は、ＧＰＳ衛星２が送信するエフェメリス・データを取得し、取得したエフェメリス・データを衛星位置速度算出部１３５に出力する。ステップＳ１終了後、ステップＳ２に進む。
【００３３】
（ステップＳ２）衛星位置速度算出部１３５は、エフェメリス取得部１３０が出力するエフェメリス・データから各衛星の状態を抽出する。衛星位置速度算出部１３５は、抽出した各衛星の状態が異常を示す値ではないか否か且つエポック時刻ｔｏｅが更新されたデータであるか否かを判別する。
各衛星の状態が異常を示す値ではない、すなわち健康であり且つエポック時刻ｔｏｅが更新されたデータであると判別した場合（ステップＳ２；Ｙｅｓ）、ステップＳ３に進む。各衛星の状態が異常を示す値である、または、エポック時刻ｔｏｅが更新されたデータではないと判別した場合（ステップＳ２；Ｎｏ）、ステップＳ１に戻る。
【００３４】
（ステップＳ３）各衛星の状態が異常を示す値ではない且つエポック時刻ｔｏｅが更新されたデータであると判別した場合、衛星位置速度算出部１３５は、取得したエフェメリス・データをエポック時刻ｔｏｅに近い偶数正時ＴＯＥ２に変換する。なお、取得したエフェメリス・データをエポック時刻ｔｏｅに近い偶数正時ＴＯＥ２に変換する理由は、ＧＰＳ衛星２から送信されたエポック時刻ｔｏｅが、２時間間隔ではなく、まれにずれることがあるためである。
時刻ＴＯＥ２を時刻ｔ（ｋ）（ｋ＝０，・・・，ｎ−１）と表現する。以降の処理上、位置はｒ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での位置のｘ成分がｒ（３ｋ）、位置のｙ成分がｒ（３ｋ＋１）、位置のz成分がｒ（３ｋ＋２）である。速度はｖ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での速度ｘ成分がｖ（３ｋ）、速度ｙ成分がｖ（３ｋ＋１）、速度ｚ成分がｖ（３ｋ＋２）である。
衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に記憶されている位置速度情報の内、算出した時刻ＴＯＥ２−３日より以前の情報を保存領域から削除する。ステップＳ３終了後、ステップＳ４に進む。
【００３５】
（ステップＳ４）衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されているか否かを判別する。
衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていると判別した場合（ステップＳ４；Ｙｅｓ）、ステップＳ５に進む。衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていないと判別した場合（ステップＳ４；Ｎｏ）、ステップＳ７に進む。なお、衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていない状態は、例えば、ＧＰＳ衛星２が測位装置１から見て遮蔽されていて、エフェメリス・データを受信できない場合、測位装置１の電源がオフ状態などである。
【００３６】
（ステップＳ５）衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていると判別した場合、衛星位置速度算出部１３５は、ＴＯＥ２の位置速度情報を削除後、衛星位置・速度記憶部１４０の保存領域の空きが、２個分の位置速度情報を記憶できる領域になるまで位置速度情報を削除する。衛星位置速度算出部１３５は、経過時間の古いものから、２ｎ時間（ｎは１以上の自然数）毎の隣り合う情報の内、新しい時刻の情報を削除する。なお、経過時間とは、衛星位置・速度記憶部１４０に１つも記憶されていない状態から位置速度情報の記憶を開始してからの時間の経過であり、図３における縦軸に示した時間である。
具体的には、２時間毎の隣り合う情報の内、新しい時刻の情報を削除する。２時間毎の隣り合う情報を比較して、情報を削除しても、空き領域を確保できない場合は、４時間毎に隣り合う情報の内、新しい方の情報を削除する。空き領域を確保できない場合、衛星位置速度算出部１３５は、６時間毎に隣り合う情報について比較、削除を行う。空き領域を確保できない場合、衛星位置速度算出部１３５は、８時間毎に隣り合う情報について比較、削除を行う。ステップＳ５終了後、ステップＳ６に進む。
【００３７】
（ステップＳ６）次に、衛星位置速度算出部１３５は、時刻ＴＯＥ２と時刻ＴＯＥ２＋２時間の各時刻のＧＰＳ衛星２の各位置と各速度を算出する。次に、衛星位置速度算出部１３５は、算出した２つの時刻の各位置速度情報を衛星位置・速度記憶部１４０に記憶させる。なお、記憶されるＧＰＳ衛星２の位置と速度はＥＣＥＦ（地球中心・地球固定）座標系である。ステップＳ６終了後、ステップＳ１に戻る。
【００３８】
（ステップＳ７）衛星位置・速度記憶部１４０に時刻ＴＯＥ２の位置速度情報が記憶されていないと判別した場合、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０の保存領域の空きが、３個分の位置速度情報を記憶できる領域になるまで位置速度情報を削除する。衛星位置速度算出部１３５は、経過時間の古いものから、２ｎ（ｎは１以上の自然数）時間毎の隣り合う情報の内、新しい時刻の情報を削除する。
具体的には、２時間毎の隣り合う情報を比較して、情報を削除しても、空き領域を確保できない場合は、４時間毎に隣り合う情報の内、新しい方の情報を削除する。空き領域を確保できない場合、衛星位置速度算出部１３５は、６時間毎に隣り合う情報について比較、削除を行う。空き領域を確保できない場合、衛星位置速度算出部１３５は、８時間毎に隣り合う情報について比較、削除を行う。ステップＳ７終了後、ステップＳ８に進む。
【００３９】
（ステップＳ８）次に、衛星位置速度算出部１３５は、時刻ＴＯＥ２、時刻ＴＯＥ２−２時間、及び時刻ＴＯＥ２＋２時間の各時刻のＧＰＳ衛星２の各位置と各速度を算出する。次に、衛星位置速度算出部１３５は、算出した３つの時刻の各位置速度情報を衛星位置・速度記憶部１４０に記憶させる。ステップＳ８終了後、ステップＳ１に戻る。
【００４０】
次に、衛星位置・速度記憶部１４０に記憶される位置速度情報の順番と、削除される順番の具体的な例について、図３を用いて説明する。図３は、本実施形態に係る衛星位置・速度記憶部１４０に記憶される位置速度情報と、削除される位置速度情報の一例を説明する図である。図３において、横軸は保存回数、縦軸は経過時間である。なお、経過時間において、２時間毎のデータであり、「２」は４時間を表している。また、保存回数とは、エフェメリス・データから、ＧＰＳ衛星２の位置と速度を計算し、衛星位置・速度記憶部１４０へ位置速度情報を保存した回数である。経過時間は、保存開始の時刻から各保存された情報の時刻までの経過時間である。
【００４１】
まず、保存回数１では、まだ時刻ＴＯＥ２の位置速度情報が衛星位置・速度記憶部１４０に記憶されていない。このため、図３に示すように、衛星位置速度算出部１３５は、経過時間０、１（２時間後）、及び２（４時間後）の位置速度情報を衛星位置・速度記憶部１４０に記憶する。この場合、衛星位置速度算出部１３５は、ステップＳ８で説明したように、時刻ＴＯＥ２（ｇ２）と、時刻ＴＯＥ２−２時間（ｇ１）、及び時刻ＴＯＥ２＋２時間（ｇ３）の３つの位置速度情報を衛星位置・速度記憶部１４０に記憶する。
保存回数２において、ステップＳ６で説明したように、衛星位置速度算出部１３５は、２つの位置速度情報が衛星位置・速度記憶部１４０に記憶する。この場合、衛星位置速度算出部１３５は、この時点で最新の位置速度情報である位置速度情報ｇ３を削除した後、時刻ＴＯＥ２の位置速度情報を位置速度情報ｇ３’として衛星位置・速度記憶部１４０に記憶させる。位置速度情報ｇ３を削除する理由は、保存回数２における時刻ＴＯＥ２時間は、保存回数１における時刻ＴＯＥ２＋２時間であるためである。そして、衛星位置速度算出部１３５は、時刻ＴＯＥ２＋２時間の位置速度情報を位置速度情報ｇ４として記憶させる。
同様に、保存回数３において、衛星位置速度算出部１３５は、位置速度情報ｇ４を削除した後、時刻ＴＯＥ２の位置速度情報を位置速度情報ｇ４’として記憶させる。そして、時刻ＴＯＥ２＋２時間の位置速度情報を、位置速度情報ｇ５として記憶する。
【００４２】
経過時間４から１０の期間、記憶されている情報がないのは、エフェメリス・データを軌道予測装置１０が取得できなかったことを表している。
保存回数４において、保存回数４における時刻ＴＯＥ２から直近の偶数正時、すなわち経過時間１０（２０時間）にＧＰＳ衛星２の位置と速度情報が衛星位置・速度記憶部１４０に記憶されていないため、ステップＳ８のように、衛星位置速度算出部１３５は、再び３つの時刻の位置速度情報を衛星位置・速度記憶部１４０に記憶させる。このような動作を繰り返していくと、保存回数６で、衛星位置・速度記憶部１４０に保存される位置速度情報が１０個になる。
【００４３】
保存回数６に達した後の衛星位置・速度記憶部１４０に記憶されている位置速度情報を、便宜的にｇ１１〜ｇ２０とする。
保存回数７において、ステップＳ５で説明したように、衛星位置速度算出部１３５は、２つの位置速度情報分の記憶領域を確保する。このため、まず、保存回数６までと同様に、衛星位置速度算出部１３５は、保存回数６で保存された位置速度情報ｇ２０を削除する。次に、衛星位置速度算出部１３５は、位置速度情報を２つ書き込める領域が確保できるまで、衛星位置・速度記憶部１４０に保存されている位置速度情報を削除する。
具体的には、衛星位置速度算出部１３５は、経過時間が一番古い２時間間隔のペアの位置速度情報ｇ１１とｇ１２を比較する。位置速度情報ｇ１１は、保存回数１で記憶された位置速度情報ｇ１であるので、保存回数１における時刻ＴＯＥ２−２時間の位置速度情報である。位置速度情報ｇ１２は、保存回数１で記憶された位置速度情報ｇ１であり、保存回数１における時刻ＴＯＥ２の位置速度情報である。位置速度情報ｇ１２の方が新しい情報であるので、衛星位置速度算出部１３５は、位置速度情報ｇ１２を削除する。このようにして、衛星位置速度算出部１３５は、２つの空き領域を確保する。
次に、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０に、時刻ＴＯＥ２時間の位置速度情報を新たな位置速度情報ｇ２０として記憶させ、時刻ＴＯＥ２＋２時間の位置速度情報をｇ２１として記憶させる。
【００４４】
なお、本実施形態では、ステップＳ５及び保存回数７において、２つの位置速度情報を削除した後に、時刻ＴＯＥの位置速度情報と時刻ＴＯＥ＋２時間の位置速度情報を衛星位置・速度記憶部１４０に記憶させる例を説明した。しかしながら、１つの位置速度情報を削除した後に、時刻ＴＯＥ＋２時間の位置速度情報を記憶させ、時刻ＴＯＥの位置速度情報と上書きするようにしてもよい。具体的には、保存回数７において衛星位置速度算出部１３５は、位置速度情報ｇ１１とｇ１２を比較する。衛星位置速度算出部１３５は、比較結果に基づき、時刻の新しい方の位置速度情報ｇ１２を削除する。
次に、衛星位置速度算出部１４は、時刻ＴＯＥ２時間の位置速度情報を、位置速度情報ｇ２０に上書きし、時刻ＴＯＥ２＋２時間の位置速度情報をｇ２１として記憶させるようにしてもよい。
【００４５】
このように、衛星位置速度算出部１３５は、衛星位置・速度記憶部１４０の位置速度情報に対して削除と追加を繰り返す。
保存回数１３において、衛星位置・速度記憶部１４０に記憶されている位置速度情報を、便宜上ｇ２１〜ｇ３０とする。
保存回数１４において、時刻ＴＯＥ２の位置速度情報が衛星位置・速度記憶部１４０に記憶されていないため、３つの位置速度情報を削除する必要がある。この場合、衛星位置速度算出部１３５は、まず、経過時間が一番古い２時間間隔のペアの位置速度情報ｇ２８とｇ２９を比較し、比較結果に基づき、位置速度情報ｇ２９を削除する。次に、衛星位置速度算出部１３５は、経過時間が一番古い４時間間隔のペアの位置速度情報ｇ２４とｇ２５を比較し、比較結果に基づき、位置速度情報ｇ２５を削除する。次に、衛星位置速度算出部１３５は、経過時間が二番古い４時間間隔のペアの位置速度情報ｇ２６とｇ２７を比較し、比較結果に基づき、位置速度情報ｇ２７を削除する。
以上の処理により、保存回数１３で保存されていた位置速度情報から３個の位置速度情報の削除を完了する。
また、保存回数１５回において、衛星位置速度算出部１３５は、経過時間０の情報が、ステップＳ３で説明した時刻ＴＯＥ２−３日の位置速度情報であるため削除する。
【００４６】
次に、このように、衛星位置・速度記憶部１４０に記憶されている位置速度情報を、間引きながら保存する理由について説明する。図４は、保存間隔毎の予測誤差の一例を説明する図である。図４（ａ）と図４（ｂ）は、位置速度情報を４時間毎、均等に３日間分取得した例である。図４（ｃ）と図４（ｄ）は、位置速度情報を３日間分の内、前後にのみ取取得し、３日に渡って中抜けした例である。図４（ｅ）と図４（ｆ）は、位置速度情報を３日の内、２日間に渡って中抜けした例である。
図４（ａ）、図４（ｃ）及び図４（ｅ）において、横軸は予測誤差の最大値を表し、縦軸はその割合を表している。図４（ｂ）、図（ｄ）及び図（ｆ）において、横軸は時間を表し、白丸は位置速度情報を取得できなかった「中抜け」状態を表し、黒丸は位置速度情報を取得できた状態を表している。
４時間間隔の場合、このように取得された位置速度情報を用いて、ＧＰＳ衛星２の位置を予測すると、図４（ａ）に示すように、予測誤差の最大値は７［ｍ］以下である。中抜けが３日に渡る場合、このように取得された位置速度情報を用いて、ＧＰＳ衛星２の位置を予測すると、図４（ｃ）に示すように、予測誤差の最大値は２０［ｍ］以上もあり、その割合もランダムである。中抜けが２日に渡る場合、このように取得された位置速度情報を用いて、ＧＰＳ衛星２の位置を予測すると、図４（ｅ）に示すように、予測誤差の最大値は２０［ｍ］以上が５０％もある。
このように、予測誤差を抑えるには、中抜けの期間が長い状態は好ましくなく、位置速度情報が均等に保存されていることが好ましい。
【００４７】
このため、本実施形態では、単に１０個の位置速度情報を時系列に記憶させているのではなく、限られた記憶容量内で３日分の位置速度情報を時間間隔が均等に記憶するように、保存した位置速度情報を削除して、新たな位置速度情報を記録するようにした。この結果、予測誤差を抑えることができるので、精度良く軌道予測を行うことができる。
【００４８】
次に、軌道予測を行う処理手順について、図５を用いて説明する。図５は、本実施形態に係る軌道予測を行う処理手順のフローチャートである。
【００４９】
（ステップＳ１０１）軌道予測装置１０は、ステップＳ１〜Ｓ８の過去の位置速度情報を記憶させる処理を行う。
ここで、エポック時刻ｔｏｅに近い偶数正時ＴＯＥ２を時刻ｔ（ｋ）（ｋ＝０，・・・，ｎ−１）と表現する。以降の処理上、位置はｒ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での位置のｘ成分がｒ（３ｋ）、位置のｙ成分がｒ（３ｋ＋１）、位置のz成分がｒ（３ｋ＋２）である。速度はｖ（ｉ）（ｉ＝０，・・・，３ｎ−１）と表現する。時刻ｔ（ｋ）での速度ｘ成分がｖ（３ｋ）、速度ｙ成分がｖ（３ｋ＋１）、速度ｚ成分がｖ（３ｋ＋２）である。
ステップＳ１０１終了後、ステップＳ１０２に進む。
【００５０】
（ステップＳ１０２）演算開始判定部１４５は、衛星位置・速度記憶部１４０から過去の位置速度情報、衛星型番号を読み出し、記憶部１１５から電離層遅延補正係数、ＧＰＳ時刻、及びＵＴＣ関係情報を読み出す。演算開始判定部１４５は、読み出した情報に基づき、演算開始条件を満たしているか否かを判定する。演算開始条件とは、以下の５つである。
１．過去の位置速度情報が９個以上記憶されていること
２．記憶部１１５に記憶されている衛星型番号が有効で、衛星型がＢＬＯＣＫ ＩＩＡ、ＢＬＯＣＫ ＩＩＲ、ＢＬＯＣＫ ＩＩＲ−Ｍであること
３．軌道予測が未完了であること
４．ＲＴＣ（実時間のクロック）の補正が完了していること
５．アルマナック・データが有効であること
演算開始条件を満たしていないと判別した場合、ステップＳ１０１に戻る。
演算開始条件を満たしていると判別した場合、演算開始判定部１４５は、パラメータ推定部１５０に読み出した情報と判定結果を出力する。
ステップＳ１０２終了後、ステップＳ１０３に進む。
【００５１】
（ステップＳ１０３）パラメータ推定部１５０は、演算開始判定部１４５が出力する読み出された情報と判定結果を用いて、以下の１１個の成分を有するベクトルである推定パラメータｓを推定する。
成分１ ｓ（０）；時刻ｔ（０）での衛星位置のｘ成分（初期値ｒ（０））
成分２ ｓ（１）；時刻ｔ（０）での衛星位置のｙ成分（初期値ｒ（１））
成分３ ｓ（２）；時刻ｔ（０）での衛星位置のｚ成分（初期値ｒ（２））
成分４ ｓ（３）；時刻ｔ（０）での衛星速度のｘ成分（初期値ｖ（０））
成分５ ｓ（４）；時刻ｔ（０）での衛星速度のｙ成分（初期値ｖ（１））
成分６ ｓ（５）；時刻ｔ（０）での衛星速度のｚ成分（初期値ｖ（２））
成分７ ｓ（６）；太陽輻射圧パラメータのワイバイアス（Ｙ−ｂｉａｓ）（初期値０）
成分８ ｓ（７）；太陽輻射圧パラメータのスケールファクタ（ｓｃａｌｅ Ｆａｃｔｏｒ）（初期値１）
成分９ ｓ（８）；地球自転パラメータの極運動ｘ方向（初期値０）
成分１０ ｓ（９）；地球自転パラメータの極運動ｙ方向（初期値０）
成分１１ ｓ（１０）；地球自転パラメータの地球自転速度変化（初期値０）
パラメータ推定部１５０は、推定した推定パラメータを軌道予測演算部１５５に出力する。ステップＳ１０３終了後、ステップＳ１０４に進む。なお、ワイバイアスとは、軌道におけるＧＰＳ衛星２の太陽パネル（ソーラーパネル）の梁方向（Ｂｅａｍ Ｄｉｒｅｃｔｉｏｎｓ）であるｙ軸に対するバイアス成分であり、スケールファクタとは、スケール係数である。太陽輻射圧、ワイバイアス及びスケースファクタについては、例えば、非特許文献１（” New Empirically Derived Solar Radiation Pressure Model for Global Positioning System Satellites”, Y.Bar-Sever and D.Kuang, IPN Progress Report 42-159, Nov 15, 2004）参照。
【００５２】
（ステップＳ１０４）軌道予測演算部１５５は、パラメータ推定部１５０が出力する推定パラメータｓを使用し、後述するように、例えば将来３日間２時間ごとの時刻ｔe（ｋ）（ｋ＝０，・・・，３６）での衛星位置ｐe（ｉ）（ｉ＝０，・・・，１１０）を軌道計算する。ステップＳ１０４終了後、ステップＳ１０５に進む。
【００５３】
（ステップＳ１０５）軌道予測演算部１５５は、算出した予測軌道値を予測結果記憶部１６０に記憶させる。ステップＳ１０５終了後、ステップＳ１０６に進む。
【００５４】
（ステップＳ１０６）測位算出部２０は、予測結果記憶部１６０から予測軌道値を読み出し、読み出した予測軌道値を外挿し衛星位置を算出し、公知の手法で測位を行う。
【００５５】
次に、ステップＳ１０３〜Ｓ１０４で行う衛星位置ｐ０（ｉ）の軌道計算の方法について、図６を用いて説明する。図６は、本実施形態に係る衛星位置ｐ０（ｉ）の軌道計算の処理手順のフローチャートである。
（ステップＳ２０１）まず、パラメータ推定部１５０は、演算開始判定部１４５が出力する読み出された情報と判定結果を用いて、推定パラメータｓで軌道計算する（ステップＳ１０３）。
時刻ｔ（０）での軌道における衛星位置において、半径方向、進行方向、軌道面垂直方向の各単位ベクトルｄ１、ｄ２、ｄ３は、次式（２）〜（４）のように表される。
【００５６】
【数２】
【００５７】
【数３】
【００５８】
【数４】
【００５９】
なお、式（２）において、ｒ（０）、ｒ（１）、ｒ（２）はＧＰＳ衛星２の位置を表し、√（（ｒ（０））２＋（ｒ（１））２＋（ｒ（２））２）は地球からＧＰＳ衛星２との距離を表している。
【００６０】
（ステップＳ２０２）パラメータ推定部１５０は、最小二乗法による推定を行うために、推定パラメータの値をΔｓずらす。パラメータ推定部１５０は、ｊ＝０に対して推定パラメータのｓ（０），ｓ（１），ｓ（２）をｓ（０）＋ｄ１（０）Δｓ（０），ｓ（１）＋ｄ１（１）Δｓ（０），ｓ（２）＋ｄ１（２）Δｓ（０）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。
【００６１】
Ａ（ｉ，０）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（０） ・・・（５）
【００６２】
式（５）において、ｐ０（ｉ）は、Δｓ（０）ずらす前のＧＰＳ衛星２の位置であり、ｐ（ｉ）は、Δｓ（０）ずらした後のＧＰＳ衛星２の位置である。このため、式（５）において、行列Ａは、各パラメータをΔｓ（０）ずらした時に、ＧＰＳ衛星２の位置がどれだけずれたかの変化量である偏微分値の行列を表している。また、ステップＳ２０２の計算は、ＧＰＳ衛星２の位置を半径方向ｄ１にΔｓ（０）ずらした衛星位置の計算を行っている。また、位置をずらす量Δｓ（０）は、予め定められている量を用いる。そして、Δｓ（０）は、例えば、実験により算出する。この場合、パラメータ推定部１５０での桁落ちを考慮し、ｓの値に応じて、Δｓｊを設定するようにしてもよい。
【００６３】
（ステップＳ２０３）軌道予測演算部１５５は、ｊ＝１に対して推定パラメータのｓ（０），ｓ（１），ｓ（２）をｓ（０）＋ｄ３（０）Δｓ（１），ｓ（１）＋ｄ３（１）Δｓ（１），ｓ（２）＋ｄ３（２）Δｓ（１）に代え、時刻ｔ（ｋ）でのＧＰＳ衛星２の位置ｐ（ｉ）を軌道計算する。また、ステップＳ２０３の計算は、ＧＰＳ衛星２の位置を垂直方向ｄ３にΔｓ（１）ずらしたＧＰＳ衛星２の位置の計算を行っている。
【００６４】
Ａ（ｉ，１）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（１） ・・・（６）
【００６５】
（ステップＳ２０４）パラメータ推定部１５０は、ｊ＝２に対して推定パラメータのｓ（３），ｓ（４），ｓ（５）をｓ（３）＋ｄ１（０）Δｓ（２），ｓ（４）＋ｄ１（１）Δｓ（２），ｓ（５）＋ｄ１（２）Δｓ（２）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。また、ステップＳ２０４の計算は、ＧＰＳ衛星２の位置を半径方向ｄ１にΔｓ（２）ずらした位置の計算を行っている。
【００６６】
Ａ（ｉ，２）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（２） ・・・（７）
【００６７】
（ステップＳ２０５）パラメータ推定部１５０は、ｊ＝３に対して推定パラメータのｓ（３），ｓ（４），ｓ（５）をｓ（３）＋ｄ２（０）Δｓ（３），ｓ（４）＋ｄ２（１）Δｓ（３），ｓ（５）＋ｄ２（２）Δｓ（３）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。また、ステップＳ２０５の計算は、ＧＰＳ衛星２の位置を進行方向ｄ２にΔｓ（３）ずらした位置の計算を行っている。
【００６８】
Ａ（ｉ，３）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（３） ・・・（８）
【００６９】
すなわち、本実施形態では、ＧＰＳ衛星２の位置の半径方向と垂直方向、ＧＰＳ衛星２の速度の半径方向と進行方向の計４つのパラメータを推定する。なお、ステップＳ２０２とＳ２０３において、半径方向と垂直方向にのみずらし、進行方向にずらした場合の位置の計算を行わない理由は、進行方向にずらしても、将来のＧＰＳ衛星２の位置への影響が少ないためである。すなわち、地球からの重力の位置への変化量の影響が少ないため、進行方向は推定しない。同様に、ステップＳ２０４とＳ２０５において、半径方向と進行方向にのみずらし、垂直方向にずらした場合の位置の計算を行わない理由は、垂直方向にずらしてもＧＰＳ衛星２の速度への影響が少ないためである。すなわち、地球からの重力の加速度への変化量の影響が少ないため、進行方向は推定しない。
【００７０】
（ステップＳ２０６）パラメータ推定部１５０は、ｊ＝４，・・・，８（成分７〜１１の太陽輻射圧と地球自転との５つパラメータ）に対して推定パラメータの１つの成分ｓ（ｊ＋２）をｓ（ｊ＋２）＋Δｓ（ｊ）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を軌道計算する。
【００７１】
Ａ（ｉ，ｊ）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（ｊ） ・・・（９）
【００７２】
ただし、式（５）〜式（７）において、ｉ＝０，・・・，３ｎ−１であり、式（９）において、ｊ＝４，・・・，８である。
【００７３】
（ステップＳ２０７）パラメータ推定部１５０は、次式（１０）を用いて、最小二乗法により推定パラメータｓを算出する。
【００７４】
ｓ＝ｓ＋Ｂ（Ａ×ＡＴ）−１×ＡＴ（ｒ−ｐ０） ・・・（１０）
【００７５】
式（１０）において、Ｔは転置行列を表し、ｒはエフェメリス・データから求めたＧＰＳ衛星２の位置を表し、ｐ０は、推定パラメータから軌道計算したＧＰＳ衛星２の位置を表し、「−１」は逆行列を表している。推定したＧＰＳ衛星２のずれは、推定パラメータのずれによるものである。このため、パラメータ推定部１５０は、推定パラメータのずれ量を式（５）〜（９）で求めて更新する。
また、式（１０）において、Ｔは転置行列を表し、行列Ｂは、次式（１１）である。なお、行列Ａは、ＧＰＳ衛星２の位置の半径方向と垂直方向の成分、ＧＰＳ衛星２の速度の半径方向と進行方向の成分である。行列Ｂは、この行列Ａの成分を、元のｘｙｚ方向の成分に座標系を変換する役割をしている。
【００７６】
【数５】
【００７７】
（ステップＳ２０８）パラメータ推定部１５０は、推定パラメータｓが収束したか否かを判別する。例えば、ステップＳ２０１〜Ｓ２０８を繰り返すことで、Δｓ変化させても、予め定められている範囲内のずれしか発生しなくなった場合、収束したとする。
推定パラメータｓが収束したと判別した場合（ステップＳ２０８；Ｙｅｓ）、ステップＳ２０９に進む。推定パラメータの１つの成分ｓが収束していないと判別した場合（ステップＳ２０８；Ｎｏ）、ステップＳ２０１に戻る。
（ステップＳ２０９）軌道予測演算部１５５は、この推定パラメータの１つの成分sを使用し、例えば将来３日間２時間ごとの時刻ｔe（ｋ）（ｋ＝０，・・・，３６）での衛星位置ｐe（ｉ）（ｉ＝０，・・・，１１０）を軌道計算する。
【００７８】
次に、パラメータ推定部１５０および軌道予測演算部１５５が行う軌道計算の例について説明する。
（手順１）推定パラメータｓの最初の６個のパラメータ衛星位置速度をＥＣＥＦ座標系
からＥＣＩ（地心慣性；Earth Centered Inertial）座標系に変換する。
（手順２）ＧＰＳ時刻をＵＴＣ時刻に変換する。
（手順３）ＧＰＳ衛星２に加わる力（地球からの重力、月からの重力、太陽からの重力、太陽輻射圧）を計算する。
（手順４）ＧＰＳ衛星２に加わる力を加速度に変換し数値積分し、将来のＧＰＳ衛星２の位置と速度を算出する。用いる数値積分法は、例えば、ＲＫＦ４５及びＣｏｗｅｌｌである。座標系、時刻系、衛星に加わる太陽輻射圧以外の力、数値積分については、例えば、非特許文献２（‘’Satellite Orbits, Models, Methods and Application’’, O.Montenbruck and E.Gill, Springer,2005）の２章「Ｉｎｔｒｏｄｕｃｔｏｒｙ Ａｓｔｒｏｄｙｎａｍｉｃｓ」、３章「Ｆｏｒｃｅ Ｍｏｄｅｌ」、４章「Ｎｕｍｅｒｉｃａｌ Ｉｎｔｅｇｒａｔｉｏｎ」、５章「Ｔｉｍｅ ａｎｄ Ｒｅｆｅｒｅｎｃｅ Ｓｙｓｔｅｍｓ」、及び８章の「Ｏｒｂｉｔ Ｄｅｔｅｒｍｉｎａｔｉｏｎ ａｎｄ Ｐａｒａｍｅｔｅｒ Ｅｓｔｉｍａｔｉｏｎ」参照。
（手順５）算出したＧＰＳ衛星２の位置と速度をＥＣＩ座標系からＥＣＥＦ座標系に変換する。
【００７９】
このように算出された位置速度情報、及び実測により取得された位置速度情報、ＵＴＣ関係情報、ＧＰＳ時刻、電離層遅延補正係数を用いて、測位算出部２０は、公知の手法で位置測定を行う。
【００８０】
なお、従来行われている軌道予測処理は、例えば以下である。
（手順１０１）本実施形態と同様に、従来の軌道予測装置は、各種の情報をＧＰＳ衛星２から受信する。
（手順１０２）従来の軌道予測速装置は、本実施形態と同様に１１個の成分を有するベクトルである推定パラメータｓを推定する。
（手順１０３）従来の軌道予測装置は、最小二乗法による予測推定を行うために、推定パラメータの値をΔｓずらす。従来の軌道予測装置は、推定パラメータｓの１つの成分ｓ（ｊ）をｓ（ｊ）＋Δｓ（ｊ）に代え、時刻ｔ（ｋ）（ｋ＝０，・・・，ｎ−１）での衛星位置ｐ（ｉ）（ｉ＝０，・・・，３ｎ−１）を軌道計算する。なお、ｉは、ＧＰＳ衛星２の位置の何番目がずれるのかを表している。ｊは、０から１０であり、どのパラメータをずらしたかを表している。
【００８１】
Ａ（ｉ，ｊ）＝（ｐ（ｉ）−ｐ０（ｉ））／Δｓ（ｊ） ・・・（１２）
【００８２】
式（１２）において、ｐ０（ｉ）は、Δｓ（ｊ）ずらす前のＧＰＳ衛星２の位置であり、ｐ（ｉ）は、Δｓ（ｊ）ずらした後のＧＰＳ衛星２の位置である。
（手順１０４）従来の軌道予測装置は、式（１２）で求めた誤差量を用いて、次式（１３）を用いて、推定パラメータｓを算出する。
【００８３】
ｓ＝ｓ＋（Ａ×ＡＴ）−１×ＡＴ（ｒ−ｐ０） ・・・（１３）
【００８４】
式（１３）において、Ｔは転置行列を表し、ｒはエフェメリス・データから求めたＧＰＳ衛星２の位置を表し、ｐ０は、推定パラメータから軌道計算したＧＰＳ衛星２の位置を表し、「−１」は逆行列を表している。推定したＧＰＳ衛星２の位置に対するずれは、推定パラメータのずれによるものであるため、推定パラメータのずれ量を式（１２）で求めて更新する。
従来の軌道予測装置は、推定パラメータｓが収束するまで、手順１０２から１０４を繰り返す。従来の軌道予測装置は、このように算出された推定パラメータｓを用いて将来のＧＰＳ衛星２の位置ｐｅを軌道計算していた。
また、従来技術では、位置と速度について、衛星位置の半径方向と進行方向及び垂直方向、衛星速度の半径方向と進行方向及び垂直方向の計６つのパラメータについて推測していた。
【００８５】
一方、本実施形態では、上述したように、推定パラメータｓの算出において、半径方向と垂直方向に対してΔｓずらして衛星の位置の推定し、半径方向と進行方向に対してΔｓずらして衛星の速度の推定を行うようにした。このように、本実施形態では、１１個の成分の内、９つのパラメータを推定するようにしたので、従来技術の１１個のパラメータを推定する場合と比較して推測にかかる演算コストを下げることができる。すなわち、パラメータ数を減らしたので、最小二乗法を行う行列を小さくできるため、演算に要する記憶部の容量を小さくでき、さらに演算に係る時間を短くできる。このように推定するパラメータを削減しても、従来技術と同等に、精度良く軌道予測を行うことができる。
【００８６】
以上のように、本実施形態の軌道予測装置１０は、エフェメリス・データから算出した位置速度情報を、時間間隔が均等になるように衛星位置・速度記憶部１４０に記憶させる。また、軌道予測装置１０は、衛星位置・速度記憶部１４０に記憶されている位置速度情報が、予め定めた個数以上の場合、時間間隔が均等になるように衛星位置・速度記憶部１４０から位置速度情報を削除する。このように均等に記憶された位置速度情報を用いて、軌道予測を行うので、従来技術と比較して、予測誤差を小さくできる効果がある。
また、推定パラメータの内、ＧＰＳ衛星２の位置の進行方向、ＧＰＳ衛星２の速度の垂直方向の計２つのパラメータについて推定しないようにした。このため、最小二乗法で用いる行列を従来技術より小さくできるので、推定に必要なメモリ容量を削減できる。また、最小二乗法で用いる行列を小さくできるので、演算に係る時間を短くできる効果もある。このため、軌道予測装置１０のコストを下げることができ、さらに演算時間を短くできるため、起動後に従来技術より短時間で測位を開始できる。この結果、利用者にとって利便性が増すので、使用頻度を増す効果もある。
【００８７】
また、本実施形態の軌道予測装置１０は、衛星位置の半径方向と進行方向及び軌道面垂直方向の各単位ベクトルｄ１とｄ２及びｄ３を用いて衛星の位置と速度との各偏微分値である行列Ａを算出し、算出した行列Ａを元のｘｙｚ方向の成分に座標系を変換する行列Ｂを用いて推定パラメータを算出するようにした。このため、本実施形態の軌道予測装置１０は、推定に使用するパラメータを従来技術より削減でき、未知数が少なくなるので、推定に必要な観測方程式も削減できる。この結果、本実施形態の軌道予測装置１０は、従来技術と比較して少ないエフェメリス・データで軌道予測を行えるので、軌道予測に要する時間を削減することができる。この結果、起動後に従来技術より短時間で位置測位を開始できるので、利用者にとって利便性が増すので、使用頻度を増す効果もある。
例えば、一般的なＧＰＳ受信装置では、電源オフ状態から電源オン状態になった後、エフェメリス・データの取得を開始する。そして、エフェメリス・データ取得後、パラメータの推定及び軌道推定を開始する。パラメータの推定にかかる時間が長い場合、推定途中で測位装置の電源をオフ状態にされてしまうと、それまでの演算結果が無効なままである。このため、電源オン状態になった後、再度、エフェメリス・データの取得からやり直す必要が生じる場合があった。
一方、本実施形態の軌道予測装置１０は、軌道推定にかかる時間が短いため、利用者が電源をオフ状態にする前に推定を完了し、推定結果をＧＰＳ受信部内の記憶部に記憶させておくことが可能になる。そして、このように記憶された推定結果を、電源オフ状態からオン状態になった場合に用いて推定を行えるので、本実施形態の軌道予測装置１０は、電源オン状態後に位置測位開始までの時間を短縮できる。
また、従来技術と同数のエフェメリス・データを用いて軌道予測を行った場合、推定に用いるパラメータの数を削減したため、観測方程式の数に対して推定パラメータの数が減ったため、パラメータ間の共分散が減少し，精度の良いパラメータ値を推定できる。この結果、本実施形態の軌道予測装置１０は、推定精度を改善することができる。なお、推定しないパラメータに関しては、実測により、軌道予測に十分な精度のデータであることを確認している。
【００８８】
［第２実施形態］
第１実施形態では、推定パラメータの内、９つのパラメータを推定する例を説明した。本実施形態では、第１実施形態で用いた９個の成分のうち、さらに、太陽輻射圧パラメータのワイバイアスとスケールファクタの２つを推定せず、７個の成分のベクトルによる推定パラメータｓを推定する例を説明する。
なお、太陽輻射圧パラメータのワイバイアスとスケールファクタは、ＧＰＳ衛星２毎に固有の固定値を、予めＧＰＳ衛星２毎に算出して軌道予測装置１０内のいずれかの記憶部か機能部に、予め記憶させておく。例えば、所定の期間のデータを算出して、平均値を記憶させるようにしてもよい。
【００８９】
次に、過去のＧＰＳ衛星２の位置と速度データを記憶させる処理について、図２、図５及び図６を用いて説明する。
ＧＰＳ衛星２の位置と速度情報の記憶と削除は、第１実施形態と同様に行う（図２、ステップＳ１〜Ｓ８）。
第１実施形態と異なるのは、衛星位置ｐ０（ｉ）の軌道計算の方法である。
（ステップＳ２０１’）軌道予測演算部１５５は、推定パラメータｓを使用し、時刻ｔ（ｋ）での衛星位置ｐ０（ｉ）（ｉ、・・・、３ｎ−１）を軌道計算する。
推定パラメータｓは、以下の９個の成分のベクトルである。
成分１ ｓ（０）；時刻ｔ（０）での衛星位置のｘ成分（初期値ｒ（０））
成分２ ｓ（１）；時刻ｔ（０）での衛星位置のｙ成分（初期値ｒ（１））
成分３ ｓ（２）；時刻ｔ（０）での衛星位置のｚ成分（初期値ｒ（２））
成分４ ｓ（３）；時刻ｔ（０）での衛星速度のｘ成分（初期値ｖ（０））
成分５ ｓ（４）；時刻ｔ（０）での衛星速度のｙ成分（初期値ｖ（１））
成分６ ｓ（５）；時刻ｔ（０）での衛星速度のｚ成分（初期値ｖ（２））
成分７ ｓ（６）；地球自転パラメータの極運動ｘ方向（初期値０）
成分８ ｓ（７）；地球自転パラメータの極運動ｙ方向（初期値０）
成分９ ｓ（８）；地球自転パラメータの地球自転速度変化（初期値０）
【００９０】
時刻ｔ（０）での衛星位置において、半径方向、進行方向、軌道面垂直方向の各単位ベクトルｄ１、ｄ２、ｄ３は、第１実施形態と同様に式（２）〜（４）のように表される。
【００９１】
ステップＳ２０２〜ステップＳ２０５は、第１実施形態と同様に行う。
（ステップＳ２０６’）軌道予測演算部１５５は、ｊ＝４，５，６に対して推定パラメータの１つの成分ｓ（ｊ＋２）をｓ（ｊ＋２）＋Δｓ（ｊ）に代え、時刻ｔ（ｋ）での衛星位置ｐ（ｉ）を式（９）のように軌道計算する。ただし、式（９）おいて、ｊ＝４，５，６である。
【００９２】
（ステップＳ２０７’）軌道予測演算部１５５は、次式（１４）を用いて、推定パラメータｓを算出する。
【００９３】
ｓ＝ｓ＋Ｂ（Ａ×ＡＴ）−１×ＡＴ（ｒ−ｐ０） ・・・（１４）
【００９４】
式（１４）は、右辺の演算を実行して、演算結果を新たなｓとして更新することを表している。なお、式（１４）において、Ｔは転置行列を表し、行列Ｂは、次式（１５）である。また、行列Ｂは、行列Ａの成分を、元のｘｙｚ方向の成分に座標系を変換する役割をしている。
【００９５】
【数６】
【００９６】
（ステップＳ２０８）軌道予測演算部１５５は、推定パラメータｓが収束したか否かを判別する。
推定パラメータｓが収束したと判別した場合（ステップＳ２０８；Ｙｅｓ）、ステップＳ２０９に進む。推定パラメータの１つの成分ｓが収束していないと判別した場合（ステップＳ２０８；Ｎｏ）、ステップＳ２０１’に戻る。
（ステップＳ２０９）軌道予測演算部１５５は、この推定パラメータの１つの成分sを使用し、将来の衛星位置ｐｅを軌道計算する。
軌道予測演算部１５５は、第１実施形態と同様に軌道計算を行う。
【００９７】
以上のように、第２実施形態によれば、第１実施形態と比較して、推定パラメータの成分の内、２つについて固定値を用いたため、パラメータを推定する際に最小二乗法で用いる行列が小さくなる。このように最小二乗法で用いる行列を小さくできた結果、本実施形態の軌道予測装置１０では、演算に用いる必要メモリの容量が小さくて済み、さらに演算時間を短くできる効果がある。演算に必要なメモリ容量を少なくすることができたので、本実施形態の軌道予測装置１０を用いた測位装置１では、コストダウンを実現できる。また、演算時間を短くできることにより、測位装置１の電源をオフ状態からオン状態にしたような場合、従来技術より早く起動後に位置測定を開始できるので、本実施形態の軌道予測装置１０では、利用者にとって利便性が増すため、利用者の使用頻度が向上する効果もある。
また、第２実施形態によれば、第１実施形態と比較して、さらに推定パラメータにおける未知数が少なくなるので、本実施形態の軌道予測装置１０では、必要な観測方程式の数を減らせる効果もある。この結果、少数のエフェメリス・データを用いて、軌道予測を行うことが可能になるので、従来技術と比較して、本実施形態の軌道予測装置１０では、利用者にとって利便性が増すため、利用者の使用頻度が向上する効果もある。
【００９８】
なお、本実施形態では、３日分の位置速度情報を衛星位置・速度記憶部１４０に記憶させる例を説明したが、本実施形態の軌道予測装置１０を用いる測位装置１の使用目的に応じて、記憶させる期間を３日より多くしても３日より少なくしてもよい。また、エフェメリス・データを記憶するタイミングとして２時間間隔の例を説明したが、測位装置１の用途に応じて、エフェメリス取得部１３０は、ＧＰＳ衛星２から送信されたエフェメリス・データを２時間毎に受信し、受信した後に、４時間間隔に記憶させるようにして間引くようにしてもいい。また、衛星位置・速度記憶部１４０に記憶させる位置速度情報が１０個の例を説明したが、測位装置１の用途に応じて、１０個以上であってもよい。
【符号の説明】
【００９９】
１・・・測位装置、２・・・ＧＰＳ衛星、１０・・・軌道予測装置、２０・・・測位算出部、１１０・・・情報取得部、１１５・・・記憶部、１２０・・・アルマナック取得部、１２５・・・衛星型抽出部、１３０・・・エフェメリス取得部、１３５・・・衛星位置速度算出部、１４０・・・衛星位置・速度記憶部、１４５・・・演算開始判定部、１５０・・・パラメータ推定部、１５５・・・軌道予測演算部、１６０・・・予測結果記憶部
【特許請求の範囲】
【請求項１】
測位衛星からの信号を受信して受信データを取得する情報取得部と、
前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出部と、
軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽から前記測位衛星への光の圧力である太陽輻射圧の影響を示す値、及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ部と、
前記推定パラメータ部で推定された推定パラメータを用いて前記測位衛星の位置を推定する軌道予測部と、
を備えることを特徴とする軌道予測装置。
【請求項２】
前記推定パラメータ部は、
前記測位衛星の半径方向と進行方向及び軌道面垂直方向との成分による座標系に変換する行列を用いて前記推定パラメータを推定する
ことを特徴とする請求項１に記載の軌道予測装置。
【請求項３】
前記推定パラメータ部は、
前記測位衛星の位置の進行方向成分に対して予め算出された値を用い、前記測位衛星の速度の軌道面垂直方向成分に対して予め算出された値を用い、前記推定パラメータを推定する
ことを特徴とする請求項１または請求項２に記載の軌道予測装置。
【請求項４】
前記推定パラメータ部は、
前記測位衛星に対する太陽輻射圧のワイバイアスとスケール係数とに対して各々予め算出した値を用い、前記推定パラメータを推定する
ことを特徴とする請求項１から請求項３のいずれか１項に記載の軌道予測装置。
【請求項５】
軌道予測装置における軌道予測方法であって、
情報取得部が、測位衛星からの信号を受信して受信データを取得する情報取得過程と、
衛星位置速度算出部が、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出過程と、
推定パラメータ部が、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽輻射圧の影響を示す値及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ過程と、
軌道予測部が、前記推定パラメータ手順で推定された前記推定パラメータを用いて前記測位衛星の位置を推定する軌道予測過程と、
を含むことを特徴とする軌道予測方法。
【請求項１】
測位衛星からの信号を受信して受信データを取得する情報取得部と、
前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出部と、
軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽から前記測位衛星への光の圧力である太陽輻射圧の影響を示す値、及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ部と、
前記推定パラメータ部で推定された推定パラメータを用いて前記測位衛星の位置を推定する軌道予測部と、
を備えることを特徴とする軌道予測装置。
【請求項２】
前記推定パラメータ部は、
前記測位衛星の半径方向と進行方向及び軌道面垂直方向との成分による座標系に変換する行列を用いて前記推定パラメータを推定する
ことを特徴とする請求項１に記載の軌道予測装置。
【請求項３】
前記推定パラメータ部は、
前記測位衛星の位置の進行方向成分に対して予め算出された値を用い、前記測位衛星の速度の軌道面垂直方向成分に対して予め算出された値を用い、前記推定パラメータを推定する
ことを特徴とする請求項１または請求項２に記載の軌道予測装置。
【請求項４】
前記推定パラメータ部は、
前記測位衛星に対する太陽輻射圧のワイバイアスとスケール係数とに対して各々予め算出した値を用い、前記推定パラメータを推定する
ことを特徴とする請求項１から請求項３のいずれか１項に記載の軌道予測装置。
【請求項５】
軌道予測装置における軌道予測方法であって、
情報取得部が、測位衛星からの信号を受信して受信データを取得する情報取得過程と、
衛星位置速度算出部が、前記受信データに基づいて、前記測位衛星の位置と前記測位衛星の速度である位置速度情報を算出する衛星位置速度算出過程と、
推定パラメータ部が、軌道における前記測位衛星の半径方向、進行方向及び軌道面垂直方向の各単位ベクトルと、前記位置速度情報とを用いて、前記測位衛星の位置と、前記測位衛星の速度と、太陽輻射圧の影響を示す値及び地球自転の影響を示す値とを含む推定パラメータを算出する推定パラメータ過程と、
軌道予測部が、前記推定パラメータ手順で推定された前記推定パラメータを用いて前記測位衛星の位置を推定する軌道予測過程と、
を含むことを特徴とする軌道予測方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２０１３−３０９３（Ｐ２０１３−３０９３Ａ）
【公開日】平成２５年１月７日（２０１３．１．７）
【国際特許分類】
【出願番号】特願２０１１−１３７３２６（Ｐ２０１１−１３７３２６）
【出願日】平成２３年６月２１日（２０１１．６．２１）
【出願人】（０００００４３３０）日本無線株式会社 (1,186)
【Ｆターム（参考）】
【公開日】平成２５年１月７日（２０１３．１．７）
【国際特許分類】
【出願日】平成２３年６月２１日（２０１１．６．２１）
【出願人】（０００００４３３０）日本無線株式会社 (1,186)
【Ｆターム（参考）】
[ Back to top ]