(57)【要約】【課題】  人工衛星の軌道を推定するにあたり、観測ノイズの分散を小さくし、高精度な位置推定を行うこと。【解決手段】  人工衛星から周波数ｆi（波長λi）の搬送波で送信される測位信号を観測して得た擬似距離観測データＰＲｉs r(t)及び搬送波位相観測データＣＰｉs r(t)と、前記測位信号と同時に周波数ｆj（波長λj）の搬送波で送られる測位信号とを観測して得た擬似距離観測データＰＲｊs r(t)及び搬送波位相観測データＣＰｊs r(t)と、を下記数１に代入して第１の合成値及び第２の合成値を算出し、これら合成値に基づいて人工衛星の位置を推定する。【選択図】  図１
2005062046.tif 000002
【特許請求の範囲】
【請求項１】
人工衛星から送信され、その先頭からの位置に対応する位置コードが配列された所定長さのパルス列を搬送波に載せた測位信号を地上のモニタ局にて受信し、モニタ局にて検出した位置コードと光速とに基づいて人工衛星とモニタ局との擬似距離観測データを求めると共に、
モニタ局にて検出した前記搬送波の波数を所定の時刻ｔ0から前記位置コードをサンプリングした時刻ｔkまでの時間カウントし、そのカウント値と前記搬送波の周波数と前記カウントした時間とに基づいて、時刻ｔ0におけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値と時刻ｔ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における人工衛星の位置を推定する推定部と、
を備えたことを特徴とする人工衛星の位置推定システム。
記
【数１】
2005062046.tif 000013
ただし、λi及びλjは夫々周波数ｆiの搬送波の波長及び周波数ｆjの搬送波の波長である。
【数２】
2005062046.tif 000014
ただし、ｐs r  (t)は人工衛星及びモニタ局の距離、Ｔs r(t)は電波の対流圏遅延に起因する誤差、Ｄｔr(t)はモニタ局の時計誤差と当該モニタ局の受信機における周波数ｆiの搬送波及び周波数ｆjの搬送波の受信のタイミングの誤差との加算値に起因する誤差、Ｄｔs(t)は人工衛星の搭載時計の誤差と当該人工衛星における周波数ｆiの搬送波及び周波数ｆjの搬送波の発信のタイミングの誤差との加算値に起因する誤差、Ｎs r  (t)は搬送波位相アンビギュイティ、ｅs r  (t)は搬送波位相観測ノイズ、Ｅs r  (t)は擬似距離観測ノイズである。
【請求項２】
前記推定部は、時刻ｔkのサンプリングよりも一つ前のサンプリングのタイミングｔk−1にて観測した観測データを用いて計算された第１の合成値ＩＦＣＰｉｊs r(t)及び第２の合成値ＩＦＰＡｉｊs r(t)に基づいて推定した前記Ｄｔr(t)及びＤｔs(t)を更に用いて、時刻ｔkにおける人工衛星の位置を推定することを特徴とする請求項１記載の人工衛星の位置検出システム。
【請求項３】
前記推定部は、第１の合成値ＩＦＣＰｉｊs r(t)及び第２の合成値ＩＦＰＡｉｊs r(t)の分散と、人工衛星の位置を含む状態量の各分散と、を更に用いて、時刻ｔkにおける人工衛星の位置を推定することを特徴とする請求項１記載の人工衛星の位置推定システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、人工衛星から送信される測位信号の観測データから測位衛星等の人工衛星の位置を推定する人工衛星の位置推定システムに関する。
【背景技術】
【０００２】
米国のＧＰＳ（Ｇｌｏｂａｌ  Ｐｏｓｉｔｉｏｎｉｎｇ  Ｓｙｓｔｅｍ）、ロシアのＧＬＯＮＡＳＳ、ヨーロッパのＧＡＬＩＲＥＯ、更に日本で計画中の準天頂衛星システム等の測位衛星システムを利用して精密な測定位置及び時刻同期を行うには、測位衛星が地球の周りを周回する軌道及び測位衛星に搭載した時計の誤差を高精度に推定することが必要不可欠である。
【０００３】
測位衛星から送信される測位信号を利用して測位衛星の軌道及び測位衛星に搭載した時計誤差を推定する際に必要となるモデルは、測位衛星の運動を表す運動（力学）モデルと測位信号を受信した観測データに含まれる各種の誤差を表現する観測モデルとが必要である。この測位信号は測位衛星と観測点間との距離を測るために使用されるが、測位信号の観測データの中には測位衛星と観測点間との距離データの他に測位衛星に搭載された時計の時間誤差、電離層や対流圏による電波遅延、受信機ノイズなど様々なデータ誤差が含まれる。これらのデータ誤差を補正及び除去することが高精度の測位衛星の軌道及び測位衛星に搭載した時計誤差の推定では必要となる。
【０００４】
ここで図６を用いて測位衛星の軌道を推定する原理について説明する。地球の大気圏外領域を回っている測位衛星１からは測位信号Ｚと呼ばれる一定の信号形式を持つ電波が常時、測位衛星１に搭載された衛星時計１１に同期して送信機１２によって送信される。この測位信号Ｚは、地上のモニタ局２Ａ〜２Ｃに設置された受信機３Ａ〜３Ｃによって受信され、以下の原理で測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を測定する。
【０００５】
先ず、測位衛星１が測位信号Ｚを送信した時刻は、各モニタ局２Ａ〜２Ｃで受信した信号を解読することにより知ることができる（送信時刻（ｔ））。モニタ局２Ａ〜２Ｃで測位信号Ｚを受信した時刻はモニタ局２Ａ〜２Ｃに設置されたモニタ局時計４Ａ〜４Ｃで知ることができる（受信時刻（Ｔ））。このことから測位信号Ｚが測位衛星１からモニタ局２Ａ〜２Ｃまでの空間を伝播した時間（伝播時間（Δｔ））はΔｔ＝Ｔ−ｔとなり、よって測位信号Ｚは光速ｃで空間を伝わることから測位衛星１とモニタ局２Ａ〜２Ｃ間との距離ＤはＤ＝Δｔ×ｃで表すことができる。
【０００６】
このように予め正確な位置が分っている複数のモニタ局２Ａ〜２Ｃからモニタ局２Ａ〜２Ｃと測位衛星１との各距離データＤ１〜Ｄ３を主管制局５に設けられた演算部５１で演算することにより測位衛星１の位置を推定することができる。このシステムにより得られた測位衛星１の位置は、例えば地上の移動体のナビゲーションシステムにも適用される。
【０００７】
上述した測位衛星１の運動を表す運動（力学）モデルは、地球の引力等の物理法則に支配されており、ある時刻での位置、速度が分れば一定時間後の位置、速度を予測することができるので、測位衛星１の軌道推定では、測位信号Ｚを正確に測定することが重要となる。
【０００８】
そこで、測位信号Ｚを利用した測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を測定する方法には、測位信号Ｚの信号形式（一定の長さを持ったパターン）を利用した擬似距離を用いる方法と、測位信号Ｚの搬送波を数える搬送波位相を用いる方法とがある。
【０００９】
この擬似距離を用いる方法では、測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を直接測ることができるが、計測精度が信号パターンを構成するデータの単位長さに制限されるため距離の誤差が大きくなる。一方前記搬送位相を用いる方法では、計測精度は測位信号Ｚの波長以下となり測定精度を上げることができるが、受信電波の波を数えるため波長の整数倍の不確定性（アンビギュイティ）が残ってしまう。
【００１０】
また、前記測位衛星１からは、詳細には、異なった周波数で複数の測位信号Ｚが同時に送信されており、これら複数の測位信号Ｚにおいて、各々の測位信号Ｚ毎に擬似距離及び搬送波位相が求められる。しかしながら、測位衛星１から送信された測位信号Ｚを地上のモニタ局２Ａ〜２Ｃに設置された受信機３Ａ〜３Ｃで受信すると測位衛星１と地上のモニタ局２Ａ〜２Ｃ間との距離データＤ１〜Ｄ３の中には（表１）に示す各種誤差が混入する。
【００１１】
【表１】
2005062046.tif 000003
これらのモニタ局２Ａ〜２Ｃで観測された擬似距離観測データと、搬送波位相観測データと、に含まれる観測誤差を表した観測モデルは、一般的に数３のように表される。
【００１２】
【数３】
2005062046.tif 000004
ここで、上述の数３を用いて測位衛星１の軌道決定のための推測処理をどのようにして行うのかを説明する。数３の式の左辺は観測値であり、例えば数３中のＰＲｉｓｒ（ｔ）は擬似距離観測値である。この擬似距離観測値は次のようにして求められる。測位衛星１は、例えば地上から２万ｋｍの高度の所を移動しており、この測位衛星１からパルスパターンを含む例えば長さ３００〜５００ｋｍ程度のパルス列を地上に送る。このパルスパターンは、パルス列の先頭からの距離に対応したパターンとして配列され、例えばパルス列の先頭からＡｃｍの位置にあればＡｃｍを表すパルスパターンが形成される。つまり、このパルス列は測位衛星１からいわば長さ３００〜５００ｋｍの物差しが飛来してくるのと同じである。
【００１３】
一方測位衛星１と地上のモニタ局２Ａ〜２Ｃとにおいては、正確に時計合わせされたモニタ局時計４Ａ〜４Ｃが夫々設けられ、測位衛星１からパルス列を発信するタイミング及びモニタ局２Ａ〜２Ｃでこのパルス列を検出するタイミングが予め決められている。従ってモニタ局２Ａ〜２Ｃの検出のタイミングでパルス列のあるパルスパターンを把握すれば、その物差しのいくつの値をどのタイミングで検出したかが分る。測位衛星１とモニタ局２Ａ〜２Ｃとの距離は運動（力学）モデルにより概ね分っているので（例えば数ｍの誤差で分っているので）、この測位衛星１とモニタ局２Ａ〜２Ｃとの距離に比べてパルス列の長さが短くても、当該パルス列の中のパルスパターンをある時刻で検出すれば、そのパルスパターンはパルス列を物差しと見立てたときのスケール値に相当するので、当該パルス列の先頭からどれだけの時間だけ遅れて送信されているかが分り、結局パルス列の先頭をモニタ局２Ａ〜２Ｃが受信した時刻が分る。この結果、パルス列の先頭が飛来してくる時間（（モニタ局２Ａ〜２Ｃで受信した時刻）−（測位衛星１から送信された時刻））が検出され、この時間に光速ｃを掛けることで測位衛星１とモニタ局２Ａ〜２Ｃとの距離を算出することができる。しかしながら、実際には数３から分るようにこの距離は、種々の誤差を含んでいる。
【００１４】
また上記の数３中のＣＰｉｓｒ（ｔ）を用いても同様にして測位衛星１とモニタ局２Ａ〜２Ｃとの距離ρｓｒ（ｔ）を用いることができる。この場合には、擬似距離の観測データを用いる代わりに搬送波の位相を観測した観測データを用いている。この搬送波とは、測位衛星１から送られてくる既述のパルス列の搬送波であり、搬送波位相観測値ＣＰｉｓｒ（ｔ）とは、ある時刻ｔ０からある時刻ｔまでに地上のモニタ局２Ａ〜２Ｃで搬送波をカウントし、搬送波をカウントした値から時刻ｔ０における測位衛星１の位置と時刻ｔにおける測位衛星１の位置とが変わることに基づいて搬送波の波数値を求めたものである。この搬送波位相による距離計測の原理について図７と数４を参照しながら詳しく説明する。
【００１５】
先ず図７で、モニタ局２に設けられた受信機３で搬送波を受信開始したモニタ局時計４の時刻をｔ０とする。この時点で受信した搬送波が測位衛星１に設けられた送信機１２から送信されたときの衛星時計１の時刻は、測位衛星１とモニタ局２間との距離をｒ（ｔ０）とするとｔ０−｛ｒ（ｔ０）／ｃ｝となり、この時点では受信機３で受信した受信波数は０個である。この時点から同時に受信機３内部の図示されない局部発振機の波数積算も開始される。ここで搬送波の計測時刻をｔとする。この時点で受信機３で受信した搬送波が測位衛星１に設けられた送信機１２から送信されたときの衛星時計１１の時刻は、測位衛星１とモニタ局２間との距離をｒ（ｔ）とするとｔ−｛ｒ（ｔ）／ｃ｝となり、受信開始時刻から計測時刻までに受信した波数Φ（ｔ）は数４で表される。
【００１６】
【数４】
2005062046.tif 000005
一般に静止系であれば、所定時間に受信する電波の数は一定であり、数４の右辺の前段の項である１／２πω・（ｔ−ｔ０）だけで決定されるが、測位衛星１は移動しているため、後段の項である{ｒ（ｔ）／λ}−{ｒ（ｔ０）／λ}が含まれてくる。この値が搬送波の波数値であり、ＣＰｉｓｒ（ｔ）に相当するものである。しかしながら、上述したようにこの搬送波の位相を利用する方法では、計測精度は測位信号Ｚの波長以下となり高精度であるが、受信電波の波数を数えるため波長の整数倍の不確定性（アンビギュイティ）が残ってしまう。
【００１７】
測位衛星１の軌道推定は、測位衛星１と地上のモニタ局２Ａ〜２Ｃ間との距離データＤ１〜Ｄ３を利用して測位衛星１の位置、速度を求めることであるが、測位信号Ｚの観測データ中には測位衛星１とモニタ局２Ａ〜２Ｃ間との距離以外の誤差が含まれているため、これらの誤差の値も同時に推定する必要が生じてくる。しかし、前記数３で表した観測モデルのままでは未知数が多く夫々の観測誤差を分離して推定することができない。そこで擬似距離観測データと搬送波位相観測データとを次のように合成した擬似的な観測データを利用して未知数を減らす方法が一般的な手法である（例えば、非特許文献１参照。）。この一般的な手法を図８と数５を参照しながら詳しく説明する。
【００１８】
【数５】
2005062046.tif 000006
この方法は、図８に示すように２つのモニタ局ｒ（ｒ´）と二つの測位衛星ｓ（ｓ´）とを用い、更に各測位衛星ｓ（ｓ´）のいずれからも送信される周波数ｆｉの搬送波と周波数ｆｊの搬送波とを利用する方法（二重差観測データ）であり、先の数３で表される擬似距離観測値と搬送波位相の波数値とを数５のように定義する。これらの式においてＣＰｉｓｒ（ｔ）は、モニタ局ｒにおいて、測位衛星ｓから送信される周波数ｆｉの搬送波について数４のように時刻ｔ〜ｔ０間で搬送波をカウントし、数４中のωは既値であることから、後段の項を求めた値である。また、ＣＰｊｓｒ´（ｔ）は、モニタ局ｒ´において測位衛星ｓから送信される周波数ｆｊの搬送波について求めた値である。更にＰＲｉｓｒ（ｔ）は、測位衛星ｓから周波数ｆｉの搬送波に乗せられて送られる既述のパルス列中のパルスパターンをモニタ局ｒについて検出し、その検出結果に基づいて得られた擬似距離観測値である。
【００１９】
そして、数５の右辺の各項に、数３で定義される値（右辺）を代入することにより測位衛星ｓ（ｓ´）と地上のモニタ局ｒ（ｒ´）間との距離以外の誤差要因（観測誤差）は時計誤差、電離層遅延及び周波数間バイアスが合成により打消されて、対流圏遅延、搬送波位相アンビギュイティ及び観測ノイズのみが残ることになる（数６）。
【００２０】
【数６】
2005062046.tif 000007
もう少し詳しく説明すると、数５の右辺は、各モニタ局ｒ（ｒ´）における観測値を演算した値であるから、各モニタ局ｒ（ｒ´）における観測値{ＰＲｉｓｒ（ｔ）等}を地上回線で主管制局のコンピュータに取り込むことによって、数５の左辺であるＩＣＰｉｊｒｓｓ´ｒｒ´（ｔ）及びＩＡＭｉｊｒｓｓ´ｒｒ´（ｔ）が求まる。
【００２１】
今、時刻ｔｋにおける擬似距離観測値と搬送波の波数値との各観測値を用いて数５を演算した結果を夫々ＩＣＰｉｊｒｓｓ´ｒｒ´（ｔｋ）及びＩＡＭｉｊｒｓｓ´ｒｒ´（ｔｋ）とする。なお、時刻ｔｋにおける搬送波の波数値は、時刻ｔ０からの搬送波の積算値に基づいて求められるが、ω（ｔｋ−ｔ０）の値は既にわかっているので、この値をキャンセルすることにより、結局、数４の{ｒ（ｔｋ）／λ}−{ｒ（ｔ０）／λ}のみが波数値として求まることになる。
【００２２】
一方コンピュータは、時刻ｔｋよりも前の時刻ｔｋ−１において推定した数６の右辺の測位衛星ｓ（ｓ´）とモニタ局ｒ（ｒ´）間との距離及び推定状態量に基づいて時刻ｔｋの右辺の値を求める。ここで測位衛星ｓ（ｓ´）とモニタ局ｒ（ｒ´）間との距離であるρｓｒ（ｔ）については初期値としては適当な値を用いる。但し、測位衛星ｓ（ｓ´）の軌道は予めかなりの精度で把握されているので、適当な値とはいっても高い精度を持つ値である。また、他の推定状態量のうちＴｓｒ（ｔ）は信頼できる理論的モデル（時間の関数）を用い、時刻ｔｋ−１のときの対流圏遅延を用いて、時刻ｔｋのときの対流圏遅延（Ｔｓｒ（ｔ））を求めることができる。またＮｉｓｒは連続して電波受信中は定数として考えることができ、更にＥｉｓｒ（ｔ）、ｅｉｓｒ（ｔ）はいずれも平均値は０と仮定して取り扱うことができる。こうして数６を用いて、数５の右辺が求まる。先に求めたＩＣＰｉｊｓｓ´ｒｒ´（ｔ）及びＩＡＭｉｊｓｓ´ｒｒ´（ｔ）は数５の左辺であるため、本来は左辺の観測値と右辺の演算値とは一致するはずであるが、実際にはモデルの誤差分や種々の仮定が含まれているため、また時刻ｔｋの値があくまでも推定値であるため、観測値と演算値とは異なってくる。
【００２３】
そこで、コンピュータは数５の夫々について左辺と右辺との差（偏差）を求める。そして、これら偏差に応じて時刻ｔにおけるρｓｒ（ｔ）、Ｔｓｒ（ｔ）及びＮｉｓｒ（ｔ）の値を修正する。実際には、主管制局では各モニタ局毎でこれらの推定値が求まり、時刻ｔｋにおける推定値としてメモリ内に格納され、次の時刻ｔｋ＋１における推定に用いられる。また、推定初期値においては、誤差が大きい場合でも時間が経過するにつれて精度が高くなり前記誤差は小さくなっていく。時刻ｔｋにおける測位衛星の位置、速度は３ヶ所以上のモニタ局と測位衛星との距離から求まることになる。
【００２４】
【非特許文献１】Global Positioning System:Theory and Applications Volume1 Part2(American Institute of Aernautics and Astronautics 1996)
【発明の開示】
【発明が解決しようとする課題】
【００２５】
しかしながら、擬似距離観測データと搬送波位相観測データとを上述のように合成して未知数を減らす一般的な手法（二重差観測データ）では、二重差の使用により測位衛星の軌道は決定可能となる８つの観測データを合成することにより観測ノイズの分散は１つの観測データの場合に比べて８倍となり精密な軌道決定には不利である。また近年、この二重差の不利を補うためにＧＰＳの軌道決定において、一般的に全世界の観測網を利用して長時間連続的に測位衛星の観測及び軌道推定を行うことで除々に高精度な軌道に収束させる方法が用いられているが、地上に多数のモニタ局が必要となるため、多大なコストが掛かるという問題がある。
【００２６】
本発明は係る事情に鑑みてなされたものであって、本発明の目的は、二重差を使用せずに、測位衛星から送信される測位信号に含まれる観測誤差（観測ノイズ）を小さく抑え、測位衛星等の人工衛星の位置の推定を可能とする人工衛星の位置推定システムを提供することにある。
【課題を解決するための手段】
【００２７】
本発明は、人工衛星から送信され、その先頭からの位置に対応する位置コードが配列された所定の長さのパルス列を搬送波に載せた測位信号を地上のモニタ局にて受信し、モニタ局にて検出した位置コードと光速とに基づいて人工衛星とモニタ局との擬似距離観測データを求めると共に、
モニタ局にて検出した前記搬送波の波数を所定の時刻ｔ0から前記位置コードをサンプリングした時刻ｔkまでの時間をカウントし、そのカウント値と前記搬送波の周波数と前記カウントした時間とに基づいて、時刻ｔ0におけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値と時刻ｔkにおけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値との差分である搬送波位相観測データを求め、
前記擬似距離観測データと前記搬送波位相観測データとに基づいて、人工衛星とモニタ局との距離を求め、当該人工衛星と複数のモニタ局の各々との距離に基づいて当該人工衛星の位置を求めるシステムにおいて、
イ．  人工衛星から送信される周波数ｆiの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｉs r(t)及び搬送波位相観測データＣＰｉs r(t)と、
前記測位信号と同時に送られる、周波数ｆjの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｊs r(t)及び搬送波位相観測データＣＰｊs r(t)と、
を下記数７に代入して第１の合成値及び第２の合成値を算出する合成値演算部と、
ロ．  時刻ｔkに観測した観測データに基づいて計算された第１の合成値及び第２の合成値と、
数７を変形して得られた下記数８と、
時刻ｔkのサンプリングよりも一つ前のサンプリングのタイミングｔk-1にて観測した観測データを用いて計算された第１の合成値及び第２の合成値に基づいて推定した人工衛星の位置と、
に基づいて、時刻ｔkにおける人工衛星の位置を推定する推定部と、を備えたことを特徴とする人工衛星の位置推定システムである。
【００２８】
記
【００２９】
【数７】
2005062046.tif 000008
ただし、λi及びλjは夫々周波数ｆiの搬送波の波長及び周波数ｆjの搬送波の波長である。
【００３０】
【数８】
2005062046.tif 000009
ただし、ρ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における人工衛星の位置を推定することが好ましい。
【００３２】
ここで数７と数８との関係について説明しておく。数７の右辺の各項に数３の右辺を代入することにより数９が得られる。
【００３３】
【数９】
2005062046.tif 000010
数９の｛  ｝内を数１０のように合成された１つの未知数として取り扱うと数８が得られる。
【００３４】
【数１０】
2005062046.tif 000011
数８において擬似距離観測ノイズＥ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に観測した観測データに基づいて計算された第１の合成値及び第２の合成値とこれら合成値の分散を用いて、今回の演算サイクルにおける状態量及び分散を推定する。なお第１の合成値及び第２の合成値の分散は定数として扱うことができる。
【００３６】
また前記推定部は、人工衛星とモニタ局との組毎に各状態量を推定するので、推定演算の過程では例えば人工衛星とモニタ局との距離を推定し、各人工衛星の位置についてはその人工衛星と３つ以上のモニタ局との距離から計算することで推定している。こうして各人工衛星の位置を推定することができ、例えば推定された位置と一つ前の演算サイクルで求めた位置とから人工衛星の速度も推定できる。
【発明の効果】
【００３７】
ＧＰＳで行われている二重差による手法（数５の手法）では、擬似距離及び搬送波位相の各々について８つの観測データを合成していたが、本発明では２つの観測データを合成した合成値に基づいて測位衛星等の人工衛星の位置が推定でき、従って観測ノイズの分散は二重差の手法に比べて４分の１となり、高精度な位置推定を行うことができる。更に測位衛星等の人工衛星の時計誤差及びモニタ局の時計誤差の推定も高精度に行うことができる。
【発明を実施するための最良の形態】
【００３８】
本発明の人工衛星の位置推定システムについて、図１〜図３に基づいて説明を行う。図１は本発明の実施の形態に係る人工衛星である測位衛星の位置測定システムを含む測位衛星の軌道推定システムの説明図である。図１中６は測位衛星であり、６Ａ、６Ｂ、６Ｃは３機の測位衛星の夫々を指している。前記測位衛星６は搭載時計６１、受信機６２、測位信号生成部６３及び送信機６４を備えている。前記測位信号生成部６３は、搭載時計６１からクロック信号に基づいて、周波数ｆｉ（例えば１５０ＧＨｚ）の搬送波で送られる測位信号であるパルス列と周波数ｆｊ（例えば１７０ＧＨｚ）の搬送波で送られる測位信号であるパルス列とを作成して前記送信機６４に受け渡す機能を備えている。これらパルス列は例えば３００ｋｍ〜５００ｋｍ程度の長さに作られており、パルス列からの先頭からの位置に対応するコードが配列されている。図２はパルス列を示すイメージ図であり、先頭には当該パルス列を発信する測位衛星６の識別コード（ＩＤ）が配置されると共に、その後には各位置に対応したコード（Ｐ１、Ｐ２、…Ｐｎ）が配置されている。各コードの長さは例えば２０ｃｍ程度とされる。
【００３９】
なお、測位衛星６は後述の主管制局８から、Ｘ、Ｙ、Ｚ座標における自己の正確な位置（Ｘ、Ｙ、Ｚ）と速度（Ｖｘ、Ｖｙ、Ｖｚ）とを含む航法メッセージを受け取り、その航法メッセージを測位信号と共に一般ユーザ１００の受信機１０１に配信する機能を備えている。前記一般ユーザ１００は、例えば車両等の移動体であり、３機の測位衛星６Ａ〜６Ｃから送られてくる航法メッセージと測位信号とに基づいて、各測位衛星６Ａ〜６Ｃまでの各距離が分るので、当該一般ユーザ１００（移動体）の位置を知ることができる。
【００４０】
また、図１中７は地上に設けられたモニタ局であり、７Ａ〜７Ｃは３個のモニタ局を夫々示している。モニタ局７は受信機時計７１及び受信機７２を備えており、受信機時計７１から出力されるクロック信号に基づいて測位衛星６から送られてくる測位信号（パルス列及び搬送波）を検出し、周波数ｆｉの測位信号に基づいて擬似距離観測データＰＲｉｓｒ（ｔ）及び搬送波位相観測データＣＰｉｓｒ（ｔ）を求め、また、周波数ｆｊの測位信号に基づいて擬似距離観測データＰＲｊｓｒ（ｔ）及び搬送波位相観測データＣＰｊｓｒ（ｔ）を求める機能を備えている。更に各測位衛星６の搭載時計６１及び各モニタ局７の受信機時計７１は正確に時刻が合わせられている。
【００４１】
擬似距離観測データＰＲｉｓｒ（ｔ）の求め方については、背景技術で既に述べたように、サンプリングのタイミングで検出したパルス列のコードに基づいて測位衛星６とモニタ局７間との距離が分る。即ち各測位衛星６の軌道はかなり正確に分っていることから、パルス列の長さが測位衛星６とモニタ局７間との距離（例えば２万ｋｍ程度）に比べて短くても、当該測位信号が何時に送信されたものであるかは把握されているので、パルス列の先頭からの長さを表すコードを検出した時刻と当該測位信号が送信された時刻との差分（時間）と光速とに基づいて測位衛星６とモニタ局７間との距離が計算できる。この距離は、既述のように数３で表されるいわば擬似距離である。これは種々の誤差を含むため、主管制局８に送られてデータ処理され、正確な距離が算出されることになる。
【００４２】
搬送波位相観測データＣＰｉｓｒ（ｔ）は、背景技術で既に述べたよう所定の時刻ｔ０からカウントした周波数ｆｉの搬送波のカウント値を求め、このカウント値から{周波数ｆｉ×（ｔ−ｔ０）}を差し引いた値である。先の数４に示すように搬送波位相観測データＣＰｉｓｒ（ｔ）は、その差分であるｒ（ｔ）／λ−ｒ（ｔ０）／λであり、仮に測位衛星６とモニタ局７との相対位置が変わらなければゼロである。時刻ｔ０は受信機７１が搬送波と同調して検出を開始したタイミングであってもよいが、通常は測位信号は連続して発信されるため、受信機７１側は連続して搬送波を検出しており、従ってある所定時刻をｔ０として、そこから積算した値を用いてもよい。更に搬送波位相観測データＣＰｉｓｒ（ｔ）は受信機７２内の図示されない局部発信器を用いて、前記局部発信器から出力される信号を利用してもよい。この場合、搬送波受信開始（カウント開始）から計測時刻までの間に局部発信器が発信する波数は次のΦ´（ｔ）で表される。
【００４３】
Φ´＝ｆ´（ｔ−ｔ０）＝｛ω´（ｔ−ｔ０）｝／２πｆ´＝ω´／２π但しω´は受信機７２の局部発振器の角周波数、ｆ´は受信機
７２の局部発振器の周波数である。
【００４４】
また、計測時刻ｔにおける受信搬送波と受信機７２の局部発振器が発生した信号との積は次のように表される。
【００４５】
【数１１】
2005062046.tif 000012
受信機７２は下式の右辺の高周波部がフィルタで除去し、ω、ω´、ｔ、ｔ０は既知であることから次の項を記録する。下式の右側の項であるｒ（ｔ）／λ−ｒ（ｔ０）／λが搬送波位相のアンビギュイティとなる。
【００４６】
また、カウント値を読み出すタイミングは、上記のパルス列のコードを検出するタイミングであり、例えばそのタイミングを時刻ｔｋとすると、時刻ｔｋにおける擬似距離観測データＰＲｉｓｒ（ｔ）及び搬送波位相観測データＣＰｉｓｒ（ｔ）が観測されたことになる。
【００４７】
図１中８は主管制局であり、この主管制局８は、観測データ記憶部８１と、合成値演算部８２と、合成値記憶部８３と、状態量推定部８４と、推定状態量記憶部８５と、航法メッセージ作成部８６と、送信機８７と、を備えている。前記観測データ記憶部８１は、各モニタ局７から送られてくる観測データを記憶する記憶部であり。各モニタ局７Ａ〜７Ｃ毎に各測位衛星６Ａ〜６Ｃからの観測データ（３×３＝９組の観測データ）が記憶される。なお、ここでは説明の便宜上測位衛星とモニタ局との数を夫々「３つ」としているが、実際のシステムはこの数に限定されるものではない。前記合成値演算部８２は、各観測データの組毎に、ＣＰｉｓｒ（ｔ）、ＣＰｊｓｒ（ｔ）、ＰＲｉｓｒ（ｔ）、ＰＲｊｓｒ（ｔ）を先の数７に代入して合成値ＩＦＣＰｉｊｒｓｒ（ｔ）及びＩＦＰＡｉｊｒｓｒ（ｔ）を求め、合成値記憶部８３に記憶される。
【００４８】
前記状態量推定部８４は前記合成値に基づいて例えばカルマンフィルタにより状態量を求めるものであり、その入力値と出力値とは図３（ａ）に示す通りである。
【００４９】
ここでいう状態量とは、数７を変形して得られた既述の数８の右辺の項目である。この実施の形態では、当該右辺の項目のうち搬送波位相観測ノイズｅｓｒ（ｔ）及び擬似距離観測ノイズＥｓｒについては平均値をゼロと仮定することによりゼロとして取り扱っている。このためＮｓｒは数８の２段目の式から第２の合成値であるＩＦＰＡｉｊｒｓｒ（ｔ）の値になる。また、光速の対流圏遅延に基づく誤差であるＴｓｒ（ｔ）は地上のモニタ局７に送られる測位信号の入射角に応じて変わってくるので、例えばモニタ局７に対する測位衛星６の相対位置に応じて決まる値として取り扱うことができ、観測の時間帯毎に割り当てられた定数として扱うようにしてもよい。このようなことから、演算サイクルで推定する状態量とは、図３（ｂ）に示すように測位衛星６の位置（例えば、慣性座標系におけるＸ、Ｙ、Ｚ方向の座標）、測位衛星６の速度（Ｘ、Ｙ、Ｚ方向の各速度）、既述のモニタ局ｒの受信機時計誤差＋周波数間バイアス（ｆｉ、ｆｊの測位信号を検出するタイミングの僅かなずれ）に起因する観測データの誤差（Ｄｔｓ（ｔ））、測位衛星ｓの搭載時計誤差＋周波数間バイアス（ｆｉ、ｆｊの測位信号を発信するタイミングの僅かなずれ）に起因する観測データの誤差Ｄｔｓ（ｔ）である。測位衛星６Ａ〜６Ｃの位置については、各モニタ局７Ａ〜７Ｃ及び各測位衛星６Ａ〜６Ｃの組み合わせ毎の観測データが主管制局８に取り込まれているので、各モニタ局７Ａ〜７Ｃ及び各測位衛星６Ａ〜６Ｃ間の距離の推定値が分り、これらの推定値から各測位衛星６Ａ〜６Ｃの位置が推定できる。また、各測位衛星６Ａ〜６Ｃの速度は、モニタ局７における観測のサンプリングの時間間隔と、相前後するサンプリングにおける測位衛星６の各位置とに基づいて推定することができる。
【００５０】
そして、例えば、モニタ局７において時刻ｔｋにて擬似距離及び搬送波位相を観測したとすると、この観測データに基づいて演算された時刻ｔｋにおける合成値ＩＦＣＰｉｊｒｓｒ（ｔ）及びＩＦＰＡｉｊｒｓｒ（ｔ）が状態量推定部８４に入力されるが、これら合成値の分散も入力される。この分散とは受信機等の性能に基づいて決められる合成値の正確性を示す値であり、この例では定数として取り扱っている。
【００５１】
更に前記状態量推定部８４には、当該演算サイクルよりも前の演算サイクル、例えば１つ前の演算サイクルで推定した状態量とこれら状態量の推定値の各々に対する共分散（推定値の正確性を示す値）が入力される。一つ前の演算サイクルで推定した状態量とは、地上のモニタ局７において時刻ｔｋよりも一つ前のサンプリング時刻ｔｋ−１において観測した観測データに基づいて合成値を求め、その合成値と更に前の演算サイクルで求めた状態量推定値等に基づいて推定した状態量、つまり時刻ｔｋ−１における状態量推定値である。
【００５２】
こうして今回（時刻ｔｋ）の合成値及びその分散と、前回（ｔｋ−１）の状態量の推定値及び各推定値の共分散とが状態量推定部８４に入力され、カルマンフィルタによる演算が行われて、今回（時刻ｔｋ）の状態量の推定値と各推定値の共分散とが求められて推定状態量記憶部８５に記憶される。
【００５３】
ここでカルマンフィルタによる演算について簡単に述べておくと、前回（時刻ｔｋ−１）の演算サイクルで推定した測位衛星６の位置と速度とに基づいて今回（時刻ｔｋ）の演算サイクル時における測位衛星６の位置を推定した測位衛星６の位置を推定し、また前回（時刻ｔｋ−１）のＤｔｒ（ｔ）と所定のモデルとに基づいて、今回（時刻ｔｋ）のＤｔｓ（ｔ）を推定する。上述の所定のモデルとは、周波数間バイアスは一定であって定数であることから、定数＋時計誤差のモデルで決まる値になるが、この時計誤差のモデルは時計の種別に応じて時間の関数として表される。そして前回（時刻ｔｋ−１）の状態量推定値に基づいて推定された今回（時刻ｔｋ）の各状態量を数８の上段の式の右辺に代入して、今回（時刻ｔｋ）の第１の合成値ＩＦＣＰｉｊｒｓｒ（ｔ）を求め、得られた第１の合成値と観測データを数７に代入して得た第１の合成値とを比較してその差分を求めると共に、第２の合成値ＩＦＰＡｉｊｒｓｒ（ｔ）についても同様にして差分を求め、これら差分に応じて今回（時刻ｔｋ）の各状態量推定値を修正する。このよにして修正された各状態量が推定状態量記憶部８５に記憶され、次の演算サイクルにて使用されると共に、航法メッサージ作成部８６により航法メッセージの中に取り込まれ、この航法メッセージは送信機８７によって各測位衛星６Ａ〜６Ｃに送信される。
【００５４】
以上の説明から分るように、各モニタ局７Ａ〜７Ｃにおける各測位衛星６Ａ〜６Ｃから観測データに基づいて、そのときの各測位衛星６Ａ〜６Ｃの位置、速度を含む既述の状態が推定できる。
【００５５】
上述の実施の形態によれば、擬似距離及び搬送波位相の観測データを用いて合成値を得るにあたって、擬似距離及び搬送波位相のいずれについても２つの観測データ（周波数ｆｉ、ｆｊに対応する観測データ）を用いているので、ＧＰＳで使用されている二重差の方式（数５）に比べて観測ノイズの分散が４分の１となるので、各状態値を高精度に推定できる。従って測位衛星の軌道を高精度に推定でき、測位衛星を活用したアプリケーション、例えばナビゲーションシステムを信頼性の高いものとすることができる。
【実施例】
【００５６】
日本で計画中の準天頂衛星及びＧＰＳ衛星に本発明の測位衛星の位置検出システムを適用し、準天頂衛星及びＧＰＳ衛星の軌道をシュミレーションにより推定した結果例を図４及び図５に示す。図４の（ａ）、（ｂ）及び（ｃ）は、準天頂衛星の軌道推定において，（ａ）はＸ座標位置の推定誤差、（ｂ）はＹ座標位置の推定誤差、（ｃ）はＺ座標位置の推定誤差を表したものである。ここでいう推定誤差とは、既述のように前回（ｔｋ−１）の状態値に基づいて今回（ｔｋ）推定した状態値と、数８の左辺と右辺との誤差に基づいて状態値を修正して得られた今回（ｔｋ）の状態値との差分である。縦軸はその座標での推定誤差（ｃｍ）であり、横軸は観測時間（ｈ）である。図４の（ａ）、（ｂ）及び（ｃ）において、従来の測位衛星観測システムでは推定誤差が５ｍなのに対して、測定時間を通して推定誤差が０に近い値を示すことから、準天頂衛星の軌道を高精度に推定することができる。
【００５７】
また、図５は、ＧＰＳ衛星の同様の誤差の推定を軌道推定シュミレーションにより求めたものである。図５の（ａ）はＸＹＺ位置推定誤差であり、縦軸は推定誤差（ｃｍ）、横軸は観測時間（ｈ）である。図５の（ｂ）は時計誤差によるバイアス推定誤差であり、縦軸は推定誤差（ｃｍ）、横軸は観測時間（ｈ）である。図５の（ｃ）は地上のモニタ局から見たモニタ局とＧＰＳ衛星との測定角度であり、縦軸は仰角（℃）、横軸は観測時間（ｈ）である。図５の（ａ）及び（ｂ）の各推定誤差（ｃｍ）は図５の（ｃ）に基づいて行っている。よってどの測定角度においても、（ａ）及び（ｂ）は常に０に近い測定誤差であることがわかる。
【図面の簡単な説明】
【００５８】
【図１】本発明に係る人工衛星の位置推定システムの実施の構成を示す構成図である。
【図２】本発明に係る測位衛星から送られてくる測位データのパルス列を示すイメージ図である。
【図３】合成演算部及び状態量推定部の入力、出力の関係と、状態量とを説明する説明図である。
【図４】本発明のシステムを用いて準天頂衛星位置の推定誤差を示す推定図である。
【図５】本発明のシステムを用いてＧＰＳ衛星位置の推定誤差を示す推定図である。
【図６】従来の観測方法において測位衛星とモニタ局とを模式的に示す説明図である。
【図７】従来の観測方法において搬送波位相による距離計測の原理を説明した説明図である。
【図８】従来の観測方法において二重差による観測の様子を示す説明図である。
【符号の説明】
【００５９】
６    測位衛星
６１  搭載時計
６２  受信機
６３  測位信号生成部
６４  送信機
７    モニタ局
７１  受信機時計
７２  受信機
８    主管制局
８１  観測データ記憶部
８２  合成値演算部
８３  合成値記憶部
８４  状態量推定部
８５  推定状態量記憶部
８６  航法メッセージ
８７  送信機
１００  一般ユーザ
１０１  受信機
【図１】
2005062046.tif 000015
【図２】
2005062046.tif 000016
【図３】
2005062046.tif 000017
【図４】
2005062046.tif 000018
【図５】
2005062046.tif 000019
【図６】
2005062046.tif 000020
【図７】
2005062046.tif 000021
【図８】
2005062046.tif 000022
【手続補正書】
【提出日】平成15年9月1日(2003.9.1)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、人工衛星から送信される測位信号の観測データから測位衛星等の人工衛星の位置を推定する人工衛星の位置推定システムに関する。
【背景技術】
【０００２】
米国のＧＰＳ（Ｇｌｏｂａｌ  Ｐｏｓｉｔｉｏｎｉｎｇ  Ｓｙｓｔｅｍ）、ロシアのＧＬＯＮＡＳＳ、ヨーロッパのＧＡＬＩＲＥＯ、更に日本で計画中の準天頂衛星システム等の測位衛星システムを利用して精密な測定位置及び時刻同期を行うには、測位衛星が地球の周りを周回する軌道及び測位衛星に搭載した時計の誤差を高精度に推定することが必要不可欠である。
【０００３】
測位衛星から送信される測位信号を利用して測位衛星の軌道及び測位衛星に搭載した時計誤差を推定する際に必要となるモデルは、測位衛星の運動を表す運動（力学）モデルと測位信号を受信した観測データに含まれる各種の誤差を表現する観測モデルとが必要である。この測位信号は測位衛星と観測点間との距離を測るために使用されるが、測位信号の観測データの中には測位衛星と観測点間との距離データの他に測位衛星に搭載された時計の時間誤差、電離層や対流圏による電波遅延、受信機ノイズなど様々なデータ誤差が含まれる。これらのデータ誤差を補正及び除去することが高精度の測位衛星の軌道及び測位衛星に搭載した時計誤差の推定では必要となる。
【０００４】
ここで図６を用いて測位衛星の軌道を推定する原理について説明する。地球の大気圏外領域を回っている測位衛星１からは測位信号Ｚと呼ばれる一定の信号形式を持つ電波が常時、測位衛星１に搭載された衛星時計１１に同期して送信機１２によって送信される。この測位信号Ｚは、地上のモニタ局２Ａ〜２Ｃに設置された受信機３Ａ〜３Ｃによって受信され、以下の原理で測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を測定する。
【０００５】
先ず、測位衛星１が測位信号Ｚを送信した時刻は、各モニタ局２Ａ〜２Ｃで受信した信号を解読することにより知ることができる（送信時刻（ｔ））。モニタ局２Ａ〜２Ｃで測位信号Ｚを受信した時刻はモニタ局２Ａ〜２Ｃに設置されたモニタ局時計４Ａ〜４Ｃで知ることができる（受信時刻（Ｔ））。このことから測位信号Ｚが測位衛星１からモニタ局２Ａ〜２Ｃまでの空間を伝播した時間（伝播時間（Δｔ））はΔｔ＝Ｔ−ｔとなり、よって測位信号Ｚは光速ｃで空間を伝わることから測位衛星１とモニタ局２Ａ〜２Ｃ間との距離ＤはＤ＝Δｔ×ｃで表すことができる。
【０００６】
このように予め正確な位置が分っている複数のモニタ局２Ａ〜２Ｃからモニタ局２Ａ〜２Ｃと測位衛星１との各距離データＤ１〜Ｄ３を主管制局５に設けられた演算部５１で演算することにより測位衛星１の位置を推定することができる。このシステムにより得られた測位衛星１の位置は、例えば地上の移動体のナビゲーションシステムにも適用される。
【０００７】
上述した測位衛星１の運動を表す運動（力学）モデルは、地球の引力等の物理法則に支配されており、ある時刻での位置、速度が分れば一定時間後の位置、速度を予測することができるので、測位衛星１の軌道推定では、測位信号Ｚを正確に測定することが重要となる。
【０００８】
そこで、測位信号Ｚを利用した測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を測定する方法には、測位信号Ｚの信号形式（一定の長さを持ったパターン）を利用した擬似距離を用いる方法と、測位信号Ｚの搬送波を数える搬送波位相を用いる方法とがある。
【０００９】
この擬似距離を用いる方法では、測位衛星１とモニタ局２Ａ〜２Ｃ間との距離を直接測ることができるが、計測精度が信号パターンを構成するデータの単位長さに制限されるため距離の誤差が大きくなる。一方前記搬送位相を用いる方法では、計測精度は測位信号Ｚの波長以下となり測定精度を上げることができるが、受信電波の波を数えるため波長の整数倍の不確定性（アンビギュイティ）が残ってしまう。
【００１０】
また、前記測位衛星１からは、詳細には、異なった周波数で複数の測位信号Ｚが同時に送信されており、これら複数の測位信号Ｚにおいて、各々の測位信号Ｚ毎に擬似距離及び搬送波位相が求められる。しかしながら、測位衛星１から送信された測位信号Ｚを地上のモニタ局２Ａ〜２Ｃに設置された受信機３Ａ〜３Ｃで受信すると測位衛星１と地上のモニタ局２Ａ〜２Ｃ間との距離データＤ１〜Ｄ３の中には（表１）に示す各種誤差が混入する。
【００１１】
【表１】
2005062046.tif 000023
これらのモニタ局２Ａ〜２Ｃで観測された擬似距離観測データと、搬送波位相観測データと、に含まれる観測誤差を表した観測モデルは、一般的に数３のように表される。
【００１２】
【数３】
2005062046.tif 000024
ここで、上述の数３を用いて測位衛星１の軌道決定のための推測処理をどのようにして行うのかを説明する。数３の式の左辺は観測値であり、例えば数３中のＰＲｉｓｒ（ｔ）は擬似距離観測値である。この擬似距離観測値は次のようにして求められる。測位衛星１は、例えば地上から２万ｋｍの高度の所を移動しており、この測位衛星１からパルスパターンを含む例えば長さ３００〜５００ｋｍ程度のパルス列を地上に送る。このパルスパターンは、パルス列の先頭からの距離に対応したパターンとして配列され、例えばパルス列の先頭からＡｃｍの位置にあればＡｃｍを表すパルスパターンが形成される。つまり、このパルス列は測位衛星１からいわば長さ３００〜５００ｋｍの物差しが飛来してくるのと同じである。
【００１３】
一方測位衛星１と地上のモニタ局２Ａ〜２Ｃとにおいては、正確に時計合わせされたモニタ局時計４Ａ〜４Ｃが夫々設けられ、測位衛星１からパルス列を発信するタイミング及びモニタ局２Ａ〜２Ｃでこのパルス列を検出するタイミングが予め決められており、一定周期で発信（検出）される。従ってモニタ局２Ａ〜２Ｃの検出のタイミングでパルス列のあるパルスパターンを把握すれば、その物差しのいくつの値をどのタイミングで検出したかが分る。測位衛星１とモニタ局２Ａ〜２Ｃとの距離は運動（力学）モデルにより概ね分っているので（例えば数ｍの誤差で分っているので）、この測位衛星１とモニタ局２Ａ〜２Ｃとの距離に比べてパルス列の長さが短くても、当該パルス列の中のパルスパターンをある時刻で検出すれば、そのパルスパターンはパルス列を物差しと見立てたときのスケール値に相当するので、当該パルス列の先頭からどれだけの時間だけ遅れて送信されているかが分り、結局パルス列の先頭をモニタ局２Ａ〜２Ｃが受信した時刻が分る。この結果、パルス列の先頭が飛来してくる時間（（モニタ局２Ａ〜２Ｃで受信した時刻）−（測位衛星１から送信された時刻））が検出され、この時間に光速ｃを掛けることで測位衛星１とモニタ局２Ａ〜２Ｃとの距離を算出することができる。しかしながら、実際には数３から分るようにこの距離は、種々の誤差を含んでいる。
【００１４】
また上記の数３中のＣＰｉｓｒ（ｔ）を用いても同様にして測位衛星１とモニタ局２Ａ〜２Ｃとの距離ρｓｒ（ｔ）を用いることができる。この場合には、擬似距離の観測データを用いる代わりに搬送波の位相を観測した観測データを用いている。この搬送波とは、測位衛星１から送られてくる既述のパルス列の搬送波であり、搬送波位相観測値ＣＰｉｓｒ（ｔ）とは、ある時刻ｔ０からある時刻ｔまでに地上のモニタ局２Ａ〜２Ｃで搬送波をカウントし、搬送波をカウントした値から時刻ｔ０における測位衛星１の位置と時刻ｔにおける測位衛星１の位置とが変わることに基づいて搬送波の波数値を求めたものである。この搬送波位相による距離計測の原理について図７と数４を参照しながら詳しく説明する。
【００１５】
先ず図７で、モニタ局２に設けられた受信機３で搬送波を受信開始したモニタ局時計４の時刻をｔ０とする。この時点で受信した搬送波が測位衛星１に設けられた送信機１２から送信されたときの衛星時計１の時刻は、測位衛星１とモニタ局２間との距離をｒ（ｔ０）とするとｔ０−｛ｒ（ｔ０）／ｃ｝となり、この時点では受信機３で受信した受信波数は０個である。この時点から同時に受信機３内部の図示されない局部発振機の波数積算も開始される。ここで搬送波の計測時刻をｔとする。この時点で受信機３で受信した搬送波が測位衛星１に設けられた送信機１２から送信されたときの衛星時計１１の時刻は、測位衛星１とモニタ局２間との距離をｒ（ｔ）とするとｔ−｛ｒ（ｔ）／ｃ｝となり、受信開始時刻から計測時刻までに受信した波数Φ（ｔ）は数４で表される。
【００１６】
【数４】
2005062046.tif 000025
一般に静止系であれば、所定時間に受信する電波の数は一定であり、数４の右辺の前段の項である１／２πω・（ｔ−ｔ０）だけで決定されるが、測位衛星１は移動しているため、後段の項である{ｒ（ｔ）／λ}−{ｒ（ｔ０）／λ}が含まれてくる。この値が搬送波の波数値であり、ＣＰｉｓｒ（ｔ）に相当するものである。しかしながら、上述したようにこの搬送波の位相を利用する方法では、計測精度は測位信号Ｚの波長以下となり高精度であるが、受信電波の波数を数えるため波長の整数倍の不確定性（アンビギュイティ）が残ってしまう。
【００１７】
測位衛星１の軌道推定は、測位衛星１と地上のモニタ局２Ａ〜２Ｃ間との距離データＤ１〜Ｄ３を利用して測位衛星１の位置、速度を求めることであるが、測位信号Ｚの観測データ中には測位衛星１とモニタ局２Ａ〜２Ｃ間との距離以外の誤差が含まれているため、これらの誤差の値も同時に推定する必要が生じてくる。しかし、前記数３で表した観測モデルのままでは未知数が多く夫々の観測誤差を分離して推定することができない。そこで擬似距離観測データと搬送波位相観測データとを次のように合成した擬似的な観測データを利用して未知数を減らす方法が一般的な手法である（例えば、非特許文献１参照。）。この一般的な手法を図８と数５を参照しながら詳しく説明する。
【００１８】
【数５】
2005062046.tif 000026
この方法は、図８に示すように２つのモニタ局ｒ（ｒ´）と二つの測位衛星ｓ（ｓ´）とを用い、更に各測位衛星ｓ（ｓ´）のいずれからも送信される周波数ｆｉの搬送波と周波数ｆｊの搬送波とを利用する方法（二重差観測データ）であり、先の数３で表される擬似距離観測値と搬送波位相の波数値とを数５のように定義する。これらの式においてＣＰｉｓｒ（ｔ）は、モニタ局ｒにおいて、測位衛星ｓから送信される周波数ｆｉの搬送波について数４のように時刻ｔ〜ｔ０間で搬送波をカウントし、数４中のωは既値であることから、後段の項を求めた値である。また、ＣＰｊｓｒ´（ｔ）は、モニタ局ｒ´において測位衛星ｓから送信される周波数ｆｊの搬送波について求めた値である。更にＰＲｉｓｒ（ｔ）は、測位衛星ｓから周波数ｆｉの搬送波に乗せられて送られる既述のパルス列中のパルスパターンをモニタ局ｒについて検出し、その検出結果に基づいて得られた擬似距離観測値である。
【００１９】
そして、数５の右辺の各項に、数３で定義される値（右辺）を代入することにより測位衛星ｓ（ｓ´）と地上のモニタ局ｒ（ｒ´）間との距離以外の誤差要因（観測誤差）は時計誤差、電離層遅延及び周波数間バイアスが合成により打消されて、対流圏遅延、搬送波位相アンビギュイティ及び観測ノイズのみが残ることになる（数６）。
【００２０】
【数６】
2005062046.tif 000027
もう少し詳しく説明すると、数５の右辺は、各モニタ局ｒ（ｒ´）における観測値を演算した値であるから、各モニタ局ｒ（ｒ´）における観測値{ＰＲｉｓｒ（ｔ）等}を地上回線で主管制局のコンピュータに取り込むことによって、数５の左辺であるＩＣＰｉｊｓｓ´ｒｒ´（ｔ）及びＩＡＭｉｊｓｓ´ｒｒ´（ｔ）が求まる。
【００２１】
今、時刻ｔｋにおける擬似距離観測値と搬送波の波数値との各観測値を用いて数５を演算した結果を夫々ＩＣＰｉｊｓｓ´ｒｒ´（ｔｋ）及びＩＡＭｉｊｓｓ´ｒｒ´（ｔｋ）とする。なお、時刻ｔｋにおける搬送波の波数値は、時刻ｔ０からの搬送波の積算値に基づいて求められるが、ω（ｔｋ−ｔ０）の値は既にわかっているので、この値をキャンセルすることにより、結局、数４の{ｒ（ｔｋ）／λ}−{ｒ（ｔ０）／λ}のみが波数値として求まることになる。
【００２２】
一方コンピュータは、時刻ｔｋよりも前の時刻ｔｋ−１において推定した数６の右辺の測位衛星ｓ（ｓ´）とモニタ局ｒ（ｒ´）間との距離及び推定状態量に基づいて時刻ｔｋの右辺の値を求める。ここで測位衛星ｓ（ｓ´）とモニタ局ｒ（ｒ´）間との距離であるρｓｒ（ｔ）については初期値としては適当な値を用いる。但し、測位衛星ｓ（ｓ´）の軌道は予めかなりの精度で把握されているので、適当な値とはいっても高い精度を持つ値である。また、他の推定状態量のうちＴｓｒ（ｔ）は信頼できる理論的モデル（時間の関数）を用い、時刻ｔｋ−１のときの対流圏遅延を用いて、時刻ｔｋのときの対流圏遅延（Ｔｓｒ（ｔ））を求めることができる。即ち、対流圏遅延は、モニタ局から見た測位衛星の仰角に応じて決まるため、モニタ局と測位衛星との相対位置に応じて決まる値として取り扱うことができる。またＮｉｓｒは連続して電波受信中は定数として考えることができ、更にＥｉｓｒ（ｔ）、ｅｉｓｒ（ｔ）はいずれも平均値は０と仮定して取り扱うことができる。こうして数６を用いて、数５の右辺が求まる。先に求めたＩＣＰｉｊｓｓ´ｒｒ´（ｔ）及びＩＡＭｉｊｓｓ´ｒｒ´（ｔ）は数５の左辺であるため、本来は左辺の観測値と右辺の演算値とは一致するはずであるが、実際にはモデルの誤差分や種々の仮定が含まれているため、また時刻ｔｋの値があくまでも推定値であるため、観測値と演算値とは異なってくる。
【００２３】
そこで、コンピュータは数５の夫々について左辺と右辺との差（偏差）を求める。そして、これら偏差に応じて時刻ｔにおけるρｓｒ（ｔ）、Ｔｓｒ（ｔ）及びＮｉｓｒ（ｔ）の値を修正する。実際には、主管制局では各モニタ局毎でこれらの推定値が求まり、時刻ｔｋにおける推定値としてメモリ内に格納され、次の時刻ｔｋ＋１における推定に用いられる。また、推定初期値においては、誤差が大きい場合でも時間が経過するにつれて精度が高くなり前記誤差は小さくなっていく。時刻ｔｋにおける測位衛星の位置、速度は３ヶ所以上のモニタ局と測位衛星との距離から求まることになる。
【００２４】
【非特許文献１】Global Positioning System:Theory and Applications Volume1 Part2(American Institute of Aernautics and Astronautics 1996)
【発明の開示】
【発明が解決しようとする課題】
【００２５】
しかしながら、擬似距離観測データと搬送波位相観測データとを上述のように合成して未知数を減らす一般的な手法（二重差観測データ）では、二重差の使用により測位衛星の軌道は決定可能となる８つの観測データを合成することにより観測ノイズの分散は１つの観測データの場合に比べて８倍となり精密な軌道決定には不利である。この二重差の不利を補うためにＧＰＳの軌道決定において、一般的に全世界の観測網を利用して長時間連続的に測位衛星の観測及び軌道推定を行うことで除々に高精度な軌道に収束させる方法が用いられているが、地上に多数のモニタ局が必要となるため、多大なコストが掛かるという問題がある。
【００２６】
本発明は係る事情に鑑みてなされたものであって、本発明の目的は、測位衛星等の人工衛星の位置を高精度に推定できる人工衛星の位置推定システムを提供することにある。
【課題を解決するための手段】
【００２７】
本発明は、人工衛星から送信され、その先頭からの位置に対応する位置コードが配列された所定の長さのパルス列を搬送波に載せた測位信号を地上のモニタ局にて受信し、モニタ局にて検出した位置コードと光速とに基づいて人工衛星とモニタ局との擬似距離観測データを求めると共に、
モニタ局にて検出した前記搬送波の波数を所定の時刻ｔ0から前記位置コードをサンプリングした時刻ｔkまでの時間をカウントし、そのカウント値と前記搬送波の周波数と前記カウントした時間とに基づいて、時刻ｔ0におけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値と時刻ｔkにおけるモニタ局及び人工衛星間の距離を搬送波の波長で割り算した値との差分である搬送波位相観測データを求め、
前記擬似距離観測データと前記搬送波位相観測データとに基づいて、人工衛星とモニタ局との距離を求め、当該人工衛星と複数のモニタ局の各々との距離に基づいて当該人工衛星の位置を求めるシステムにおいて、
イ．  人工衛星から送信される周波数ｆiの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｉs r(t)及び搬送波位相観測データＣＰｉs r(t)と、
前記測位信号と同時に送られる、周波数ｆjの搬送波に前記パルス列を載せた測位信号を観測して得た擬似距離観測データＰＲｊs r(t)及び搬送波位相観測データＣＰｊs r(t)と、
を下記数７に代入して第１の合成値及び第２の合成値を算出する合成値演算部と、
ロ．  時刻ｔkに観測した観測データに基づいて計算された第１の合成値及び第２の合成値と、
数７を変形して得られた下記数８と、
時刻ｔkのサンプリングよりも一つ前のサンプリングのタイミングｔk-1にて観測した観測データを用いて計算された第１の合成値及び第２の合成値に基づいて推定した人工衛星の位置と、
に基づいて、時刻ｔkにおける人工衛星の位置を推定する推定部と、を備えたことを特徴とする人工衛星の位置推定システムである。
【００２８】
記
【００２９】
【数７】
2005062046.tif 000028
ただし、λi及びλjは夫々周波数ｆiの搬送波の波長及び周波数ｆjの搬送波の波長である。
【００３０】
【数８】
2005062046.tif 000029
ただし、ρ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における人工衛星の位置を推定することが好ましい。
【００３２】
ここで数７と数８との関係について説明しておく。数７の右辺の各項に数３の右辺を代入することにより数９が得られる。
【００３３】
【数９】
2005062046.tif 000030
数９の｛  ｝内を数１０のように合成された１つの未知数として取り扱うと数８が得られる。
【００３４】
【数１０】
2005062046.tif 000031
数８において擬似距離観測ノイズＥ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に観測した観測データに基づいて計算された第１の合成値及び第２の合成値とこれら合成値の分散を用いて、今回の演算サイクルにおける状態量及び分散を推定する。なお第１の合成値及び第２の合成値の分散は定数として扱うことができる。
【００３６】
また前記推定部は、人工衛星とモニタ局との組毎に各状態量を推定するので、推定演算の過程では例えば人工衛星とモニタ局との距離を推定し、各人工衛星の位置についてはその人工衛星と３つ以上のモニタ局との距離から計算することで推定している。こうして各人工衛星の位置を推定することができ、例えば推定された位置と一つ前の演算サイクルで求めた位置とから人工衛星の速度も推定できる。
【発明の効果】
【００３７】
ＧＰＳで行われている二重差による手法（数５の手法）では、擬似距離及び搬送波位相の各々について８つの観測データを合成していたが、本発明では２つの観測データを合成した合成値に基づいて測位衛星等の人工衛星の位置が推定でき、従って観測ノイズの分散は二重差の手法に比べて４分の１となり、高精度な位置推定を行うことができる。更に測位衛星等の人工衛星の時計誤差及びモニタ局の時計誤差の推定も高精度に行うことができる。
【発明を実施するための最良の形態】
【００３８】
本発明の人工衛星の位置推定システムについて、図１〜図３に基づいて説明を行う。図１は本発明の実施の形態に係る人工衛星である測位衛星の位置測定システムを含む測位衛星の軌道推定システムの説明図である。図１中６は測位衛星であり、６Ａ、６Ｂ、６Ｃは３機の測位衛星の夫々を指している。前記測位衛星６は搭載時計６１、受信機６２、測位信号生成部６３及び送信機６４を備えている。前記測位信号生成部６３は、搭載時計６１からクロック信号に基づいて、周波数ｆｉの搬送波で送られる測位信号であるパルス列と周波数ｆｊの搬送波で送られる測位信号であるパルス列とを作成して前記送信機６４に受け渡す機能を備えている。これらパルス列は例えば３００ｋｍ〜５００ｋｍ程度の長さに作られており、パルス列からの先頭からの位置に対応するコードが配列されている。図２はパルス列を示すイメージ図であり、当該パルス列を発信する測位衛星６の識別コード（ＩＤ）が配置されると共に、その後には各位置に対応したコード（Ｐ１、Ｐ２、…Ｐｎ）が配置されている。各コードの長さは例えば２０ｃｍ程度とされる。
【００３９】
なお、測位衛星６は後述の主管制局８から、Ｘ、Ｙ、Ｚ座標における自己の正確な位置（Ｘ、Ｙ、Ｚ）と速度（Ｖｘ、Ｖｙ、Ｖｚ）とを特定可能なデータを含む航法メッセージを受け取り、その航法メッセージを測位信号と共に一般ユーザ１００の受信機１０１に配信する機能を備えている。前記一般ユーザ１００は、例えば車両等の移動体であり、３機の測位衛星６Ａ〜６Ｃから送られてくる航法メッセージと測位信号とに基づいて、各測位衛星６Ａ〜６Ｃまでの各距離が分るので、当該一般ユーザ１００（移動体）の位置を知ることができる。
【００４０】
また、図１中７は地上に設けられたモニタ局であり、７Ａ〜７Ｃは３個のモニタ局を夫々示している。モニタ局７は受信機時計７１及び受信機７２を備えており、受信機時計７１から出力されるクロック信号に基づいて測位衛星６から送られてくる測位信号（パルス列及び搬送波）を検出し、周波数ｆｉの測位信号に基づいて擬似距離観測データＰＲｉｓｒ（ｔ）及び搬送波位相観測データＣＰｉｓｒ（ｔ）を求め、また、周波数ｆｊの測位信号に基づいて擬似距離観測データＰＲｊｓｒ（ｔ）及び搬送波位相観測データＣＰｊｓｒ（ｔ）を求める機能を備えている。更に各測位衛星６の搭載時計６１及び各モニタ局７の受信機時計７１は正確に時刻が合わせられている。
【００４１】
擬似距離観測データＰＲｉｓｒ（ｔ）の求め方については、背景技術で既に述べたように、サンプリングのタイミングで検出したパルス列のコードに基づいて測位衛星６とモニタ局７間との距離が分る。即ち各測位衛星６の軌道はかなり正確に分っていることから、パルス列の長さが測位衛星６とモニタ局７間との距離（例えば２万ｋｍ程度）に比べて短くても、当該測位信号が何時に送信されたものであるかは把握されているので、パルス列の先頭からの長さを表すコードを検出した時刻と当該測位信号が送信された時刻との差分（時間）と光速とに基づいて測位衛星６とモニタ局７間との距離が計算できる。この距離は、既述のように数３で表されるいわば擬似距離である。これは種々の誤差を含むため、主管制局８に送られてデータ処理され、正確な距離が算出されることになる。
【００４２】
搬送波位相観測データＣＰｉｓｒ（ｔ）は、背景技術で既に述べたよう所定の時刻ｔ０からカウントした周波数ｆｉの搬送波のカウント値を求め、このカウント値から{周波数ｆｉ×（ｔ−ｔ０）}を差し引いた値である。先の数４に示すように搬送波位相観測データＣＰｉｓｒ（ｔ）は、その差分であるｒ（ｔ）／λ−ｒ（ｔ０）／λであり、仮に測位衛星６とモニタ局７との相対位置が変わらなければゼロである。時刻ｔ０は受信機７１が搬送波と同調して検出を開始したタイミングであってもよいが、通常は測位信号は連続して発信されるため、受信機７１側は連続して搬送波を検出しており、従ってある所定時刻をｔ０として、そこから積算した値を用いてもよい。更に搬送波位相観測データＣＰｉｓｒ（ｔ）は受信機７２内の図示されない局部発信器を用いて、前記局部発信器から出力される信号を利用してもよい。この場合、搬送波受信開始（カウント開始）から計測時刻までの間に局部発信器が発信する波数は次のΦ´（ｔ）で表される。
【００４３】
Φ´＝ｆ´（ｔ−ｔ０）＝｛ω´（ｔ−ｔ０）｝／２πｆ´＝ω´／２π但しω´は受信機７２の局部発振器の角周波数、ｆ´は受信機７２の局部発振器の周波数である。
【００４４】
また、計測時刻ｔにおける受信搬送波と受信機７２の局部発振器が発生した信号との積は次のように表される。
【００４５】
【数１１】
2005062046.tif 000032
受信機７２は下式の右辺の高周波部がフィルタで除去し、ω、ω´、ｔ、ｔ０は既知であることから、｛ｒ（ｔ）／λ｝−｛ｒ（ｔ０）／λ｝を記録する。この式の右側の項である｛ｒ（ｔ０）／λ｝が搬送波位相のアンビギュイティとなる。
【００４６】
また、カウント値を読み出すタイミングは、上記のパルス列のコードを検出するタイミングであり、例えばそのタイミングを時刻ｔｋとすると、時刻ｔｋにおける擬似距離観測データＰＲｉｓｒ（ｔ）及び搬送波位相観測データＣＰｉｓｒ（ｔ）が観測されたことになる。
【００４７】
図１中８は主管制局であり、この主管制局８は、観測データ記憶部８１と、合成値演算部８２と、合成値記憶部８３と、状態量推定部８４と、推定状態量記憶部８５と、航法メッセージ作成部８６と、送信機８７と、を備えている。前記観測データ記憶部８１は、各モニタ局７から送られてくる観測データを記憶する記憶部であり。各モニタ局７Ａ〜７Ｃ毎に各測位衛星６Ａ〜６Ｃからの観測データ（３×３＝９組の観測データ）が記憶される。なお、ここでは説明の便宜上測位衛星とモニタ局との数を夫々「３つ」としているが、実際のシステムはこの数に限定されるものではない。前記合成値演算部８２は、各観測データの組毎に、ＣＰｉｓｒ（ｔ）、ＣＰｊｓｒ（ｔ）、ＰＲｉｓｒ（ｔ）、ＰＲｊｓｒ（ｔ）を先の数７に代入して合成値ＩＦＣＰｉｊｓｒ（ｔ）及びＩＦＰＡｉｊｓｒ（ｔ）を求め、合成値記憶部８３に記憶される。
【００４８】
前記状態量推定部８４は前記合成値に基づいて例えばカルマンフィルタにより状態量を求めるものであり、その入力値と出力値とは図３（ａ）に示す通りである。
【００４９】
ここでいう状態量とは、数７を変形して得られた既述の数８の右辺の項目である。この実施の形態では、当該右辺の項目のうち搬送波位相観測ノイズｅｓｒ（ｔ）及び擬似距離観測ノイズＥｓｒについては平均値をゼロと仮定することによりゼロとして取り扱っている。このためＮｓｒは数８の２段目の式から第２の合成値であるＩＦＰＡｉｊｓｒ（ｔ）の値になる。また、光速の対流圏遅延に基づく誤差であるＴｓｒ（ｔ）は地上のモニタ局７に送られる測位信号の入射角に応じて変わってくるので、例えばモニタ局７に対する測位衛星６の相対位置に応じて決まる値として取り扱うことができ、観測の時間帯毎に割り当てられた定数として扱うようにしてもよい。このようなことから、演算サイクルで推定する状態量とは、図３（ｂ）に示すように測位衛星６の位置（例えば、慣性座標系におけるＸ、Ｙ、Ｚ方向の座標）、測位衛星６の速度（Ｘ、Ｙ、Ｚ方向の各速度）、既述のモニタ局ｒの受信機時計誤差＋周波数間バイアス（ｆｉ、ｆｊの測位信号を検出するタイミングの僅かなずれ）に起因する観測データの誤差（Ｄｔｒ（ｔ））、測位衛星ｓの搭載時計誤差＋周波数間バイアス（ｆｉ、ｆｊの測位信号を発信するタイミングの僅かなずれ）に起因する観測データの誤差Ｄｔｓ（ｔ）、搬送波位相アンビギュイティに起因する誤差である。測位衛星６Ａ〜６Ｃの位置については、各モニタ局７Ａ〜７Ｃ及び各測位衛星６Ａ〜６Ｃの組み合わせ毎の観測データが主管制局８に取り込まれているので、各モニタ局７Ａ〜７Ｃ及び各測位衛星６Ａ〜６Ｃ間の距離の推定値が分り、これらの推定値から各測位衛星６Ａ〜６Ｃの位置が推定できる。
【００５０】
そして、例えば、モニタ局７において時刻ｔｋにて擬似距離及び搬送波位相を観測したとすると、この観測データに基づいて演算された時刻ｔｋにおける合成値ＩＦＣＰｉｊｓｒ（ｔ）及びＩＦＰＡｉｊｒｓｒ（ｔ）が状態量推定部８４に入力されるが、これら合成値の分散も入力される。この分散とは受信機等の性能に基づいて決められる合成値の正確性を示す値であり、この例では定数として取り扱っている。
【００５１】
更に前記状態量推定部８４には、当該演算サイクルよりも前の演算サイクル、例えば１つ前の演算サイクルで推定した状態量とこれら状態量の推定値の各々に対する共分散（推定値の正確性を示す値）が入力される。一つ前の演算サイクルで推定した状態量とは、地上のモニタ局７において時刻ｔｋよりも一つ前のサンプリング時刻ｔｋ−１において観測した観測データに基づいて合成値を求め、その合成値と更に前の演算サイクルで求めた状態量推定値等に基づいて推定した状態量、つまり時刻ｔｋ−１における状態量推定値である。
【００５２】
こうして今回（時刻ｔｋ）の合成値及びその分散と、前回（ｔｋ−１）の状態量の推定値及び各推定値の共分散とが状態量推定部８４に入力され、カルマンフィルタによる演算が行われて、今回（時刻ｔｋ）の状態量の推定値と各推定値の共分散とが求められて推定状態量記憶部８５に記憶される。
【００５３】
ここでカルマンフィルタによる演算について簡単に述べておくと、前回（時刻ｔｋ−１）の演算サイクルで推定した測位衛星６の位置と速度とに基づいて今回（時刻ｔｋ）の演算サイクル時における測位衛星６の位置と速度とを推定し、また前回（時刻ｔｋ−１）のＤｔｒ（ｔ）と所定のモデルとに基づいて、今回（時刻ｔｋ）のＤｔｒ（ｔ）を推定し、前回（時刻ｔｋ−１）のＤｔｓ（ｔ）と所定のモデルとに基づいて今回（時刻ｔｋ）のＤｔｓ（ｔ）を推定する。上述の所定のモデルとは、周波数間バイアスは一定であって定数であることから、定数＋時計誤差のモデルで決まる値になるが、この時計誤差のモデルは時計の種別に応じて時間の関数として表される。そして前回（時刻ｔｋ−１）の状態量推定値に基づいて推定された今回（時刻ｔｋ）の各状態量を数８の上段の式の右辺に代入して、今回（時刻ｔｋ）の第１の合成値ＩＦＣＰｉｊｓｒ（ｔ）を求め、得られた第１の合成値と観測データを数７に代入して得た第１の合成値とを比較してその差分を求めると共に、第２の合成値ＩＦＰＡｉｊｓｒ（ｔ）についても同様にして差分を求め、これら差分に応じて今回（時刻ｔｋ）の各状態量推定値を修正する。このよにして修正された各状態量が推定状態量記憶部８５に記憶され、次の演算サイクルにて使用されると共に、航法メッサージ作成部８６により航法メッセージの中に取り込まれ、この航法メッセージは送信機８７によって各測位衛星６Ａ〜６Ｃに送信される。
【００５４】
以上の説明から分るように、各モニタ局７Ａ〜７Ｃにおける各測位衛星６Ａ〜６Ｃから観測データに基づいて、そのときの各測位衛星６Ａ〜６Ｃの位置、速度を含む既述の状態が推定できる。
【００５５】
上述の実施の形態によれば、擬似距離及び搬送波位相の観測データを用いて合成値を得るにあたって、擬似距離及び搬送波位相のいずれについても２つの観測データ（周波数ｆｉ、ｆｊに対応する観測データ）を用いているので、ＧＰＳで使用されている二重差の方式（数５）に比べて観測ノイズの分散が４分の１となるので、各状態値を高精度に推定できる。従って測位衛星の軌道を高精度に推定でき、測位衛星を活用したアプリケーション、例えばナビゲーションシステムを信頼性の高いものとすることができる。
【実施例】
【００５６】
日本で計画中の準天頂衛星及びＧＰＳ衛星に本発明の測位衛星の位置検出システムを適用し、準天頂衛星及びＧＰＳ衛星の軌道をシミュレーションにより推定した結果例を図４及び図５に示す。図４の（ａ）、（ｂ）及び（ｃ）は、準天頂衛星の軌道推定において，（ａ）はＸ座標位置の推定誤差、（ｂ）はＹ座標位置の推定誤差、（ｃ）はＺ座標位置の推定誤差を表したものである。この例ではシミュレーションにより測位衛星の真の軌道は把握されているので、ここでいう推定誤差とは、指定した位置と真の位置との誤差である。縦軸はその座標での推定誤差（ｃｍ）であり、横軸は観測時間（ｈ）である。図４の（ａ）、（ｂ）及び（ｃ）において、従来の測位衛星観測システムでは推定誤差が５ｍなのに対して、測定時間を通して推定誤差が０に近い値を示すことから、準天頂衛星の軌道を高精度に推定することができる。
【００５７】
また、図５は、ＧＰＳ衛星の同様の誤差の推定を軌道推定シミュレーションにより求めたものである。図５の（ａ）はＸＹＺ位置推定誤差であり、縦軸は推定誤差（ｃｍ）、横軸は観測時間（ｈ）である。図５の（ｂ）は時計誤差によるバイアス推定誤差であり、縦軸は推定誤差（ｃｍ）、横軸は観測時間（ｈ）である。図５の（ｃ）は地上のモニタ局から見たモニタ局とＧＰＳ衛星との測定角度であり、縦軸は仰角（℃）、横軸は観測時間（ｈ）である。図５の（ａ）及び（ｂ）の各推定誤差（ｃｍ）は図５の（ｃ）に基づいて行っている。よってどの測定角度においても、（ａ）及び（ｂ）は常に０に近い測定誤差であることがわかる。
【図面の簡単な説明】
【００５８】
【図１】本発明に係る人工衛星の位置推定システムの実施の構成を示す構成図である。
【図２】本発明に係る測位衛星から送られてくる測位データのパルス列を示すイメージ図である。
【図３】合成演算部及び状態量推定部の入力、出力の関係と、状態量とを説明する説明図である。
【図４】本発明のシステムを用いて準天頂衛星位置の推定誤差を示す推定図である。
【図５】本発明のシステムを用いてＧＰＳ衛星位置の推定誤差を示す推定図である。
【図６】従来の観測方法において測位衛星とモニタ局とを模式的に示す説明図である。
【図７】従来の観測方法において搬送波位相による距離計測の原理を説明した説明図である。
【図８】従来の観測方法において二重差による観測の様子を示す説明図である。
【符号の説明】
【００５９】
６    測位衛星
６１  搭載時計
６２  受信機
６３  測位信号生成部
６４  送信機
７    モニタ局
７１  受信機時計
７２  受信機
８    主管制局
８１  観測データ記憶部
８２  合成値演算部
８３  合成値記憶部
８４  状態量推定部
８５  推定状態量記憶部
８６  航法メッセージ
８７  送信機
１００  一般ユーザ
１０１  受信機
【手続補正２】
【補正対象書類名】図面
【補正対象項目名】全図
【補正方法】変更
【補正の内容】
【図１】
2005062046.tif 000033
【図２】
2005062046.tif 000034
【図３】
2005062046.tif 000035
【図４】
2005062046.tif 000036
【図５】
2005062046.tif 000037
【図６】
2005062046.tif 000038
【図７】
2005062046.tif 000039
【図８】
2005062046.tif 000040