JP2005062046A - 人工衛星の位置推定システム - Google Patents

Abstract

【課題】 人工衛星の軌道を推定するにあたり、観測ノイズの分散を小さくし、高精度な位置推定を行うこと。
【解決手段】 人工衛星から周波数ｆ_i（波長λ_i）の搬送波で送信される測位信号を観測して得た擬似距離観測データＰＲｉ^s _r(t)及び搬送波位相観測データＣＰｉ^s _r(t)と、前記測位信号と同時に周波数ｆ_j（波長λ_j）の搬送波で送られる測位信号とを観測して得た擬似距離観測データＰＲｊ^s _r(t)及び搬送波位相観測データＣＰｊ^s _r(t)と、を下記数１に代入して第１の合成値及び第２の合成値を算出し、これら合成値に基づいて人工衛星の位置を推定する。
【選択図】 図１

Description

本発明は、人工衛星から送信される測位信号の観測データから測位衛星等の人工衛星の位置を推定する人工衛星の位置推定システムに関する。

米国のＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）、ロシアのＧＬＯＮＡＳＳ、ヨーロッパのＧＡＬＩＲＥＯ、更に日本で計画中の準天頂衛星システム等の測位衛星システムを利用して精密な測定位置及び時刻同期を行うには、測位衛星が地球の周りを周回する軌道及び測位衛星に搭載した時計の誤差を高精度に推定することが必要不可欠である。

測位衛星から送信される測位信号を利用して測位衛星の軌道及び測位衛星に搭載した時計誤差を推定する際に必要となるモデルは、測位衛星の運動を表す運動（力学）モデルと測位信号を受信した観測データに含まれる各種の誤差を表現する観測モデルとが必要である。この測位信号は測位衛星と観測点間との距離を測るために使用されるが、測位信号の観測データの中には測位衛星と観測点間との距離データの他に測位衛星に搭載された時計の時間誤差、電離層や対流圏による電波遅延、受信機ノイズなど様々なデータ誤差が含まれる。これらのデータ誤差を補正及び除去することが高精度の測位衛星の軌道及び測位衛星に搭載した時計誤差の推定では必要となる。

ここで図６を用いて測位衛星の軌道を推定する原理について説明する。地球の大気圏外領域を回っている測位衛星１からは測位信号Ｚと呼ばれる一定の信号形式を持つ電波が常時、測位衛星１に搭載された衛星時計１１に同期して送信機１２によって送信される。この測位信号Ｚは、地上のモニタ局２Ａ〜２Ｃに設置された受信機３Ａ〜３Ｃによって受信され、以下の原理で測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を測定する。

先ず、測位衛星１が測位信号Ｚを送信した時刻は、各モニタ局２Ａ〜２Ｃで受信した信号を解読することにより知ることができる（送信時刻（ｔ））。モニタ局２Ａ〜２Ｃで測位信号Ｚを受信した時刻はモニタ局２Ａ〜２Ｃに設置されたモニタ局時計４Ａ〜４Ｃで知ることができる（受信時刻（Ｔ））。このことから測位信号Ｚが測位衛星１からモニタ局２Ａ〜２Ｃまでの空間を伝播した時間（伝播時間（Δｔ））はΔｔ＝Ｔ−ｔとなり、よって測位信号Ｚは光速ｃで空間を伝わることから測位衛星１とモニタ局２Ａ〜２Ｃ間との距離ＤはＤ＝Δｔ×ｃで表すことができる。

このように予め正確な位置が分っている複数のモニタ局２Ａ〜２Ｃからモニタ局２Ａ〜２Ｃと測位衛星１との各距離データＤ１〜Ｄ３を主管制局５に設けられた演算部５１で演算することにより測位衛星１の位置を推定することができる。このシステムにより得られた測位衛星１の位置は、例えば地上の移動体のナビゲーションシステムにも適用される。

上述した測位衛星１の運動を表す運動（力学）モデルは、地球の引力等の物理法則に支配されており、ある時刻での位置、速度が分れば一定時間後の位置、速度を予測することができるので、測位衛星１の軌道推定では、測位信号Ｚを正確に測定することが重要となる。

そこで、測位信号Ｚを利用した測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を測定する方法には、測位信号Ｚの信号形式（一定の長さを持ったパターン）を利用した擬似距離を用いる方法と、測位信号Ｚの搬送波を数える搬送波位相を用いる方法とがある。

この擬似距離を用いる方法では、測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を直接測ることができるが、計測精度が信号パターンを構成するデータの単位長さに制限されるため距離の誤差が大きくなる。一方前記搬送位相を用いる方法では、計測精度は測位信号Ｚの波長以下となり測定精度を上げることができるが、受信電波の波を数えるため波長の整数倍の不確定性（アンビギュイティ）が残ってしまう。

また、前記測位衛星１からは、詳細には、異なった周波数で複数の測位信号Ｚが同時に送信されており、これら複数の測位信号Ｚにおいて、各々の測位信号Ｚ毎に擬似距離及び搬送波位相が求められる。しかしながら、測位衛星１から送信された測位信号Ｚを地上のモニタ局２Ａ〜２Ｃに設置された受信機３Ａ〜３Ｃで受信すると測位衛星１と地上のモニタ局２Ａ〜２Ｃ間との距離データＤ１〜Ｄ３の中には（表１）に示す各種誤差が混入する。

これらのモニタ局２Ａ〜２Ｃで観測された擬似距離観測データと、搬送波位相観測データと、に含まれる観測誤差を表した観測モデルは、一般的に数３のように表される。

ここで、上述の数３を用いて測位衛星１の軌道決定のための推測処理をどのようにして行うのかを説明する。数３の式の左辺は観測値であり、例えば数３中のＰＲｉ^ｓ _ｒ（ｔ）は擬似距離観測値である。この擬似距離観測値は次のようにして求められる。測位衛星１は、例えば地上から２万ｋｍの高度の所を移動しており、この測位衛星１からパルスパターンを含む例えば長さ３００〜５００ｋｍ程度のパルス列を地上に送る。このパルスパターンは、パルス列の先頭からの距離に対応したパターンとして配列され、例えばパルス列の先頭からＡｃｍの位置にあればＡｃｍを表すパルスパターンが形成される。つまり、このパルス列は測位衛星１からいわば長さ３００〜５００ｋｍの物差しが飛来してくるのと同じである。

一方測位衛星１と地上のモニタ局２Ａ〜２Ｃとにおいては、正確に時計合わせされたモニタ局時計４Ａ〜４Ｃが夫々設けられ、測位衛星１からパルス列を発信するタイミング及びモニタ局２Ａ〜２Ｃでこのパルス列を検出するタイミングが予め決められている。従ってモニタ局２Ａ〜２Ｃの検出のタイミングでパルス列のあるパルスパターンを把握すれば、その物差しのいくつの値をどのタイミングで検出したかが分る。測位衛星１とモニタ局２Ａ〜２Ｃとの距離は運動（力学）モデルにより概ね分っているので（例えば数ｍの誤差で分っているので）、この測位衛星１とモニタ局２Ａ〜２Ｃとの距離に比べてパルス列の長さが短くても、当該パルス列の中のパルスパターンをある時刻で検出すれば、そのパルスパターンはパルス列を物差しと見立てたときのスケール値に相当するので、当該パルス列の先頭からどれだけの時間だけ遅れて送信されているかが分り、結局パルス列の先頭をモニタ局２Ａ〜２Ｃが受信した時刻が分る。この結果、パルス列の先頭が飛来してくる時間（（モニタ局２Ａ〜２Ｃで受信した時刻）−（測位衛星１から送信された時刻））が検出され、この時間に光速ｃを掛けることで測位衛星１とモニタ局２Ａ〜２Ｃとの距離を算出することができる。しかしながら、実際には数３から分るようにこの距離は、種々の誤差を含んでいる。

また上記の数３中のＣＰｉ^ｓ _ｒ（ｔ）を用いても同様にして測位衛星１とモニタ局２Ａ〜２Ｃとの距離ρ^ｓ _ｒ（ｔ）を用いることができる。この場合には、擬似距離の観測データを用いる代わりに搬送波の位相を観測した観測データを用いている。この搬送波とは、測位衛星１から送られてくる既述のパルス列の搬送波であり、搬送波位相観測値ＣＰｉ^ｓ _ｒ（ｔ）とは、ある時刻ｔ_０からある時刻ｔまでに地上のモニタ局２Ａ〜２Ｃで搬送波をカウントし、搬送波をカウントした値から時刻ｔ_０における測位衛星１の位置と時刻ｔにおける測位衛星１の位置とが変わることに基づいて搬送波の波数値を求めたものである。この搬送波位相による距離計測の原理について図７と数４を参照しながら詳しく説明する。

先ず図７で、モニタ局２に設けられた受信機３で搬送波を受信開始したモニタ局時計４の時刻をｔ_０とする。この時点で受信した搬送波が測位衛星１に設けられた送信機１２から送信されたときの衛星時計１の時刻は、測位衛星１とモニタ局２間との距離をｒ（ｔ_０）とするとｔ_０−｛ｒ（ｔ_０）／ｃ｝となり、この時点では受信機３で受信した受信波数は０個である。この時点から同時に受信機３内部の図示されない局部発振機の波数積算も開始される。ここで搬送波の計測時刻をｔとする。この時点で受信機３で受信した搬送波が測位衛星１に設けられた送信機１２から送信されたときの衛星時計１１の時刻は、測位衛星１とモニタ局２間との距離をｒ（ｔ）とするとｔ−｛ｒ（ｔ）／ｃ｝となり、受信開始時刻から計測時刻までに受信した波数Φ（ｔ）は数４で表される。

一般に静止系であれば、所定時間に受信する電波の数は一定であり、数４の右辺の前段の項である１／２πω・（ｔ−ｔ_０）だけで決定されるが、測位衛星１は移動しているため、後段の項である{ｒ（ｔ）／λ}−{ｒ（ｔ_０）／λ}が含まれてくる。この値が搬送波の波数値であり、ＣＰｉ^ｓ _ｒ（ｔ）に相当するものである。しかしながら、上述したようにこの搬送波の位相を利用する方法では、計測精度は測位信号Ｚの波長以下となり高精度であるが、受信電波の波数を数えるため波長の整数倍の不確定性（アンビギュイティ）が残ってしまう。

測位衛星１の軌道推定は、測位衛星１と地上のモニタ局２Ａ〜２Ｃ間との距離データＤ１〜Ｄ３を利用して測位衛星１の位置、速度を求めることであるが、測位信号Ｚの観測データ中には測位衛星１とモニタ局２Ａ〜２Ｃ間との距離以外の誤差が含まれているため、これらの誤差の値も同時に推定する必要が生じてくる。しかし、前記数３で表した観測モデルのままでは未知数が多く夫々の観測誤差を分離して推定することができない。そこで擬似距離観測データと搬送波位相観測データとを次のように合成した擬似的な観測データを利用して未知数を減らす方法が一般的な手法である（例えば、非特許文献１参照。）。この一般的な手法を図８と数５を参照しながら詳しく説明する。

この方法は、図８に示すように２つのモニタ局ｒ（ｒ´）と二つの測位衛星ｓ（ｓ´）とを用い、更に各測位衛星ｓ（ｓ´）のいずれからも送信される周波数ｆ_ｉの搬送波と周波数ｆ_ｊの搬送波とを利用する方法（二重差観測データ）であり、先の数３で表される擬似距離観測値と搬送波位相の波数値とを数５のように定義する。これらの式においてＣＰｉ^ｓ _ｒ（ｔ）は、モニタ局ｒにおいて、測位衛星ｓから送信される周波数ｆ_ｉの搬送波について数４のように時刻ｔ〜ｔ_０間で搬送波をカウントし、数４中のωは既値であることから、後段の項を求めた値である。また、ＣＰｊ^ｓ _ｒ´（ｔ）は、モニタ局ｒ´において測位衛星ｓから送信される周波数ｆ_ｊの搬送波について求めた値である。更にＰＲｉ^ｓ _ｒ（ｔ）は、測位衛星ｓから周波数ｆ_ｉの搬送波に乗せられて送られる既述のパルス列中のパルスパターンをモニタ局ｒについて検出し、その検出結果に基づいて得られた擬似距離観測値である。

そして、数５の右辺の各項に、数３で定義される値（右辺）を代入することにより測位衛星ｓ（ｓ´）と地上のモニタ局ｒ（ｒ´）間との距離以外の誤差要因（観測誤差）は時計誤差、電離層遅延及び周波数間バイアスが合成により打消されて、対流圏遅延、搬送波位相アンビギュイティ及び観測ノイズのみが残ることになる（数６）。

もう少し詳しく説明すると、数５の右辺は、各モニタ局ｒ（ｒ´）における観測値を演算した値であるから、各モニタ局ｒ（ｒ´）における観測値{ＰＲｉ^ｓ _ｒ（ｔ）等}を地上回線で主管制局のコンピュータに取り込むことによって、数５の左辺であるＩＣＰｉｊｒ^ｓｓ´ _ｒｒ´（ｔ）及びＩＡＭｉｊｒ^ｓｓ´ _ｒｒ´（ｔ）が求まる。

今、時刻ｔ_ｋにおける擬似距離観測値と搬送波の波数値との各観測値を用いて数５を演算した結果を夫々ＩＣＰｉｊｒ^ｓｓ´ _ｒｒ´（ｔ_ｋ）及びＩＡＭｉｊｒ^ｓｓ´ _ｒｒ´（ｔ_ｋ）とする。なお、時刻ｔ_ｋにおける搬送波の波数値は、時刻ｔ_０からの搬送波の積算値に基づいて求められるが、ω（ｔ_ｋ−ｔ_０）の値は既にわかっているので、この値をキャンセルすることにより、結局、数４の{ｒ（ｔ_ｋ）／λ}−{ｒ（ｔ_０）／λ}のみが波数値として求まることになる。

一方コンピュータは、時刻ｔ_ｋよりも前の時刻ｔ_ｋ−１において推定した数６の右辺の測位衛星ｓ（ｓ´）とモニタ局ｒ（ｒ´）間との距離及び推定状態量に基づいて時刻ｔ_ｋの右辺の値を求める。ここで測位衛星ｓ（ｓ´）とモニタ局ｒ（ｒ´）間との距離であるρ^ｓ _ｒ（ｔ）については初期値としては適当な値を用いる。但し、測位衛星ｓ（ｓ´）の軌道は予めかなりの精度で把握されているので、適当な値とはいっても高い精度を持つ値である。また、他の推定状態量のうちＴ^ｓ _ｒ（ｔ）は信頼できる理論的モデル（時間の関数）を用い、時刻ｔ_ｋ−１のときの対流圏遅延を用いて、時刻ｔ_ｋのときの対流圏遅延（Ｔ^ｓ _ｒ（ｔ））を求めることができる。またＮｉ^ｓ _ｒは連続して電波受信中は定数として考えることができ、更にＥｉ^ｓ _ｒ（ｔ）、ｅｉ^ｓ _ｒ（ｔ）はいずれも平均値は０と仮定して取り扱うことができる。こうして数６を用いて、数５の右辺が求まる。先に求めたＩＣＰｉｊ^ｓｓ´ _ｒｒ´（ｔ）及びＩＡＭｉｊ^ｓｓ´ _ｒｒ´（ｔ）は数５の左辺であるため、本来は左辺の観測値と右辺の演算値とは一致するはずであるが、実際にはモデルの誤差分や種々の仮定が含まれているため、また時刻ｔ_ｋの値があくまでも推定値であるため、観測値と演算値とは異なってくる。

そこで、コンピュータは数５の夫々について左辺と右辺との差（偏差）を求める。そして、これら偏差に応じて時刻ｔにおけるρ^ｓ _ｒ（ｔ）、Ｔ^ｓ _ｒ（ｔ）及びＮｉ^ｓ _ｒ（ｔ）の値を修正する。実際には、主管制局では各モニタ局毎でこれらの推定値が求まり、時刻ｔ_ｋにおける推定値としてメモリ内に格納され、次の時刻ｔ_ｋ＋１における推定に用いられる。また、推定初期値においては、誤差が大きい場合でも時間が経過するにつれて精度が高くなり前記誤差は小さくなっていく。時刻ｔ_ｋにおける測位衛星の位置、速度は３ヶ所以上のモニタ局と測位衛星との距離から求まることになる。

Global Positioning System:Theory and Applications Volume1 Part2(American Institute of Aernautics and Astronautics 1996)

しかしながら、擬似距離観測データと搬送波位相観測データとを上述のように合成して未知数を減らす一般的な手法（二重差観測データ）では、二重差の使用により測位衛星の軌道は決定可能となる８つの観測データを合成することにより観測ノイズの分散は１つの観測データの場合に比べて８倍となり精密な軌道決定には不利である。また近年、この二重差の不利を補うためにＧＰＳの軌道決定において、一般的に全世界の観測網を利用して長時間連続的に測位衛星の観測及び軌道推定を行うことで除々に高精度な軌道に収束させる方法が用いられているが、地上に多数のモニタ局が必要となるため、多大なコストが掛かるという問題がある。

本発明は係る事情に鑑みてなされたものであって、本発明の目的は、二重差を使用せずに、測位衛星から送信される測位信号に含まれる観測誤差（観測ノイズ）を小さく抑え、測位衛星等の人工衛星の位置の推定を可能とする人工衛星の位置推定システムを提供することにある。

本発明は、人工衛星から送信され、その先頭からの位置に対応する位置コードが配列された所定の長さのパルス列を搬送波に載せた測位信号を地上のモニタ局にて受信し、モニタ局にて検出した位置コードと光速とに基づいて人工衛星とモニタ局との擬似距離観測データを求めると共に、
モニタ局にて検出した前記搬送波の波数を所定の時刻ｔ₀から前記位置コードをサンプリングした時刻ｔ_kまでの時間をカウントし、そのカウント値と前記搬送波の周波数と前記カウントした時間とに基づいて、時刻ｔ₀におけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値と時刻ｔ_kにおけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値との差分である搬送波位相観測データを求め、
前記擬似距離観測データと前記搬送波位相観測データとに基づいて、人工衛星とモニタ局との距離を求め、当該人工衛星と複数のモニタ局の各々との距離に基づいて当該人工衛星の位置を求めるシステムにおいて、
イ． 人工衛星から送信される周波数ｆ_iの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｉ^s _r(t)及び搬送波位相観測データＣＰｉ^s _r(t)と、
前記測位信号と同時に送られる、周波数ｆ_jの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｊ^s _r(t)及び搬送波位相観測データＣＰｊ^s _r(t)と、
を下記数７に代入して第１の合成値及び第２の合成値を算出する合成値演算部と、
ロ． 時刻ｔ_kに観測した観測データに基づいて計算された第１の合成値及び第２の合成値と、
数７を変形して得られた下記数８と、
時刻ｔ_kのサンプリングよりも一つ前のサンプリングのタイミングｔ_k-1にて観測した観測データを用いて計算された第１の合成値及び第２の合成値に基づいて推定した人工衛星の位置と、
に基づいて、時刻ｔ_kにおける人工衛星の位置を推定する推定部と、を備えたことを特徴とする人工衛星の位置推定システムである。

記

ただし、λ_i及びλ_jは夫々周波数ｆ_iの搬送波の波長及び周波数ｆ_jの搬送波の波長である。

ただし、ρ^s _r (t)は人工衛星及びモニタ局の距離、Ｔ^s _r(t)は電波の対流圏遅延に起因する誤差、Ｄｔ_r(t)はモニタ局の時計誤差と当該モニタ局の受信機における周波数ｆ_iの搬送波及び周波数ｆ_jの搬送波の受信のタイミングの誤差との加算値に起因する誤差、Ｄｔ^s(t)は人工衛星の搭載時計の誤差と当該人工衛星における周波数ｆ_iの搬送波及び周波数ｆ_jの搬送波の発信のタイミングの誤差との加算値に起因する誤差、Ｎ^s _r は搬送波位相アンビギュイティ、ｅ^s _r (t)は搬送波位相観測ノイズ、Ｅ^s _r (t)は擬似距離観測ノイズである。

前記推定部は、時刻ｔ_kのサンプリングよりも一つ前のサンプリングのタイミングｔ_k−1にて観測した観測データを用いて計算された第１の合成値及び第２の合成値に基づいて推定した前記Ｄｔ_r(t)及びＤｔ^s(t)を更に用いて、時刻ｔ_kにおける人工衛星の位置を推定することが好ましい。また前記推定部は、第１の合成値及び第２の合成値の分散と、人工衛星の位置を含む状態量の各共分散と、を更に用いて、時刻ｔ_kにおける人工衛星の位置を推定することが好ましい。

ここで数７と数８との関係について説明しておく。数７の右辺の各項に数３の右辺を代入することにより数９が得られる。

数９の｛ ｝内を数１０のように合成された１つの未知数として取り扱うと数８が得られる。

数８において擬似距離観測ノイズＥ^s _r （ｔ）は平均値がゼロと仮定することによりゼロとして取り扱うことができ、従って下段の式から搬送波位相アンビギュイティＮ^s _r の推定値が求まる。また搬送波位相観測ノイズｅ^s _r (t)も平均値がゼロと仮定することによりゼロとして取り扱うことができる。このため残る状態量は、人工衛星とモニタ局との距離ｐ^s _r (t)、電波の対流圏遅延に起因する誤差Ｔ^s _r(t)、モニタ局の時計誤差と当該モニタ局の受信機における周波数ｆ_iの搬送波及び周波数ｆ_jの搬送波の受信のタイミングの誤差（周波数間バイアス）との加算値に起因する誤差Ｄｔ_r(t)、人工衛星の搭載時計の誤差と当該人工衛星における周波数ｆ_iの搬送波及び周波数ｆ_jの搬送波の発信のタイミングの誤差（周波数間バイアス）との加算値に起因する誤差Ｄｔ^s(t)である。

これら状態量の中で対流圏遅延は例えば種々の周知の理論的モデルを用い、モニタ局と人工衛星との相対位置に応じた定数として取り扱うことができる。従って時刻ｔ_kにおける人工衛星の位置を推定する推定部においては、例えば人工衛星の位置（Ｘ，Ｙ，Ｚの座標軸の各位置）、Ｄｔ_r(t)及びＤｔ^s(t)とこれら状態量の各々の分散（共分散）とにおいて、時刻ｔ_kのサンプリングよりも一つ前のサンプリングのタイミングｔ_k-1に対応する値、つまりこれから状態量を推定するために行う今回の演算のサイクルの一つ前の演算サイクルにて推定された推定値を用い、更に時刻ｔ_kに観測した観測データに基づいて計算された第１の合成値及び第２の合成値とこれら合成値の分散を用いて、今回の演算サイクルにおける状態量及び分散を推定する。なお第１の合成値及び第２の合成値の分散は定数として扱うことができる。

また前記推定部は、人工衛星とモニタ局との組毎に各状態量を推定するので、推定演算の過程では例えば人工衛星とモニタ局との距離を推定し、各人工衛星の位置についてはその人工衛星と３つ以上のモニタ局との距離から計算することで推定している。こうして各人工衛星の位置を推定することができ、例えば推定された位置と一つ前の演算サイクルで求めた位置とから人工衛星の速度も推定できる。

ＧＰＳで行われている二重差による手法（数５の手法）では、擬似距離及び搬送波位相の各々について８つの観測データを合成していたが、本発明では２つの観測データを合成した合成値に基づいて測位衛星等の人工衛星の位置が推定でき、従って観測ノイズの分散は二重差の手法に比べて４分の１となり、高精度な位置推定を行うことができる。更に測位衛星等の人工衛星の時計誤差及びモニタ局の時計誤差の推定も高精度に行うことができる。

本発明の人工衛星の位置推定システムについて、図１〜図３に基づいて説明を行う。図１は本発明の実施の形態に係る人工衛星である測位衛星の位置測定システムを含む測位衛星の軌道推定システムの説明図である。図１中６は測位衛星であり、６Ａ、６Ｂ、６Ｃは３機の測位衛星の夫々を指している。前記測位衛星６は搭載時計６１、受信機６２、測位信号生成部６３及び送信機６４を備えている。前記測位信号生成部６３は、搭載時計６１からクロック信号に基づいて、周波数ｆ_ｉ（例えば１５０ＧＨｚ）の搬送波で送られる測位信号であるパルス列と周波数ｆ_ｊ（例えば１７０ＧＨｚ）の搬送波で送られる測位信号であるパルス列とを作成して前記送信機６４に受け渡す機能を備えている。これらパルス列は例えば３００ｋｍ〜５００ｋｍ程度の長さに作られており、パルス列からの先頭からの位置に対応するコードが配列されている。図２はパルス列を示すイメージ図であり、先頭には当該パルス列を発信する測位衛星６の識別コード（ＩＤ）が配置されると共に、その後には各位置に対応したコード（Ｐ_１、Ｐ_２、…Ｐ_ｎ）が配置されている。各コードの長さは例えば２０ｃｍ程度とされる。

なお、測位衛星６は後述の主管制局８から、Ｘ、Ｙ、Ｚ座標における自己の正確な位置（Ｘ、Ｙ、Ｚ）と速度（Ｖ_ｘ、Ｖ_ｙ、Ｖ_ｚ）とを含む航法メッセージを受け取り、その航法メッセージを測位信号と共に一般ユーザ１００の受信機１０１に配信する機能を備えている。前記一般ユーザ１００は、例えば車両等の移動体であり、３機の測位衛星６Ａ〜６Ｃから送られてくる航法メッセージと測位信号とに基づいて、各測位衛星６Ａ〜６Ｃまでの各距離が分るので、当該一般ユーザ１００（移動体）の位置を知ることができる。

また、図１中７は地上に設けられたモニタ局であり、７Ａ〜７Ｃは３個のモニタ局を夫々示している。モニタ局７は受信機時計７１及び受信機７２を備えており、受信機時計７１から出力されるクロック信号に基づいて測位衛星６から送られてくる測位信号（パルス列及び搬送波）を検出し、周波数ｆ_ｉの測位信号に基づいて擬似距離観測データＰＲｉ^ｓ _ｒ（ｔ）及び搬送波位相観測データＣＰｉ^ｓ _ｒ（ｔ）を求め、また、周波数ｆ_ｊの測位信号に基づいて擬似距離観測データＰＲｊ^ｓ _ｒ（ｔ）及び搬送波位相観測データＣＰｊ^ｓ _ｒ（ｔ）を求める機能を備えている。更に各測位衛星６の搭載時計６１及び各モニタ局７の受信機時計７１は正確に時刻が合わせられている。

擬似距離観測データＰＲｉ^ｓ _ｒ（ｔ）の求め方については、背景技術で既に述べたように、サンプリングのタイミングで検出したパルス列のコードに基づいて測位衛星６とモニタ局７間との距離が分る。即ち各測位衛星６の軌道はかなり正確に分っていることから、パルス列の長さが測位衛星６とモニタ局７間との距離（例えば２万ｋｍ程度）に比べて短くても、当該測位信号が何時に送信されたものであるかは把握されているので、パルス列の先頭からの長さを表すコードを検出した時刻と当該測位信号が送信された時刻との差分（時間）と光速とに基づいて測位衛星６とモニタ局７間との距離が計算できる。この距離は、既述のように数３で表されるいわば擬似距離である。これは種々の誤差を含むため、主管制局８に送られてデータ処理され、正確な距離が算出されることになる。

搬送波位相観測データＣＰｉ^ｓ _ｒ（ｔ）は、背景技術で既に述べたよう所定の時刻ｔ_０からカウントした周波数ｆ_ｉの搬送波のカウント値を求め、このカウント値から{周波数ｆ_ｉ×（ｔ−ｔ_０）}を差し引いた値である。先の数４に示すように搬送波位相観測データＣＰｉ^ｓ _ｒ（ｔ）は、その差分であるｒ（ｔ）／λ−ｒ（ｔ_０）／λであり、仮に測位衛星６とモニタ局７との相対位置が変わらなければゼロである。時刻ｔ_０は受信機７１が搬送波と同調して検出を開始したタイミングであってもよいが、通常は測位信号は連続して発信されるため、受信機７１側は連続して搬送波を検出しており、従ってある所定時刻をｔ_０として、そこから積算した値を用いてもよい。更に搬送波位相観測データＣＰｉ^ｓ _ｒ（ｔ）は受信機７２内の図示されない局部発信器を用いて、前記局部発信器から出力される信号を利用してもよい。この場合、搬送波受信開始（カウント開始）から計測時刻までの間に局部発信器が発信する波数は次のΦ´（ｔ）で表される。

Φ´＝ｆ´（ｔ−ｔ_０）＝｛ω´（ｔ−ｔ_０）｝／２πｆ´＝ω´／２π但しω´は受信機７２の局部発振器の角周波数、ｆ´は受信機
７２の局部発振器の周波数である。

また、計測時刻ｔにおける受信搬送波と受信機７２の局部発振器が発生した信号との積は次のように表される。

受信機７２は下式の右辺の高周波部がフィルタで除去し、ω、ω´、ｔ、ｔ_０は既知であることから次の項を記録する。下式の右側の項であるｒ（ｔ）／λ−ｒ（ｔ_０）／λが搬送波位相のアンビギュイティとなる。

また、カウント値を読み出すタイミングは、上記のパルス列のコードを検出するタイミングであり、例えばそのタイミングを時刻ｔ_ｋとすると、時刻ｔ_ｋにおける擬似距離観測データＰＲｉ^ｓ _ｒ（ｔ）及び搬送波位相観測データＣＰｉ^ｓ _ｒ（ｔ）が観測されたことになる。

図１中８は主管制局であり、この主管制局８は、観測データ記憶部８１と、合成値演算部８２と、合成値記憶部８３と、状態量推定部８４と、推定状態量記憶部８５と、航法メッセージ作成部８６と、送信機８７と、を備えている。前記観測データ記憶部８１は、各モニタ局７から送られてくる観測データを記憶する記憶部であり。各モニタ局７Ａ〜７Ｃ毎に各測位衛星６Ａ〜６Ｃからの観測データ（３×３＝９組の観測データ）が記憶される。なお、ここでは説明の便宜上測位衛星とモニタ局との数を夫々「３つ」としているが、実際のシステムはこの数に限定されるものではない。前記合成値演算部８２は、各観測データの組毎に、ＣＰｉ^ｓ _ｒ（ｔ）、ＣＰｊ^ｓ _ｒ（ｔ）、ＰＲｉ^ｓ _ｒ（ｔ）、ＰＲｊ^ｓ _ｒ（ｔ）を先の数７に代入して合成値ＩＦＣＰｉｊｒ^ｓ _ｒ（ｔ）及びＩＦＰＡｉｊｒ^ｓ _ｒ（ｔ）を求め、合成値記憶部８３に記憶される。

前記状態量推定部８４は前記合成値に基づいて例えばカルマンフィルタにより状態量を求めるものであり、その入力値と出力値とは図３（ａ）に示す通りである。

ここでいう状態量とは、数７を変形して得られた既述の数８の右辺の項目である。この実施の形態では、当該右辺の項目のうち搬送波位相観測ノイズｅ^ｓ _ｒ（ｔ）及び擬似距離観測ノイズＥ^ｓ _ｒについては平均値をゼロと仮定することによりゼロとして取り扱っている。このためＮ^ｓ _ｒは数８の２段目の式から第２の合成値であるＩＦＰＡｉｊｒ^ｓ _ｒ（ｔ）の値になる。また、光速の対流圏遅延に基づく誤差であるＴ^ｓ _ｒ（ｔ）は地上のモニタ局７に送られる測位信号の入射角に応じて変わってくるので、例えばモニタ局７に対する測位衛星６の相対位置に応じて決まる値として取り扱うことができ、観測の時間帯毎に割り当てられた定数として扱うようにしてもよい。このようなことから、演算サイクルで推定する状態量とは、図３（ｂ）に示すように測位衛星６の位置（例えば、慣性座標系におけるＸ、Ｙ、Ｚ方向の座標）、測位衛星６の速度（Ｘ、Ｙ、Ｚ方向の各速度）、既述のモニタ局ｒの受信機時計誤差＋周波数間バイアス（ｆ_ｉ、ｆ_ｊの測位信号を検出するタイミングの僅かなずれ）に起因する観測データの誤差（Ｄｔ_ｓ（ｔ））、測位衛星ｓの搭載時計誤差＋周波数間バイアス（ｆ_ｉ、ｆ_ｊの測位信号を発信するタイミングの僅かなずれ）に起因する観測データの誤差Ｄｔ^ｓ（ｔ）である。測位衛星６Ａ〜６Ｃの位置については、各モニタ局７Ａ〜７Ｃ及び各測位衛星６Ａ〜６Ｃの組み合わせ毎の観測データが主管制局８に取り込まれているので、各モニタ局７Ａ〜７Ｃ及び各測位衛星６Ａ〜６Ｃ間の距離の推定値が分り、これらの推定値から各測位衛星６Ａ〜６Ｃの位置が推定できる。また、各測位衛星６Ａ〜６Ｃの速度は、モニタ局７における観測のサンプリングの時間間隔と、相前後するサンプリングにおける測位衛星６の各位置とに基づいて推定することができる。

そして、例えば、モニタ局７において時刻ｔ_ｋにて擬似距離及び搬送波位相を観測したとすると、この観測データに基づいて演算された時刻ｔ_ｋにおける合成値ＩＦＣＰｉｊｒ^ｓ _ｒ（ｔ）及びＩＦＰＡｉｊｒ^ｓ _ｒ（ｔ）が状態量推定部８４に入力されるが、これら合成値の分散も入力される。この分散とは受信機等の性能に基づいて決められる合成値の正確性を示す値であり、この例では定数として取り扱っている。

更に前記状態量推定部８４には、当該演算サイクルよりも前の演算サイクル、例えば１つ前の演算サイクルで推定した状態量とこれら状態量の推定値の各々に対する共分散（推定値の正確性を示す値）が入力される。一つ前の演算サイクルで推定した状態量とは、地上のモニタ局７において時刻ｔ_ｋよりも一つ前のサンプリング時刻ｔ_ｋ−１において観測した観測データに基づいて合成値を求め、その合成値と更に前の演算サイクルで求めた状態量推定値等に基づいて推定した状態量、つまり時刻ｔ_ｋ−１における状態量推定値である。

こうして今回（時刻ｔ_ｋ）の合成値及びその分散と、前回（ｔ_ｋ−１）の状態量の推定値及び各推定値の共分散とが状態量推定部８４に入力され、カルマンフィルタによる演算が行われて、今回（時刻ｔ_ｋ）の状態量の推定値と各推定値の共分散とが求められて推定状態量記憶部８５に記憶される。

ここでカルマンフィルタによる演算について簡単に述べておくと、前回（時刻ｔ_ｋ−１）の演算サイクルで推定した測位衛星６の位置と速度とに基づいて今回（時刻ｔ_ｋ）の演算サイクル時における測位衛星６の位置を推定した測位衛星６の位置を推定し、また前回（時刻ｔ_ｋ−１）のＤｔ_ｒ（ｔ）と所定のモデルとに基づいて、今回（時刻ｔ_ｋ）のＤｔ^ｓ（ｔ）を推定する。上述の所定のモデルとは、周波数間バイアスは一定であって定数であることから、定数＋時計誤差のモデルで決まる値になるが、この時計誤差のモデルは時計の種別に応じて時間の関数として表される。そして前回（時刻ｔ_ｋ−１）の状態量推定値に基づいて推定された今回（時刻ｔ_ｋ）の各状態量を数８の上段の式の右辺に代入して、今回（時刻ｔ_ｋ）の第１の合成値ＩＦＣＰｉｊｒ^ｓ _ｒ（ｔ）を求め、得られた第１の合成値と観測データを数７に代入して得た第１の合成値とを比較してその差分を求めると共に、第２の合成値ＩＦＰＡｉｊｒ^ｓ _ｒ（ｔ）についても同様にして差分を求め、これら差分に応じて今回（時刻ｔ_ｋ）の各状態量推定値を修正する。このよにして修正された各状態量が推定状態量記憶部８５に記憶され、次の演算サイクルにて使用されると共に、航法メッサージ作成部８６により航法メッセージの中に取り込まれ、この航法メッセージは送信機８７によって各測位衛星６Ａ〜６Ｃに送信される。

以上の説明から分るように、各モニタ局７Ａ〜７Ｃにおける各測位衛星６Ａ〜６Ｃから観測データに基づいて、そのときの各測位衛星６Ａ〜６Ｃの位置、速度を含む既述の状態が推定できる。

上述の実施の形態によれば、擬似距離及び搬送波位相の観測データを用いて合成値を得るにあたって、擬似距離及び搬送波位相のいずれについても２つの観測データ（周波数ｆ_ｉ、ｆ_ｊに対応する観測データ）を用いているので、ＧＰＳで使用されている二重差の方式（数５）に比べて観測ノイズの分散が４分の１となるので、各状態値を高精度に推定できる。従って測位衛星の軌道を高精度に推定でき、測位衛星を活用したアプリケーション、例えばナビゲーションシステムを信頼性の高いものとすることができる。

日本で計画中の準天頂衛星及びＧＰＳ衛星に本発明の測位衛星の位置検出システムを適用し、準天頂衛星及びＧＰＳ衛星の軌道をシュミレーションにより推定した結果例を図４及び図５に示す。図４の（ａ）、（ｂ）及び（ｃ）は、準天頂衛星の軌道推定において，（ａ）はＸ座標位置の推定誤差、（ｂ）はＹ座標位置の推定誤差、（ｃ）はＺ座標位置の推定誤差を表したものである。ここでいう推定誤差とは、既述のように前回（ｔ_ｋ−１）の状態値に基づいて今回（ｔ_ｋ）推定した状態値と、数８の左辺と右辺との誤差に基づいて状態値を修正して得られた今回（ｔ_ｋ）の状態値との差分である。縦軸はその座標での推定誤差（ｃｍ）であり、横軸は観測時間（ｈ）である。図４の（ａ）、（ｂ）及び（ｃ）において、従来の測位衛星観測システムでは推定誤差が５ｍなのに対して、測定時間を通して推定誤差が０に近い値を示すことから、準天頂衛星の軌道を高精度に推定することができる。

また、図５は、ＧＰＳ衛星の同様の誤差の推定を軌道推定シュミレーションにより求めたものである。図５の（ａ）はＸＹＺ位置推定誤差であり、縦軸は推定誤差（ｃｍ）、横軸は観測時間（ｈ）である。図５の（ｂ）は時計誤差によるバイアス推定誤差であり、縦軸は推定誤差（ｃｍ）、横軸は観測時間（ｈ）である。図５の（ｃ）は地上のモニタ局から見たモニタ局とＧＰＳ衛星との測定角度であり、縦軸は仰角（℃）、横軸は観測時間（ｈ）である。図５の（ａ）及び（ｂ）の各推定誤差（ｃｍ）は図５の（ｃ）に基づいて行っている。よってどの測定角度においても、（ａ）及び（ｂ）は常に０に近い測定誤差であることがわかる。

本発明に係る人工衛星の位置推定システムの実施の構成を示す構成図である。 本発明に係る測位衛星から送られてくる測位データのパルス列を示すイメージ図である。 合成演算部及び状態量推定部の入力、出力の関係と、状態量とを説明する説明図である。 本発明のシステムを用いて準天頂衛星位置の推定誤差を示す推定図である。 本発明のシステムを用いてＧＰＳ衛星位置の推定誤差を示す推定図である。 従来の観測方法において測位衛星とモニタ局とを模式的に示す説明図である。 従来の観測方法において搬送波位相による距離計測の原理を説明した説明図である。 従来の観測方法において二重差による観測の様子を示す説明図である。

符号の説明

６ 測位衛星
６１ 搭載時計
６２ 受信機
６３ 測位信号生成部
６４ 送信機
７ モニタ局
７１ 受信機時計
７２ 受信機
８ 主管制局
８１ 観測データ記憶部
８２ 合成値演算部
８３ 合成値記憶部
８４ 状態量推定部
８５ 推定状態量記憶部
８６ 航法メッセージ
８７ 送信機
１００ 一般ユーザ
１０１ 受信機

Claims (3)

1. 人工衛星から送信され、その先頭からの位置に対応する位置コードが配列された所定長さのパルス列を搬送波に載せた測位信号を地上のモニタ局にて受信し、モニタ局にて検出した位置コードと光速とに基づいて人工衛星とモニタ局との擬似距離観測データを求めると共に、
   モニタ局にて検出した前記搬送波の波数を所定の時刻ｔ₀から前記位置コードをサンプリングした時刻ｔ_kまでの時間カウントし、そのカウント値と前記搬送波の周波数と前記カウントした時間とに基づいて、時刻ｔ₀におけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値と時刻ｔ_kにおけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値との差分である搬送波位相観測データを求め、
   前記擬似距離観測データと前記搬送波位相観測データとに基づいて、人工衛星とモニタ局との距離を求め、当該人工衛星と複数のモニタ局の各々との距離に基づいて当該人工衛星の位置を求めるシステムにおいて、
   イ． 人工衛星から送信される周波数ｆ_iの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｉ^s _r(t)及び搬送波位相観測データＣＰｉ^s _r(t)と、
   前記測位信号と同時に送られる、周波数ｆ_jの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｊ^s _r(t)及び搬送波位相観測データＣＰｊ^s _r(t)と、
   を下記数１に代入して第１の合成値ＩＦＣＰｉｊ^s _r(t)及び第２の合成値ＩＦＰＡｉｊ^s _r(t)を算出する合成値演算部と、
   ロ． 時刻ｔ_kに観測した観測データに基づいて計算された第１の合成値ＩＦＣＰｉｊ^s _r(t)及び第２の合成値ＩＦＰＡｉｊ^s _r(t)と、
   数１を変形して得られた下記数２と、
   時刻ｔ_kのサンプリングよりも一つ前のサンプリングのタイミングｔ_k−1にて観測した観測データを用いて計算された第１の合成値ＩＦＣＰｉｊ^s _r(t)及び第２の合成値ＩＦＰＡｉｊ^s _r(t)に基づいて推定した人工衛星の位置と、
   に基づいて、時刻ｔ_kにおける人工衛星の位置を推定する推定部と、
   を備えたことを特徴とする人工衛星の位置推定システム。
   記
   

   

   ただし、λ_i及びλ_jは夫々周波数ｆ_iの搬送波の波長及び周波数ｆ_jの搬送波の波長である。
   

   

   ただし、ｐ^s _r (t)は人工衛星及びモニタ局の距離、Ｔ^s _r(t)は電波の対流圏遅延に起因する誤差、Ｄｔ_r(t)はモニタ局の時計誤差と当該モニタ局の受信機における周波数ｆ_iの搬送波及び周波数ｆ_jの搬送波の受信のタイミングの誤差との加算値に起因する誤差、Ｄｔ^s(t)は人工衛星の搭載時計の誤差と当該人工衛星における周波数ｆ_iの搬送波及び周波数ｆ_jの搬送波の発信のタイミングの誤差との加算値に起因する誤差、Ｎ^s _r (t)は搬送波位相アンビギュイティ、ｅ^s _r (t)は搬送波位相観測ノイズ、Ｅ^s _r (t)は擬似距離観測ノイズである。
2. 前記推定部は、時刻ｔ_kのサンプリングよりも一つ前のサンプリングのタイミングｔ_k−1にて観測した観測データを用いて計算された第１の合成値ＩＦＣＰｉｊ^s _r(t)及び第２の合成値ＩＦＰＡｉｊ^s _r(t)に基づいて推定した前記Ｄｔ_r(t)及びＤｔ^s(t)を更に用いて、時刻ｔ_kにおける人工衛星の位置を推定することを特徴とする請求項１記載の人工衛星の位置検出システム。
3. 前記推定部は、第１の合成値ＩＦＣＰｉｊ^s _r(t)及び第２の合成値ＩＦＰＡｉｊ^s _r(t)の分散と、人工衛星の位置を含む状態量の各分散と、を更に用いて、時刻ｔ_kにおける人工衛星の位置を推定することを特徴とする請求項１記載の人工衛星の位置推定システム。
   

Images