搬送波位相式測位装置
【課題】サイクルスリップが発生した場合であっても、整数値バイアスの確定処理を一からやり直す必要の無い搬送波位相式測位装置の提供。
【解決手段】移動局及び既知点で測定される衛星信号の搬送波位相の積算値又はその位相差を観測量とし、移動局の位置と搬送波位相の積算値に含まれる整数値バイアス又はその位相差とを状態変数とし、複数の観測周期での衛星データに基づいて前記状態変数を推定して移動局の測位を行う搬送波位相式測位装置であって、位相積算値のサイクルスリップの有無を検出するサイクルスリップ検出手段を備え、サイクルスリップが検出された場合は、サイクルスリップの検出された衛星に係る状態変数だけを初期化して測位を継続する。
【解決手段】移動局及び既知点で測定される衛星信号の搬送波位相の積算値又はその位相差を観測量とし、移動局の位置と搬送波位相の積算値に含まれる整数値バイアス又はその位相差とを状態変数とし、複数の観測周期での衛星データに基づいて前記状態変数を推定して移動局の測位を行う搬送波位相式測位装置であって、位相積算値のサイクルスリップの有無を検出するサイクルスリップ検出手段を備え、サイクルスリップが検出された場合は、サイクルスリップの検出された衛星に係る状態変数だけを初期化して測位を継続する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、衛星信号の搬送波位相の積算値を測定して移動局の測位を行う搬送波位相式測位装置に関する。
【背景技術】
【０００２】
近年、測量の分野では、搬送波位相によるＧＰＳ測量が広く利用されている。この搬送波位相によるＧＰＳ測量では、基準側の受信機と測位側の受信機とが、複数の衛星から送られる衛星信号を同時に受信し、基準側と測位側とで各衛星信号の搬送波位相の積算値をそれぞれ独立に算出する。この搬送波位相の積算値（以下、単に「位相積算値」という）には、搬送波の波長の整数倍に相当する不確定な要素（以下、「整数値バイアス」という）が含まれているが、この整数値バイアスは、時計誤差等とは異なり、位相積算値の一重位相差や二重位相差を取ることによっても消去することができない。このため、ＧＰＳ測量の分野において、位相積算値の三重位相差を取ることで不確定要素である整数値バイアスを消去することや、整数値バイアスそのものを求める技術が提案されている。
【０００３】
ここで、整数値バイアスを確定する技術として、カルマンフィルタを用いる技術が知られている（例えば、特許文献１参照）。この技術では、測位側の位置と整数値バイアスを状態変数とし、基準側に対する測位側の位相積算値の一重位相差を観測量として、観測を重ねる毎に前記状態変数を更新する追尾フィルタが構成される。また、整数値バイアスを確定するその他の技術として、整数値バイアスを含んだ搬送波の二重位相差を用いて、最小二乗法により所定の条件で二重位相差の整数値バイアスを求める技術が知られている（例えば、特許文献２参照）。
【特許文献１】特開２００４−７７２２８号公報
【特許文献２】特開２００３−９８２４５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
ところで、上述のような従来技術は、少なくとも５つの衛星からの衛星信号が、整数値バイアスの確定処理中に継続して（サイクルスリップ無しに）受信されることを前提としている。しかしながら、実際には、整数値バイアス確定後若しくは確定中にある衛星からの衛星信号（電波）の途切れが生じる場合（サイクルスリップ）がある。このサイクルスリップが起こった場合、測定した位相積算値に整数波数分のジャンプが生じうるので、整数値バイアスを確定するための処理を一からやり直さなければならないという不都合がある。
【０００５】
また、上述の従来技術では、整数値バイアスの確定処理中に測位側が移動することを想定していないが（整数値バイアスの確定処理中は測位側を固定することを前提としているが）、仮に測位側が車両のような移動体であることを想定し、整数値バイアスの確定処理中に測位側が移動することを許容すると、周辺建物やトンネル等の影響をよりサイクルスリップが生じ易くなることが予想される。
【０００６】
そこで、本発明は、上記の問題点に鑑み、サイクルスリップが発生した場合であっても、整数値バイアスの確定処理を一からやり直す必要の無い搬送波位相式測位装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記課題を解決するため、本発明の一局面によれば、移動局及び既知点で測定される衛星信号の搬送波位相の積算値又はその位相差を観測量とし、移動局の位置と搬送波位相の積算値に含まれる整数値バイアス又はその位相差とを状態変数とし、複数の観測周期での衛星データに基づいて前記状態変数を推定して移動局の測位を行う搬送波位相式測位装置であって、
位相積算値のサイクルスリップの有無を検出するサイクルスリップ検出手段を備え、
サイクルスリップが検出された場合は、サイクルスリップの検出された衛星に係る状態変数だけを初期化して測位を継続することを特徴とする、搬送波位相式測位装置が提供される。
【０００８】
本局面において、サイクルスリップが検出された観測周期では、サイクルスリップの検出された衛星に係る状態変数に初期値を用い、他の衛星に係る状態変数に前回周期で導出した値を用いて、前記状態変数の推定を行ってよい。
【０００９】
また、各衛星に係る共分散を含む共分散行列を衛星データに基づいて観測周期毎に更新して前記状態変数を推定する搬送波位相式測位装置においては、
サイクルスリップが検出された場合、前記共分散行列において、該サイクルスリップの検出された衛星に係る共分散に初期値が入れられてよい。
【００１０】
また、サイクルスリップ検出手段は、観測値と理論値との関係値の各観測周期での変化態様に基づいてサイクルスリップを検出するものであってよい。前記状態変数の推定は、カルマンフィルタ又は最小二乗法若しくはその類の最小二乗原理の適用により実行されてよい。
【発明の効果】
【００１１】
本発明によれば、サイクルスリップが発生した場合であっても、整数値バイアスの確定処理を一からやり直す必要の無い搬送波位相式測位装置を得ることができる。
【発明を実施するための最良の形態】
【００１２】
以下、図面を参照して、本発明を実施するための最良の形態の説明を行う。
【００１３】
図１は、本発明に係る搬送波位相式測位システムの構成図である。図１に示すように、ＧＰＳ測位システムは、地球周りを周回するＧＰＳ衛星１０と、地球上の所定位置（既知点）に設置される固定型の基準局２０と、地球上に位置し地球上を移動しうる移動局３０とから構成される。
【００１４】
ＧＰＳ衛星１０は、航法メッセージを地球に向けて常時放送する。航法メッセージには、対応するＧＰＳ衛星１０に関する軌道情報、時計の補正値、電離層の補正係数が含まれている。航法メッセージは、Ｃ／Ａコードにより拡散されＬ１搬送波（周波数：1575.42ＭＨz）に乗せられて、地球に向けて常時放送されている。
【００１５】
尚、現在、２４個のＧＰＳ衛星１０が高度約２０，０００ｋｍの上空で地球を一周しており、各４個のＧＰＳ衛星１０が５５度ずつ傾いた６つの地球周回軌道面に均等に配置されている。従って、天空が開けている場所であれば、地球上のどの場所にいても、常時、少なくとも５個以上のＧＰＳ衛星１０が観測可能である。
【００１６】
図２は、図１の搬送波位相式ＧＰＳ測位システムのより詳細な構成図である。移動局３０は、ＧＰＳ受信機３２を備える。ＧＰＳ受信機３２内には、その周波数がＧＰＳ衛星１０の搬送周波数と一致する発振器（図示せず）が内蔵されている。ＧＰＳ受信機３２は、ＧＰＳアンテナ３２ａを介してＧＰＳ衛星１０から受信した電波（衛星信号）を中間周波数に変換後、ＧＰＳ受信機３２内で発生させたＣ／Ａコードを用いてＣ／Ａコード同期を行い、航法メッセージを取り出す。
【００１７】
また、ＧＰＳ受信機３２は、各ＧＰＳ衛星１０ｉからの搬送波に基づいて、搬送波位相の位相積算値Φｉｕを計測する。尚、位相積算値Φｉｕについて、添え字ｉ（＝０，１，２，・・・）は、各ＧＰＳ衛星１０ｉに割り当てられた番号を示し、添え字ｕは移動局３０側での積算値であることを示す。位相積算値Φｉｕは、次式に示すように、搬送波受信時刻ｔでの発振器の位相Θｉｕ（ｔ）と、ＧＰＳ衛星１０ｉでの衛星信号発生時の搬送波位相Θｉｕ（ｔ−τ）との差として得られる。
Φｉｕ（ｔ）＝Θｉｕ（ｔ）−Θｉｕ（ｔ−τｕ）＋Ｎｉｕ＋εｉｕ（ｔ） 式（１）
ここで、τｕは、ＧＰＳ衛星１０からＧＰＳ受信機３２までのトラベル時間を示し、εｉｕは、ノイズ（誤差）を表わす。尚、位相差の観測開始時点では、ＧＰＳ受信機３２は、搬送波位相の１波長以内の位相を正確に測定できるが、それが何波長目に相当するかを確定できない。このため、位相積算値Φｉｕ（ｔ）には、上式に示すように、不確定な要素として整数値バイアスＮｉｕが導入される。
【００１８】
移動局３０は、また、携帯電話等の通信機３３を備える。通信機３３は、後述する如く、基準局２０側の携帯電話基地局等のような通信施設２３と双方向通信を行うように構成されている。
【００１９】
基準局２０は、ＧＰＳアンテナ２２ａを備えるＧＰＳ受信機２２を有する。ＧＰＳ受信機２２は、移動局３０のＧＰＳ受信機３２と同様に、各ＧＰＳ衛星１０ｉからの搬送波に基づいて、次式に示すように、搬送波受信時刻ｔにおける搬送波位相の積算値Φｉｂ（ｔ）を計測する。
Φｉｂ（ｔ）＝Θｉｂ（ｔ）−Θｉｂ（ｔ−τｂ）＋Ｎｉｂ＋εｉｂ（ｔ） 式（２）
尚、Ｎｉｂは、整数値バイアスを示し、εｉｂは、ノイズ（誤差）を表わす。尚、位相積算値Φｉｂについて、添え字ｂは基準局２０側での積算値であることを示す。基準局２０は、計測した位相積算値Φｉｂを通信施設２３を介して移動局３０に送信する。尚、基準局２０は、所定領域に複数設置されている。各基準局２０と通信施設２３（複数も可）とは、図２に示すように、インターネット等のネットワークを介して接続されてよく、若しくは、各基準局２０毎に通信施設２３が設けられてもよい。前者の構成では、移動局３０は、通信施設２３との間で通信可能な状態である限り、各基準局２０が受信した情報を得ることができる。
【００２０】
図３は、移動局３０に搭載される本発明による搬送波位相式測位装置３４（以下、「測位装置３４」という）の一実施例を示す機能ブロック図である。本実施例の測位装置３４は、演算器４０を中心に構成され、演算器４０には、上述のＧＰＳ受信機３２及び通信機３３に接続されている。演算器４０には、更に、移動局３０に搭載される各種センサ５０が接続される。尚、演算器４０は、ＧＰＳ受信機３２に内蔵されるものであってもよい。また、移動局３０が車両の場合、ＧＰＳ受信機３２及び演算器４０及び／又は通信機３３は、ナビゲーション装置内に実装されてよい。
【００２１】
演算器４０は、マイクロコンピューターから構成されてよく、図３に示すように、衛星位置算出部４２と、動的状態量導入部４４と、整数値バイアス推定部４８とを含む。
【００２２】
衛星位置算出部４２は、ＧＰＳ受信機３２が受信した航法メッセージの軌道情報に基づいて、観測可能な各ＧＰＳ衛星１０ｉの、時刻ｔにおけるワールド座標系での位置（Ｘｉ（ｔ）、Ｙｉ（ｔ）、Ｚｉ（ｔ））を計算する。尚、ＧＰＳ衛星１０は、人工衛星の１つであるので、その運動は、地球重心を含む一定面内（軌道面）に限定される。また、ＧＰＳ衛星１０の軌道は地球重心を１つの焦点とする楕円運動であり、ケプラーの方程式を逐次数値計算することで、軌道面上でのＧＰＳ衛星１０の位置が計算できる。また、搬送波受信時刻ｔでの各ＧＰＳ衛星１０ｉの位置（Ｘｉ（ｔ）、Ｙｉ（ｔ）、Ｚｉ（ｔ））は、ＧＰＳ衛星１０の軌道面とワールド座標系の赤道面が回転関係にあることを考慮して、軌道面上でのＧＰＳ衛星１０の位置を３次元的な回転座標変換することで得られる。尚、ワールド座標系とは、図４に示すように、地球重心を原点として、赤道面内で互いに直交するＸ軸及びY軸、並びに、この両軸に直交するＺ軸により定義される。
【００２３】
動的状態量導入部４４は、所定周期で入力される各種センサ５０の出力信号に基づいて、移動局３０の所定の動的状態量を算出し、既知入力を作成する。例えば、移動局３０が車両の場合、動的状態量導入部４４は、車両の非駆動輪に設定された２つの車輪速センサや、ヨーレートセンサ、左右Ｇ加速度センサ、方位角計等の各種センサ５０の出力信号に基づいて、搬送波受信時刻ｔでの車両の前後速度Ｖｘ（ｔ）及び車両の左右速度Ｖｙ（ｔ）を算出する。
【００２４】
車両の速度ベクトル（Ｖｘ（ｔ）、Ｖｙ（ｔ））は、車体を基準としたボディ座標系（図４参照）に基づいているため、動的状態量導入部４４は、速度ベクトル（Ｖｘ（ｔ）、Ｖｙ（ｔ））を、ローカル座標系を介してワールド座標系へと座標変換する。通常、座標の回転変換は、オイラー角を用いて実現できるが、ここでは、ボディ座標系からローカル座標系への変換に関しては、ロール角及びピッチ角が小さいとしてヨー角ψ（ｔ）のみで実現することとする（但し、ロール角及びピッチ角を考慮することも、ヨー角を無視することも当然に可能である。）。また、ローカル座標系からワールド座標系への変換に関しては、車両位置の経度φ（ｔ）及び緯度λ（ｔ）を用いた変換で実現される。
【００２５】
具体的には、ワールド座標系における車両位置を（Ｘｕ、Ｙｕ、Ｚｕ）とし、車両位置の経度及び緯度を（φ、λ）とすると、ワールド座標系で表わした車両の速度ベクトルｄ／ｄｔ［Ｘｕ、Ｙｕ、Ｚｕ］は、次式の通りである。
d/dt［Ｘｕ、Ｙｕ、Ｚｕ］Ｔ＝rot（φ、λ）・rot（ψ）・［Ｖｘ、Ｖｙ］Ｔ 式（３）
ここで、［］Ｔは行列の転置を意味し、rot（φ、λ）及びrot（ψ）は、次の通りとする。
【００２６】
【数１】
尚、経度φ（ｔ）及び緯度λ（ｔ）は、測量が既に完了している所定地点の既知の経度及び緯度（固定値）であってよい。但し、経度φ（ｔ）及び緯度λ（ｔ）は、単独測位により得られる移動局３０の経度φ（ｔ）及び緯度λ（ｔ）（変動値）であってもよい。また、ヨー角度ψ（ｔ）は、ヨー角速度（ヨーレートセンサの検出信号）を積分することで算出されてよく、若しくは、方位角計を用いて決定されてもよい。
【００２７】
上記式（３）の右辺を入力Ｕ０１、Ｕ０２及びＵ０３と置いて、離散化すれば次のようになる。
Ｘｕ（ｔｎ）＝Ｘｕ（ｔｎ−１）＋ＤＴ・Ｕ０１ 式（４）
Ｙｕ（ｔｎ）＝Ｙｕ（ｔｎ−１）＋ＤＴ・Ｕ０２ 式（５）
Ｚｕ（ｔｎ）＝Ｚｕ（ｔｎ−１）＋ＤＴ・Ｕ０３ 式（６）
よって、最終的な既知入力は、次の通りとなる。
Ｕ＝［ＤＴ・Ｕ０１、ＤＴ・Ｕ０２、ＤＴ・Ｕ０３］Ｔ 式（６−１）
尚、上式において、ＤＴは、サンプル時間（データ更新間隔）であり、ｔｎ＝ｔｎ−１＋ＤＴである。尚、以下では、説明上、サンプル時間ＤＴは、上述のＧＰＳ受信機２２，３２による位相積算値の演算周期(観測周期)と同一であるとする。
【００２８】
整数値バイアス推定部４８では、各ＧＰＳ衛星１０ｉに係る観測データ（特に、移動局３０が通信機３３を介して受信する基準局２０側の位相積算値Φｉｂ、及び、移動局３０側の位相積算値Φｉｕ）に基づいて整数値バイアスが推定される。
【００２９】
具体的には、時刻ｔにおける２つのＧＰＳ衛星１０ｊ、１０ｈ（ｉ＝ｊ、ｈ、但し、ｊ≠ｈ）に関する位相積算値の２重位相差は、次式となる。
Φｊｈｂｕ＝（Φｊｂ（ｔ）−Φｊｕ（ｔ））−（Φｈｂ（ｔ）−Φｈｕ（ｔ）） 式（７）
一方、位相積算値の２重位相差Φｊｈｂｕは、（ＧＰＳ衛星１０ｉとＧＰＳ受信機２２若しくは３２との距離）＝（搬送波の波長Ｌ）×（位相積算値）という物理的な意味合いから、次のようになる。
【００３０】
【数２】
ここで、式（８）における［Ｘｂ（ｔ）、Ｙｂ（ｔ）、Ｚｂ（ｔ）］は、時刻ｔにおける基準局２０のワールド座標系における座標値（既知）であり、［Ｘｕ（ｔ）、Ｙｕ（ｔ）、Ｚｕ（ｔ）］は、時刻ｔにおける移動局３０の座標値（未知）であり、［Ｘｊ（ｔ）、Ｙｊ（ｔ）、Ｚｊ（ｔ）］及び［Ｘｈ（ｔ）、Ｙｈ（ｔ）、Ｚｈ（ｔ）］は、時刻ｔにおける各ＧＰＳ衛星１０ｊ、１０ｈの座標値（衛星位置算出部４２により算出）である。Ｎｊｈｂｕは、整数値バイアスの２重位相差である（即ち、Ｎｊｈｂｕ＝（Ｎｊｂ−Ｎｊｕ）−（Ｎｈｂ−Ｎｈｕ））。尚、時刻ｔは、例えばGPS時刻で同期が取られているものとする。
【００３１】
整数値バイアス推定部４８では、動的状態量導入部４４によって導出された既知入力（上記式（６−１）参照）を用いて、次の状態方程式が設定される。
η（ｔｎ）＝η（ｔｎ−１）＋Ｕ（ｔｎ−１）＋Ｗ（ｔｎ−１） 式（９）
ここで、η（ｔｎ）は、時刻ｔ＝ｔｎでの状態変数を表わし、移動局３０の位置［Ｘｕ（ｔｎ）、Ｙｕ（ｔｎ）、Ｚｕ（ｔｎ）］、及び、整数値バイアスの２重位相差Ｎｉｊｂｕである（但し、ｉ≠ｊ）。ここで、整数値バイアスの２重位相差Ｎｉｊｂｕは少なくとも４個以上必要であり、例えば、５つのGPS衛星１００〜４が観測可能な場合、GPS衛星１００を基準として、η＝［Ｘｕ、Ｙｕ、Ｚｕ、Ｎ０１ｂｕ、Ｎ０２ｂｕ、Ｎ０３ｂｕ、Ｎ０４ｂｕ］であってよい。また、Ｕ及びＷは、それぞれ、上述の既知入力及び外乱（システム雑音：正規性白色雑音）である。
【００３２】
整数値バイアス推定部４８では、次の観測方程式が採用される。
Ｚ（ｔｎ）＝Ｈ（ｔｎ）・η（ｔｎ）＋Ｖ（ｔｎ） 式（１０）
を用いて、ここで、Ｚ及びＶは、それぞれ、観測量及び観測ノイズ（正規性白色雑音）を示す。観測量Ｚは、位相積算値の２重位相差（上記式（７）参照）である。上記式（９）の状態方程式は線形であるが、観測量Ｚは、状態変数Ｘｕ、Ｙｕ及びＺｕに関して非線形であるため、式（８）の各項が状態変数Ｘｕ、Ｙｕ及びＺｕのそれぞれで偏微分され、式（１０）のＨが求められる。
【００３３】
従って、上記式（９）の状態方程式及び上記式（１０）の観測方程式にカルマンフィルタを適用すると、以下の式が得られる。
時間更新として、
η（ｔｎ）（−）＝η（ｔｎ−１）（＋）＋Ｕ（ｔｎ−１） 式（１１）
Ｐ（ｔｎ）（−）＝Ｐ（ｔｎ−１）（＋）＋Ｑ（ｔｎ−１） 式（１２）
また、観測更新として、
Ｋ（ｔｎ）＝Ｐ（ｔｎ）（−）・ＨＴ（ｔｎ）・（Ｈ（ｔｎ）・Ｐ（ｔｎ）（−）・ＨＴ（ｔｎ）＋Ｒ（ｔｎ））−１ 式（１３）
η（ｔｎ）（＋）＝η（ｔｎ）（−）＋Ｋ（ｔｎ）・（Ｚ（ｔｎ）−Ｈ（ｔｎ）・η（ｔｎ）（−）） 式（１４）
Ｐ（ｔｎ）（＋）＝Ｐ（ｔｎ）（−）−Ｋ（ｔｎ）・Ｈ（ｔｎ）・Ｐ（ｔｎ）（−） 式（１５）
ここで、Ｑ，Ｒは、外乱の共分散行列及び観測ノイズの共分散行列をそれぞれ表わす。尚、上記式（１１）及び式（１４）がフィルタ方程式、上記式（１３）がフィルタゲイン、上記式（１２）及び式（１５）が共分散方程式となる。η（ｉ）（−）及びη（ｉ）（＋）は、それぞれ予測値及び推定値であり、Ｐ（ｉ）（−）及びＰ（ｉ）（＋）は、それぞれ予測誤差共分散及び推定誤差共分散であり、η（ｉ）（−）及びη（ｉ）（＋）の推定精度を表している。また、上記式（１４）は、推定値η（ｉ）（＋）は、予測値η（ｉ）（−）に、観測予測誤差にカルマンゲインＫ（修正ゲイン）を乗じたものを補正項として加えることで導出されることを示している。
【００３４】
次に、図５を参照して、本実施例の測位装置３４の特徴的構成・動作について説明する。図５に示す処理ルーチンは、ＧＰＳ受信機２２，３２による位相積算値の演算周期毎（観測周期毎、１観測周期＝サンプル時間ＤＴ）に実行される。以下、時刻ｔ＝ｔｎ（＝ｔｎ−１＋ＤＴ）における処理ルーチンを例として説明する。
【００３５】
先ずステップ１００において、時刻ｔ＝ｔｎにおける衛星情報（各ＧＰＳ衛星１０ｉの位置情報、各ＧＰＳ衛星１０ｉに対する基準局２０側の位相積算値、及び、各ＧＰＳ衛星１０ｉに対する移動局３０側の位相積算値）が取得される。
【００３６】
続くステップ１１０では、取得した位相積算値のサイクルスリップがあったか否かを検査する。ここで、サイクルスリップとは、一エポック内におけるＧＰＳ衛星１０ｉからの衛星信号（電波）の瞬間的な中断（瞬断）に起因して当該エポックに係る位相積算値に整数部の繰り上がり、繰り下がりが生じる場合の他、電波の瞬間的な中断以外の要因（例えば大きなノイズ）に起因して位相積算値に同様の誤差が発生する場合をも含む。
【００３７】
サイクルスリップは、例えば今回周期以前に得た観測値と理論値（又は推定値）との関係（例えば比や残差（イノベーション））に基づいて確率統計解析により検出されてよい。このサイクルスリップの検出方法の一例については後に詳説する。
【００３８】
上記ステップ１１０でサイクルスリップがあったと判断された場合、サイクルスリップのあったＧＰＳ衛星１０ｉに係る誤差共分散行列Ｐの共分散が初期化される（ステップ１２０）。例えば、前回の周期で基準衛星１００以外に５つのＧＰＳ衛星１０１−５が受信されており、今回の周期でＧＰＳ衛星１０４からの受信が途絶えた場合、前回の周期で得られた誤差共分散行列Ｐ（ｔｎ−１）（＋）に対して、図６に示すように、ＧＰＳ衛星１０４に係る共分散が初期化され、他の４つのＧＰＳ衛星１０１−３、５に係る共分散が引き継がれる。即ち、４つのＧＰＳ衛星１０１−３、５に係る共分散は、初期化されることはなく、前回の周期で得られた共分散が引き継がれ、サイクルスリップのあったＧＰＳ衛星１０４に係る共分散のみが、所定の初期値により初期化される。尚、図６には、所定の初期値の１例として枠７０内の数値（対角成分100/3以外は０）が示されている。所定の初期値は、必ずしも所定値である必要は無く、尤もらしい可変値が採用されてよい。
【００３９】
また、本ステップ１２０では、上記式（１１）の時間更新において、誤差共分散行列Ｐ（ｔｎ−１）（＋）と同様に、サイクルスリップのあったＧＰＳ衛星１０ｉに係る状態変数η（ｔｎ−１）（＋）（上記式（１２）参照）に対する初期化が実行される。例えば、図７に示す例では、ＧＰＳ衛星１０４からの衛星信号のサイクルスリップに対応して、ＧＰＳ衛星１０４に関連する状態変数Ｎ０４ｂｕが初期値Ｎ０４ｂｕ’に初期化されている。状態変数の初期値についても同様、必ずしも所定値である必要は無く、尤もらしい可変値が採用されてよい。
【００４０】
尚、上記ステップ１２０において、誤差共分散行列Ｐ（ｔｎ−１）（＋）及び状態変数η（ｔｎ−１）（＋）に対してではなく、誤差共分散行列Ｐ（ｔｎ）（−）及び状態変数η（ｔｎ）（−）に対して初期化が実行されてもよい。
【００４１】
続くステップ１３０では、上記式（１３）乃至式（１５）が適用され、測位算出が実行される。この際、上記ステップ１２０を経由した場合には、上述の初期化が反映されたη（ｔｎ）（−）及びＰ（ｔｎ）（−）が用いられる。本ステップ１３０で得られる今回の周期の誤差共分散行列Ｐ（ｔｎ）（＋）は、状態変数（ｔｎ）（＋）と共に、次回の周期（ｔ＝ｔｎ＋１）における上記ステップ１２０の処理で利用されることになる（上記式（１１）及び（１２）参照）。
【００４２】
続くステップ１４０では、測位解を算出・出力する処理が実行される。上記ステップ１３０から得られる整数値バイアスの推定値は、実数解として求められる。しかし、整数値バイアスは、実際には整数値であるので、求めた実数解に対して最も近い整数解（即ち、波数）を求める。この手法としては、整数値バイアスの無相関化をはかり、整数解の探索空間を狭めて解を特定するＬＡ
ＭＢＤＡ法等が使用されてよい。また、ＧＰＳ受信機２２、３２が、GPS衛星１０から発射されるＬ１波及びＬ２波の双方を受信可能な２周波受信機である場合には、Ｌ１波及びＬ２波のそれぞれに対して上述と同様の推定を同時・並列的に実行し、双方の周期の和（ワイドレーン）を作成して整数値バイアスの整数解の候補を絞り込んでもよい。また、衛星信号に乗せられるＣ／ＡコードやＰコード若しくはその類のＰＲＮコード（擬似雑音符号）を用いて導出される擬似距離の1重又は2重位相差を上述の状態変数に組み込んでもよい。
【００４３】
尚、このようにして整数値バイアスの整数解が確定されると、以後、サイクルスリップが生じない限り、移動局３０の位置は、当該整数値バイアスの整数解を用いた測位により高精度に算出できる。
【００４４】
ところで、本実施例によれば、既知の外部入力Ｕ（k）をカルマンフィルタに入力することで、移動局３０が移動しながらでも状態変数（測位、整数値バイアス）の算出が可能となっている。従って、本実施例による測位装置３４は、車両のような移動体に搭載しても、高精度の測位を実現することができる。
【００４５】
しかしながら、実際の移動局３０の移動中においては、周辺建物やトンネル等の影響により各ＧＰＳ衛星１０ｉからの電波受信状態が安定せず、ＧＰＳ衛星１０ｉからの電波の中断等に起因にしたサイクルスリップが生ずる場合が多くなることが予測される。サイクルスリップが生ずると、位相積算値の整数部の繰り上がり、繰り下がりが不明となるため、従来的には、サイクルスリップが生じたときに、新たなに整数値バイアスの確定処理を一からやり直す構成（即ち、新たな５つのＧＰＳ衛星１０ｉに係る共分散をすべて初期化する構成）をとっていた。しかしながら、かかる構成では、当該初期化による測位精度の一時的な悪化が頻繁に発生するので、車両のような移動体に対して不適である。
【００４６】
これに対して、本実施例によれば、上述の如く、サイクルスリップが発生した場合であっても、それ以前に導出されていた共分散を引き継いで利用しているので、サイクルスリップ発生時に整数値バイアスの確定処理を一からやり直す必要が無く、また、測位精度が大きく悪化することも無い。従って、本実施例によれば、サイクルスリップの発生に対してロバストな測位を実現でき、車両のような移動体に搭載された場合であっても、安定した測位精度を維持することができる。
【００４７】
次に、図７のステップ１１０に関連したサイクルスリップの検出方法の一例について概説する。尚、本発明は、特に以下説明するサイクルスリップの検出方法に限定されるものでなく、如何なる適切な検出方法に対しても適用可能である。
【００４８】
先ず、イノベーションベクトル（理論値と観測値の差）の要素νｔ（ｔ＝1,2,... ,j）について、時刻ｔ＝ｒで期待値E[νｔ]が変化し、変化後の期待値μは未知であると考える。即ち、
【００４９】
【数３】
とする。尚、上述のカルマンフィルタを用いる場合、このイノベーションベクトルの要素νｔは、式（１４）における観測予測誤差（Ｚ（ｔｎ）−Ｈ（ｔｎ）・η（ｔｎ）（−））が用いられてよい。
【００５０】
このとき、H０：変化が起こらなかった場合（ｊ＜ｒ）と、H１：変化が起こった場合（１＜ｒ＜ｊ）の仮説を設定すると、それらの尤度比は、
【００５１】
【数４】
で表される。数４の対数部分をΛｊ（ｒ、μ）とすると、
【００５２】
【数５】
を得る。但し、σ２はイノベーションの分散とする。このとき、変化時刻ｒ及び変化後の未知数μの最尤推定量ｒｊ、μｊ（これらの記号の上には記号ハットが付いている、以下、上付きハットという）は、
【００５３】
【数６】
となる。ここで、
【００５４】
【数７】
として、閾値をλとすると、ｇｊ＞λにおいて仮説H１が採択され、推定変化時刻ｒｊ（上付きハット）を得る。一方、ｇｊ＜λにおいては仮説H０が採択され、仮説H１が棄却される。このようにしてサイクルスリップが推定変化時刻ｒｊで発生したと検出されると、推定変化時刻ｒｊ（周期）で導出された共分散行列や状態変数が、上述のステップ１２０の処理により遡及的に初期化・更新されることになる。尚、推定変化時刻ｒｊまで遡及して初期化をするのではなく、推定変化時刻ｒｊ（周期）の次回周期から同様の初期化が実行されてもよい。
【００５５】
図８は、本実施例による測位装置３４により実現される高い測位精度を実証する試験データである。図８（Ａ）は、本実施例のアルゴリズムが適用された構成（上述のサイクルスリップのあったＧＰＳ衛星１０ｉに係る共分散のみを初期化する構成）による試験結果を示し、横軸にエポック数（観測周期数）を示し、縦軸に測位誤差を示している。図８（Ｂ）は、対照として、サイクルスリップ発生時に誤差共分散行列Ｐ全体を初期化する構成による試験結果を示している。
【００５６】
本実施例のアルゴリズムが適用されない場合には、図８（Ｂ）に示すように、サイクルスリップ発生時に大きな精度誤差が発生し、解のばらつきも大きいのに対して、本実施例のアルゴリズムが適用された場合、図８（Ａ）に示すように、解のばらつきが小さく、良好な精度の解が求めることが可能となった。
【００５７】
以上、本発明の好ましい実施例について詳説したが、本発明は、上述した実施例に制限されることはなく、本発明の範囲を逸脱することなく、上述した実施例に種々の変形及び置換を加えることができる。
【００５８】
例えば、上述した実施例では、上記式（９）の状態方程式及び上記式（１０）の観測方程式にカルマンフィルタを適用するものであったが、本発明は、最小二乗法やベイズ法等の他の推定手法を適用する構成に対しても適用可能である。例えば、最小２乗法の場合、同様に、誤差行列の非対角成分の共分散（及び／又は対角成分の分散）に対して上述と同様の初期化が適用されてよい。
【００５９】
また、上述した実施例では、移動局３０の移動中においても短時間且つ高精度に整数値バイアスを確定できるように、上記式（９）の状態方程式に既知入力Ｕを導入しているが、他の方法により移動局３０の動的な移動が補償されてもよく、或いは、簡易的に当該補償が省略されてもよい（即ち、既知入力Ｕ無し）。
【００６０】
また、上述した実施例では、上述の如く２重位相差を取ることでＧＰＳ受信機２２，３２内での発振器の初期位相、及び、時計誤差等の影響を消去しているが、ＧＰＳ衛星１０の初期位相及びＧＰＳ時計誤差のみを消去できる一重位相差を取る構成であってもよい。また、本実施例では、電離層屈折効果、対流圏屈折効果及びマルチパスの影響を無視しているが、これらを考慮するものであってもよい。
【００６１】
また、上述の説明では、便宜上、ＧＰＳ衛星１００を参照衛星としている場合があるが、移動局３０と基準局２０の位置等に依存して、他のＧＰＳ衛星１０ｉ（＝１，２，・・・）が参照衛星となりえる。
【００６２】
また、上述の説明では、移動局３０の例として車両を挙げたが、移動局３０は、受信機３２及び／又は演算器４０が実装されたホークリフト、ロボットや、受信機３２及び／又は演算器４０を内蔵する携帯電話等の情報端末を含む。
【００６３】
また、上述の実施例では、他の共分散行列（例えば、Q、R）には定数を用いているため、誤差共分散行列Pに対するような初期化処理が不要であるが、定数を採用しない場合、当該他の共分散行列に対しても同様の初期化処理が適用されてよい。
【図面の簡単な説明】
【００６４】
【図１】本発明に係る搬送波位相式ＧＰＳ測位システムの構成図である。
【図２】図１の搬送波位相式ＧＰＳ測位システムのより詳細な構成図である。
【図３】移動局３０に搭載される本発明による測位装置３４の一実施例を示す機能ブロック図である。
【図４】ワールド座標系とローカル座標系との関係、及び、ローカル座標系とボディ座標との関係を示す図である。
【図５】本実施例の測位装置３４により実現される処理のフローチャートである。
【図６】サイクルスリップの発生に伴う誤差共分散行列Ｐの初期化処理の説明図である。
【図７】サイクルスリップの発生に伴う状態変数ηの初期化処理の説明図である。
【図８】本実施例による測位装置３４により実現される高い測位精度を実証する試験データである。
【符号の説明】
【００６５】
１０ ＧＰＳ衛星
２０ 基準局
２２ 基準局側ＧＰＳ受信機
３０ 移動局
３２ 移動局側ＧＰＳ受信機
３４ 搬送波位相式測位装置
４０ 演算器
４２ 衛星位置算出部
４４ 動的状態量導入部
４８ 整数値バイアス推定部
【技術分野】
【０００１】
本発明は、衛星信号の搬送波位相の積算値を測定して移動局の測位を行う搬送波位相式測位装置に関する。
【背景技術】
【０００２】
近年、測量の分野では、搬送波位相によるＧＰＳ測量が広く利用されている。この搬送波位相によるＧＰＳ測量では、基準側の受信機と測位側の受信機とが、複数の衛星から送られる衛星信号を同時に受信し、基準側と測位側とで各衛星信号の搬送波位相の積算値をそれぞれ独立に算出する。この搬送波位相の積算値（以下、単に「位相積算値」という）には、搬送波の波長の整数倍に相当する不確定な要素（以下、「整数値バイアス」という）が含まれているが、この整数値バイアスは、時計誤差等とは異なり、位相積算値の一重位相差や二重位相差を取ることによっても消去することができない。このため、ＧＰＳ測量の分野において、位相積算値の三重位相差を取ることで不確定要素である整数値バイアスを消去することや、整数値バイアスそのものを求める技術が提案されている。
【０００３】
ここで、整数値バイアスを確定する技術として、カルマンフィルタを用いる技術が知られている（例えば、特許文献１参照）。この技術では、測位側の位置と整数値バイアスを状態変数とし、基準側に対する測位側の位相積算値の一重位相差を観測量として、観測を重ねる毎に前記状態変数を更新する追尾フィルタが構成される。また、整数値バイアスを確定するその他の技術として、整数値バイアスを含んだ搬送波の二重位相差を用いて、最小二乗法により所定の条件で二重位相差の整数値バイアスを求める技術が知られている（例えば、特許文献２参照）。
【特許文献１】特開２００４−７７２２８号公報
【特許文献２】特開２００３−９８２４５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
ところで、上述のような従来技術は、少なくとも５つの衛星からの衛星信号が、整数値バイアスの確定処理中に継続して（サイクルスリップ無しに）受信されることを前提としている。しかしながら、実際には、整数値バイアス確定後若しくは確定中にある衛星からの衛星信号（電波）の途切れが生じる場合（サイクルスリップ）がある。このサイクルスリップが起こった場合、測定した位相積算値に整数波数分のジャンプが生じうるので、整数値バイアスを確定するための処理を一からやり直さなければならないという不都合がある。
【０００５】
また、上述の従来技術では、整数値バイアスの確定処理中に測位側が移動することを想定していないが（整数値バイアスの確定処理中は測位側を固定することを前提としているが）、仮に測位側が車両のような移動体であることを想定し、整数値バイアスの確定処理中に測位側が移動することを許容すると、周辺建物やトンネル等の影響をよりサイクルスリップが生じ易くなることが予想される。
【０００６】
そこで、本発明は、上記の問題点に鑑み、サイクルスリップが発生した場合であっても、整数値バイアスの確定処理を一からやり直す必要の無い搬送波位相式測位装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記課題を解決するため、本発明の一局面によれば、移動局及び既知点で測定される衛星信号の搬送波位相の積算値又はその位相差を観測量とし、移動局の位置と搬送波位相の積算値に含まれる整数値バイアス又はその位相差とを状態変数とし、複数の観測周期での衛星データに基づいて前記状態変数を推定して移動局の測位を行う搬送波位相式測位装置であって、
位相積算値のサイクルスリップの有無を検出するサイクルスリップ検出手段を備え、
サイクルスリップが検出された場合は、サイクルスリップの検出された衛星に係る状態変数だけを初期化して測位を継続することを特徴とする、搬送波位相式測位装置が提供される。
【０００８】
本局面において、サイクルスリップが検出された観測周期では、サイクルスリップの検出された衛星に係る状態変数に初期値を用い、他の衛星に係る状態変数に前回周期で導出した値を用いて、前記状態変数の推定を行ってよい。
【０００９】
また、各衛星に係る共分散を含む共分散行列を衛星データに基づいて観測周期毎に更新して前記状態変数を推定する搬送波位相式測位装置においては、
サイクルスリップが検出された場合、前記共分散行列において、該サイクルスリップの検出された衛星に係る共分散に初期値が入れられてよい。
【００１０】
また、サイクルスリップ検出手段は、観測値と理論値との関係値の各観測周期での変化態様に基づいてサイクルスリップを検出するものであってよい。前記状態変数の推定は、カルマンフィルタ又は最小二乗法若しくはその類の最小二乗原理の適用により実行されてよい。
【発明の効果】
【００１１】
本発明によれば、サイクルスリップが発生した場合であっても、整数値バイアスの確定処理を一からやり直す必要の無い搬送波位相式測位装置を得ることができる。
【発明を実施するための最良の形態】
【００１２】
以下、図面を参照して、本発明を実施するための最良の形態の説明を行う。
【００１３】
図１は、本発明に係る搬送波位相式測位システムの構成図である。図１に示すように、ＧＰＳ測位システムは、地球周りを周回するＧＰＳ衛星１０と、地球上の所定位置（既知点）に設置される固定型の基準局２０と、地球上に位置し地球上を移動しうる移動局３０とから構成される。
【００１４】
ＧＰＳ衛星１０は、航法メッセージを地球に向けて常時放送する。航法メッセージには、対応するＧＰＳ衛星１０に関する軌道情報、時計の補正値、電離層の補正係数が含まれている。航法メッセージは、Ｃ／Ａコードにより拡散されＬ１搬送波（周波数：1575.42ＭＨz）に乗せられて、地球に向けて常時放送されている。
【００１５】
尚、現在、２４個のＧＰＳ衛星１０が高度約２０，０００ｋｍの上空で地球を一周しており、各４個のＧＰＳ衛星１０が５５度ずつ傾いた６つの地球周回軌道面に均等に配置されている。従って、天空が開けている場所であれば、地球上のどの場所にいても、常時、少なくとも５個以上のＧＰＳ衛星１０が観測可能である。
【００１６】
図２は、図１の搬送波位相式ＧＰＳ測位システムのより詳細な構成図である。移動局３０は、ＧＰＳ受信機３２を備える。ＧＰＳ受信機３２内には、その周波数がＧＰＳ衛星１０の搬送周波数と一致する発振器（図示せず）が内蔵されている。ＧＰＳ受信機３２は、ＧＰＳアンテナ３２ａを介してＧＰＳ衛星１０から受信した電波（衛星信号）を中間周波数に変換後、ＧＰＳ受信機３２内で発生させたＣ／Ａコードを用いてＣ／Ａコード同期を行い、航法メッセージを取り出す。
【００１７】
また、ＧＰＳ受信機３２は、各ＧＰＳ衛星１０ｉからの搬送波に基づいて、搬送波位相の位相積算値Φｉｕを計測する。尚、位相積算値Φｉｕについて、添え字ｉ（＝０，１，２，・・・）は、各ＧＰＳ衛星１０ｉに割り当てられた番号を示し、添え字ｕは移動局３０側での積算値であることを示す。位相積算値Φｉｕは、次式に示すように、搬送波受信時刻ｔでの発振器の位相Θｉｕ（ｔ）と、ＧＰＳ衛星１０ｉでの衛星信号発生時の搬送波位相Θｉｕ（ｔ−τ）との差として得られる。
Φｉｕ（ｔ）＝Θｉｕ（ｔ）−Θｉｕ（ｔ−τｕ）＋Ｎｉｕ＋εｉｕ（ｔ） 式（１）
ここで、τｕは、ＧＰＳ衛星１０からＧＰＳ受信機３２までのトラベル時間を示し、εｉｕは、ノイズ（誤差）を表わす。尚、位相差の観測開始時点では、ＧＰＳ受信機３２は、搬送波位相の１波長以内の位相を正確に測定できるが、それが何波長目に相当するかを確定できない。このため、位相積算値Φｉｕ（ｔ）には、上式に示すように、不確定な要素として整数値バイアスＮｉｕが導入される。
【００１８】
移動局３０は、また、携帯電話等の通信機３３を備える。通信機３３は、後述する如く、基準局２０側の携帯電話基地局等のような通信施設２３と双方向通信を行うように構成されている。
【００１９】
基準局２０は、ＧＰＳアンテナ２２ａを備えるＧＰＳ受信機２２を有する。ＧＰＳ受信機２２は、移動局３０のＧＰＳ受信機３２と同様に、各ＧＰＳ衛星１０ｉからの搬送波に基づいて、次式に示すように、搬送波受信時刻ｔにおける搬送波位相の積算値Φｉｂ（ｔ）を計測する。
Φｉｂ（ｔ）＝Θｉｂ（ｔ）−Θｉｂ（ｔ−τｂ）＋Ｎｉｂ＋εｉｂ（ｔ） 式（２）
尚、Ｎｉｂは、整数値バイアスを示し、εｉｂは、ノイズ（誤差）を表わす。尚、位相積算値Φｉｂについて、添え字ｂは基準局２０側での積算値であることを示す。基準局２０は、計測した位相積算値Φｉｂを通信施設２３を介して移動局３０に送信する。尚、基準局２０は、所定領域に複数設置されている。各基準局２０と通信施設２３（複数も可）とは、図２に示すように、インターネット等のネットワークを介して接続されてよく、若しくは、各基準局２０毎に通信施設２３が設けられてもよい。前者の構成では、移動局３０は、通信施設２３との間で通信可能な状態である限り、各基準局２０が受信した情報を得ることができる。
【００２０】
図３は、移動局３０に搭載される本発明による搬送波位相式測位装置３４（以下、「測位装置３４」という）の一実施例を示す機能ブロック図である。本実施例の測位装置３４は、演算器４０を中心に構成され、演算器４０には、上述のＧＰＳ受信機３２及び通信機３３に接続されている。演算器４０には、更に、移動局３０に搭載される各種センサ５０が接続される。尚、演算器４０は、ＧＰＳ受信機３２に内蔵されるものであってもよい。また、移動局３０が車両の場合、ＧＰＳ受信機３２及び演算器４０及び／又は通信機３３は、ナビゲーション装置内に実装されてよい。
【００２１】
演算器４０は、マイクロコンピューターから構成されてよく、図３に示すように、衛星位置算出部４２と、動的状態量導入部４４と、整数値バイアス推定部４８とを含む。
【００２２】
衛星位置算出部４２は、ＧＰＳ受信機３２が受信した航法メッセージの軌道情報に基づいて、観測可能な各ＧＰＳ衛星１０ｉの、時刻ｔにおけるワールド座標系での位置（Ｘｉ（ｔ）、Ｙｉ（ｔ）、Ｚｉ（ｔ））を計算する。尚、ＧＰＳ衛星１０は、人工衛星の１つであるので、その運動は、地球重心を含む一定面内（軌道面）に限定される。また、ＧＰＳ衛星１０の軌道は地球重心を１つの焦点とする楕円運動であり、ケプラーの方程式を逐次数値計算することで、軌道面上でのＧＰＳ衛星１０の位置が計算できる。また、搬送波受信時刻ｔでの各ＧＰＳ衛星１０ｉの位置（Ｘｉ（ｔ）、Ｙｉ（ｔ）、Ｚｉ（ｔ））は、ＧＰＳ衛星１０の軌道面とワールド座標系の赤道面が回転関係にあることを考慮して、軌道面上でのＧＰＳ衛星１０の位置を３次元的な回転座標変換することで得られる。尚、ワールド座標系とは、図４に示すように、地球重心を原点として、赤道面内で互いに直交するＸ軸及びY軸、並びに、この両軸に直交するＺ軸により定義される。
【００２３】
動的状態量導入部４４は、所定周期で入力される各種センサ５０の出力信号に基づいて、移動局３０の所定の動的状態量を算出し、既知入力を作成する。例えば、移動局３０が車両の場合、動的状態量導入部４４は、車両の非駆動輪に設定された２つの車輪速センサや、ヨーレートセンサ、左右Ｇ加速度センサ、方位角計等の各種センサ５０の出力信号に基づいて、搬送波受信時刻ｔでの車両の前後速度Ｖｘ（ｔ）及び車両の左右速度Ｖｙ（ｔ）を算出する。
【００２４】
車両の速度ベクトル（Ｖｘ（ｔ）、Ｖｙ（ｔ））は、車体を基準としたボディ座標系（図４参照）に基づいているため、動的状態量導入部４４は、速度ベクトル（Ｖｘ（ｔ）、Ｖｙ（ｔ））を、ローカル座標系を介してワールド座標系へと座標変換する。通常、座標の回転変換は、オイラー角を用いて実現できるが、ここでは、ボディ座標系からローカル座標系への変換に関しては、ロール角及びピッチ角が小さいとしてヨー角ψ（ｔ）のみで実現することとする（但し、ロール角及びピッチ角を考慮することも、ヨー角を無視することも当然に可能である。）。また、ローカル座標系からワールド座標系への変換に関しては、車両位置の経度φ（ｔ）及び緯度λ（ｔ）を用いた変換で実現される。
【００２５】
具体的には、ワールド座標系における車両位置を（Ｘｕ、Ｙｕ、Ｚｕ）とし、車両位置の経度及び緯度を（φ、λ）とすると、ワールド座標系で表わした車両の速度ベクトルｄ／ｄｔ［Ｘｕ、Ｙｕ、Ｚｕ］は、次式の通りである。
d/dt［Ｘｕ、Ｙｕ、Ｚｕ］Ｔ＝rot（φ、λ）・rot（ψ）・［Ｖｘ、Ｖｙ］Ｔ 式（３）
ここで、［］Ｔは行列の転置を意味し、rot（φ、λ）及びrot（ψ）は、次の通りとする。
【００２６】
【数１】
尚、経度φ（ｔ）及び緯度λ（ｔ）は、測量が既に完了している所定地点の既知の経度及び緯度（固定値）であってよい。但し、経度φ（ｔ）及び緯度λ（ｔ）は、単独測位により得られる移動局３０の経度φ（ｔ）及び緯度λ（ｔ）（変動値）であってもよい。また、ヨー角度ψ（ｔ）は、ヨー角速度（ヨーレートセンサの検出信号）を積分することで算出されてよく、若しくは、方位角計を用いて決定されてもよい。
【００２７】
上記式（３）の右辺を入力Ｕ０１、Ｕ０２及びＵ０３と置いて、離散化すれば次のようになる。
Ｘｕ（ｔｎ）＝Ｘｕ（ｔｎ−１）＋ＤＴ・Ｕ０１ 式（４）
Ｙｕ（ｔｎ）＝Ｙｕ（ｔｎ−１）＋ＤＴ・Ｕ０２ 式（５）
Ｚｕ（ｔｎ）＝Ｚｕ（ｔｎ−１）＋ＤＴ・Ｕ０３ 式（６）
よって、最終的な既知入力は、次の通りとなる。
Ｕ＝［ＤＴ・Ｕ０１、ＤＴ・Ｕ０２、ＤＴ・Ｕ０３］Ｔ 式（６−１）
尚、上式において、ＤＴは、サンプル時間（データ更新間隔）であり、ｔｎ＝ｔｎ−１＋ＤＴである。尚、以下では、説明上、サンプル時間ＤＴは、上述のＧＰＳ受信機２２，３２による位相積算値の演算周期(観測周期)と同一であるとする。
【００２８】
整数値バイアス推定部４８では、各ＧＰＳ衛星１０ｉに係る観測データ（特に、移動局３０が通信機３３を介して受信する基準局２０側の位相積算値Φｉｂ、及び、移動局３０側の位相積算値Φｉｕ）に基づいて整数値バイアスが推定される。
【００２９】
具体的には、時刻ｔにおける２つのＧＰＳ衛星１０ｊ、１０ｈ（ｉ＝ｊ、ｈ、但し、ｊ≠ｈ）に関する位相積算値の２重位相差は、次式となる。
Φｊｈｂｕ＝（Φｊｂ（ｔ）−Φｊｕ（ｔ））−（Φｈｂ（ｔ）−Φｈｕ（ｔ）） 式（７）
一方、位相積算値の２重位相差Φｊｈｂｕは、（ＧＰＳ衛星１０ｉとＧＰＳ受信機２２若しくは３２との距離）＝（搬送波の波長Ｌ）×（位相積算値）という物理的な意味合いから、次のようになる。
【００３０】
【数２】
ここで、式（８）における［Ｘｂ（ｔ）、Ｙｂ（ｔ）、Ｚｂ（ｔ）］は、時刻ｔにおける基準局２０のワールド座標系における座標値（既知）であり、［Ｘｕ（ｔ）、Ｙｕ（ｔ）、Ｚｕ（ｔ）］は、時刻ｔにおける移動局３０の座標値（未知）であり、［Ｘｊ（ｔ）、Ｙｊ（ｔ）、Ｚｊ（ｔ）］及び［Ｘｈ（ｔ）、Ｙｈ（ｔ）、Ｚｈ（ｔ）］は、時刻ｔにおける各ＧＰＳ衛星１０ｊ、１０ｈの座標値（衛星位置算出部４２により算出）である。Ｎｊｈｂｕは、整数値バイアスの２重位相差である（即ち、Ｎｊｈｂｕ＝（Ｎｊｂ−Ｎｊｕ）−（Ｎｈｂ−Ｎｈｕ））。尚、時刻ｔは、例えばGPS時刻で同期が取られているものとする。
【００３１】
整数値バイアス推定部４８では、動的状態量導入部４４によって導出された既知入力（上記式（６−１）参照）を用いて、次の状態方程式が設定される。
η（ｔｎ）＝η（ｔｎ−１）＋Ｕ（ｔｎ−１）＋Ｗ（ｔｎ−１） 式（９）
ここで、η（ｔｎ）は、時刻ｔ＝ｔｎでの状態変数を表わし、移動局３０の位置［Ｘｕ（ｔｎ）、Ｙｕ（ｔｎ）、Ｚｕ（ｔｎ）］、及び、整数値バイアスの２重位相差Ｎｉｊｂｕである（但し、ｉ≠ｊ）。ここで、整数値バイアスの２重位相差Ｎｉｊｂｕは少なくとも４個以上必要であり、例えば、５つのGPS衛星１００〜４が観測可能な場合、GPS衛星１００を基準として、η＝［Ｘｕ、Ｙｕ、Ｚｕ、Ｎ０１ｂｕ、Ｎ０２ｂｕ、Ｎ０３ｂｕ、Ｎ０４ｂｕ］であってよい。また、Ｕ及びＷは、それぞれ、上述の既知入力及び外乱（システム雑音：正規性白色雑音）である。
【００３２】
整数値バイアス推定部４８では、次の観測方程式が採用される。
Ｚ（ｔｎ）＝Ｈ（ｔｎ）・η（ｔｎ）＋Ｖ（ｔｎ） 式（１０）
を用いて、ここで、Ｚ及びＶは、それぞれ、観測量及び観測ノイズ（正規性白色雑音）を示す。観測量Ｚは、位相積算値の２重位相差（上記式（７）参照）である。上記式（９）の状態方程式は線形であるが、観測量Ｚは、状態変数Ｘｕ、Ｙｕ及びＺｕに関して非線形であるため、式（８）の各項が状態変数Ｘｕ、Ｙｕ及びＺｕのそれぞれで偏微分され、式（１０）のＨが求められる。
【００３３】
従って、上記式（９）の状態方程式及び上記式（１０）の観測方程式にカルマンフィルタを適用すると、以下の式が得られる。
時間更新として、
η（ｔｎ）（−）＝η（ｔｎ−１）（＋）＋Ｕ（ｔｎ−１） 式（１１）
Ｐ（ｔｎ）（−）＝Ｐ（ｔｎ−１）（＋）＋Ｑ（ｔｎ−１） 式（１２）
また、観測更新として、
Ｋ（ｔｎ）＝Ｐ（ｔｎ）（−）・ＨＴ（ｔｎ）・（Ｈ（ｔｎ）・Ｐ（ｔｎ）（−）・ＨＴ（ｔｎ）＋Ｒ（ｔｎ））−１ 式（１３）
η（ｔｎ）（＋）＝η（ｔｎ）（−）＋Ｋ（ｔｎ）・（Ｚ（ｔｎ）−Ｈ（ｔｎ）・η（ｔｎ）（−）） 式（１４）
Ｐ（ｔｎ）（＋）＝Ｐ（ｔｎ）（−）−Ｋ（ｔｎ）・Ｈ（ｔｎ）・Ｐ（ｔｎ）（−） 式（１５）
ここで、Ｑ，Ｒは、外乱の共分散行列及び観測ノイズの共分散行列をそれぞれ表わす。尚、上記式（１１）及び式（１４）がフィルタ方程式、上記式（１３）がフィルタゲイン、上記式（１２）及び式（１５）が共分散方程式となる。η（ｉ）（−）及びη（ｉ）（＋）は、それぞれ予測値及び推定値であり、Ｐ（ｉ）（−）及びＰ（ｉ）（＋）は、それぞれ予測誤差共分散及び推定誤差共分散であり、η（ｉ）（−）及びη（ｉ）（＋）の推定精度を表している。また、上記式（１４）は、推定値η（ｉ）（＋）は、予測値η（ｉ）（−）に、観測予測誤差にカルマンゲインＫ（修正ゲイン）を乗じたものを補正項として加えることで導出されることを示している。
【００３４】
次に、図５を参照して、本実施例の測位装置３４の特徴的構成・動作について説明する。図５に示す処理ルーチンは、ＧＰＳ受信機２２，３２による位相積算値の演算周期毎（観測周期毎、１観測周期＝サンプル時間ＤＴ）に実行される。以下、時刻ｔ＝ｔｎ（＝ｔｎ−１＋ＤＴ）における処理ルーチンを例として説明する。
【００３５】
先ずステップ１００において、時刻ｔ＝ｔｎにおける衛星情報（各ＧＰＳ衛星１０ｉの位置情報、各ＧＰＳ衛星１０ｉに対する基準局２０側の位相積算値、及び、各ＧＰＳ衛星１０ｉに対する移動局３０側の位相積算値）が取得される。
【００３６】
続くステップ１１０では、取得した位相積算値のサイクルスリップがあったか否かを検査する。ここで、サイクルスリップとは、一エポック内におけるＧＰＳ衛星１０ｉからの衛星信号（電波）の瞬間的な中断（瞬断）に起因して当該エポックに係る位相積算値に整数部の繰り上がり、繰り下がりが生じる場合の他、電波の瞬間的な中断以外の要因（例えば大きなノイズ）に起因して位相積算値に同様の誤差が発生する場合をも含む。
【００３７】
サイクルスリップは、例えば今回周期以前に得た観測値と理論値（又は推定値）との関係（例えば比や残差（イノベーション））に基づいて確率統計解析により検出されてよい。このサイクルスリップの検出方法の一例については後に詳説する。
【００３８】
上記ステップ１１０でサイクルスリップがあったと判断された場合、サイクルスリップのあったＧＰＳ衛星１０ｉに係る誤差共分散行列Ｐの共分散が初期化される（ステップ１２０）。例えば、前回の周期で基準衛星１００以外に５つのＧＰＳ衛星１０１−５が受信されており、今回の周期でＧＰＳ衛星１０４からの受信が途絶えた場合、前回の周期で得られた誤差共分散行列Ｐ（ｔｎ−１）（＋）に対して、図６に示すように、ＧＰＳ衛星１０４に係る共分散が初期化され、他の４つのＧＰＳ衛星１０１−３、５に係る共分散が引き継がれる。即ち、４つのＧＰＳ衛星１０１−３、５に係る共分散は、初期化されることはなく、前回の周期で得られた共分散が引き継がれ、サイクルスリップのあったＧＰＳ衛星１０４に係る共分散のみが、所定の初期値により初期化される。尚、図６には、所定の初期値の１例として枠７０内の数値（対角成分100/3以外は０）が示されている。所定の初期値は、必ずしも所定値である必要は無く、尤もらしい可変値が採用されてよい。
【００３９】
また、本ステップ１２０では、上記式（１１）の時間更新において、誤差共分散行列Ｐ（ｔｎ−１）（＋）と同様に、サイクルスリップのあったＧＰＳ衛星１０ｉに係る状態変数η（ｔｎ−１）（＋）（上記式（１２）参照）に対する初期化が実行される。例えば、図７に示す例では、ＧＰＳ衛星１０４からの衛星信号のサイクルスリップに対応して、ＧＰＳ衛星１０４に関連する状態変数Ｎ０４ｂｕが初期値Ｎ０４ｂｕ’に初期化されている。状態変数の初期値についても同様、必ずしも所定値である必要は無く、尤もらしい可変値が採用されてよい。
【００４０】
尚、上記ステップ１２０において、誤差共分散行列Ｐ（ｔｎ−１）（＋）及び状態変数η（ｔｎ−１）（＋）に対してではなく、誤差共分散行列Ｐ（ｔｎ）（−）及び状態変数η（ｔｎ）（−）に対して初期化が実行されてもよい。
【００４１】
続くステップ１３０では、上記式（１３）乃至式（１５）が適用され、測位算出が実行される。この際、上記ステップ１２０を経由した場合には、上述の初期化が反映されたη（ｔｎ）（−）及びＰ（ｔｎ）（−）が用いられる。本ステップ１３０で得られる今回の周期の誤差共分散行列Ｐ（ｔｎ）（＋）は、状態変数（ｔｎ）（＋）と共に、次回の周期（ｔ＝ｔｎ＋１）における上記ステップ１２０の処理で利用されることになる（上記式（１１）及び（１２）参照）。
【００４２】
続くステップ１４０では、測位解を算出・出力する処理が実行される。上記ステップ１３０から得られる整数値バイアスの推定値は、実数解として求められる。しかし、整数値バイアスは、実際には整数値であるので、求めた実数解に対して最も近い整数解（即ち、波数）を求める。この手法としては、整数値バイアスの無相関化をはかり、整数解の探索空間を狭めて解を特定するＬＡ
ＭＢＤＡ法等が使用されてよい。また、ＧＰＳ受信機２２、３２が、GPS衛星１０から発射されるＬ１波及びＬ２波の双方を受信可能な２周波受信機である場合には、Ｌ１波及びＬ２波のそれぞれに対して上述と同様の推定を同時・並列的に実行し、双方の周期の和（ワイドレーン）を作成して整数値バイアスの整数解の候補を絞り込んでもよい。また、衛星信号に乗せられるＣ／ＡコードやＰコード若しくはその類のＰＲＮコード（擬似雑音符号）を用いて導出される擬似距離の1重又は2重位相差を上述の状態変数に組み込んでもよい。
【００４３】
尚、このようにして整数値バイアスの整数解が確定されると、以後、サイクルスリップが生じない限り、移動局３０の位置は、当該整数値バイアスの整数解を用いた測位により高精度に算出できる。
【００４４】
ところで、本実施例によれば、既知の外部入力Ｕ（k）をカルマンフィルタに入力することで、移動局３０が移動しながらでも状態変数（測位、整数値バイアス）の算出が可能となっている。従って、本実施例による測位装置３４は、車両のような移動体に搭載しても、高精度の測位を実現することができる。
【００４５】
しかしながら、実際の移動局３０の移動中においては、周辺建物やトンネル等の影響により各ＧＰＳ衛星１０ｉからの電波受信状態が安定せず、ＧＰＳ衛星１０ｉからの電波の中断等に起因にしたサイクルスリップが生ずる場合が多くなることが予測される。サイクルスリップが生ずると、位相積算値の整数部の繰り上がり、繰り下がりが不明となるため、従来的には、サイクルスリップが生じたときに、新たなに整数値バイアスの確定処理を一からやり直す構成（即ち、新たな５つのＧＰＳ衛星１０ｉに係る共分散をすべて初期化する構成）をとっていた。しかしながら、かかる構成では、当該初期化による測位精度の一時的な悪化が頻繁に発生するので、車両のような移動体に対して不適である。
【００４６】
これに対して、本実施例によれば、上述の如く、サイクルスリップが発生した場合であっても、それ以前に導出されていた共分散を引き継いで利用しているので、サイクルスリップ発生時に整数値バイアスの確定処理を一からやり直す必要が無く、また、測位精度が大きく悪化することも無い。従って、本実施例によれば、サイクルスリップの発生に対してロバストな測位を実現でき、車両のような移動体に搭載された場合であっても、安定した測位精度を維持することができる。
【００４７】
次に、図７のステップ１１０に関連したサイクルスリップの検出方法の一例について概説する。尚、本発明は、特に以下説明するサイクルスリップの検出方法に限定されるものでなく、如何なる適切な検出方法に対しても適用可能である。
【００４８】
先ず、イノベーションベクトル（理論値と観測値の差）の要素νｔ（ｔ＝1,2,... ,j）について、時刻ｔ＝ｒで期待値E[νｔ]が変化し、変化後の期待値μは未知であると考える。即ち、
【００４９】
【数３】
とする。尚、上述のカルマンフィルタを用いる場合、このイノベーションベクトルの要素νｔは、式（１４）における観測予測誤差（Ｚ（ｔｎ）−Ｈ（ｔｎ）・η（ｔｎ）（−））が用いられてよい。
【００５０】
このとき、H０：変化が起こらなかった場合（ｊ＜ｒ）と、H１：変化が起こった場合（１＜ｒ＜ｊ）の仮説を設定すると、それらの尤度比は、
【００５１】
【数４】
で表される。数４の対数部分をΛｊ（ｒ、μ）とすると、
【００５２】
【数５】
を得る。但し、σ２はイノベーションの分散とする。このとき、変化時刻ｒ及び変化後の未知数μの最尤推定量ｒｊ、μｊ（これらの記号の上には記号ハットが付いている、以下、上付きハットという）は、
【００５３】
【数６】
となる。ここで、
【００５４】
【数７】
として、閾値をλとすると、ｇｊ＞λにおいて仮説H１が採択され、推定変化時刻ｒｊ（上付きハット）を得る。一方、ｇｊ＜λにおいては仮説H０が採択され、仮説H１が棄却される。このようにしてサイクルスリップが推定変化時刻ｒｊで発生したと検出されると、推定変化時刻ｒｊ（周期）で導出された共分散行列や状態変数が、上述のステップ１２０の処理により遡及的に初期化・更新されることになる。尚、推定変化時刻ｒｊまで遡及して初期化をするのではなく、推定変化時刻ｒｊ（周期）の次回周期から同様の初期化が実行されてもよい。
【００５５】
図８は、本実施例による測位装置３４により実現される高い測位精度を実証する試験データである。図８（Ａ）は、本実施例のアルゴリズムが適用された構成（上述のサイクルスリップのあったＧＰＳ衛星１０ｉに係る共分散のみを初期化する構成）による試験結果を示し、横軸にエポック数（観測周期数）を示し、縦軸に測位誤差を示している。図８（Ｂ）は、対照として、サイクルスリップ発生時に誤差共分散行列Ｐ全体を初期化する構成による試験結果を示している。
【００５６】
本実施例のアルゴリズムが適用されない場合には、図８（Ｂ）に示すように、サイクルスリップ発生時に大きな精度誤差が発生し、解のばらつきも大きいのに対して、本実施例のアルゴリズムが適用された場合、図８（Ａ）に示すように、解のばらつきが小さく、良好な精度の解が求めることが可能となった。
【００５７】
以上、本発明の好ましい実施例について詳説したが、本発明は、上述した実施例に制限されることはなく、本発明の範囲を逸脱することなく、上述した実施例に種々の変形及び置換を加えることができる。
【００５８】
例えば、上述した実施例では、上記式（９）の状態方程式及び上記式（１０）の観測方程式にカルマンフィルタを適用するものであったが、本発明は、最小二乗法やベイズ法等の他の推定手法を適用する構成に対しても適用可能である。例えば、最小２乗法の場合、同様に、誤差行列の非対角成分の共分散（及び／又は対角成分の分散）に対して上述と同様の初期化が適用されてよい。
【００５９】
また、上述した実施例では、移動局３０の移動中においても短時間且つ高精度に整数値バイアスを確定できるように、上記式（９）の状態方程式に既知入力Ｕを導入しているが、他の方法により移動局３０の動的な移動が補償されてもよく、或いは、簡易的に当該補償が省略されてもよい（即ち、既知入力Ｕ無し）。
【００６０】
また、上述した実施例では、上述の如く２重位相差を取ることでＧＰＳ受信機２２，３２内での発振器の初期位相、及び、時計誤差等の影響を消去しているが、ＧＰＳ衛星１０の初期位相及びＧＰＳ時計誤差のみを消去できる一重位相差を取る構成であってもよい。また、本実施例では、電離層屈折効果、対流圏屈折効果及びマルチパスの影響を無視しているが、これらを考慮するものであってもよい。
【００６１】
また、上述の説明では、便宜上、ＧＰＳ衛星１００を参照衛星としている場合があるが、移動局３０と基準局２０の位置等に依存して、他のＧＰＳ衛星１０ｉ（＝１，２，・・・）が参照衛星となりえる。
【００６２】
また、上述の説明では、移動局３０の例として車両を挙げたが、移動局３０は、受信機３２及び／又は演算器４０が実装されたホークリフト、ロボットや、受信機３２及び／又は演算器４０を内蔵する携帯電話等の情報端末を含む。
【００６３】
また、上述の実施例では、他の共分散行列（例えば、Q、R）には定数を用いているため、誤差共分散行列Pに対するような初期化処理が不要であるが、定数を採用しない場合、当該他の共分散行列に対しても同様の初期化処理が適用されてよい。
【図面の簡単な説明】
【００６４】
【図１】本発明に係る搬送波位相式ＧＰＳ測位システムの構成図である。
【図２】図１の搬送波位相式ＧＰＳ測位システムのより詳細な構成図である。
【図３】移動局３０に搭載される本発明による測位装置３４の一実施例を示す機能ブロック図である。
【図４】ワールド座標系とローカル座標系との関係、及び、ローカル座標系とボディ座標との関係を示す図である。
【図５】本実施例の測位装置３４により実現される処理のフローチャートである。
【図６】サイクルスリップの発生に伴う誤差共分散行列Ｐの初期化処理の説明図である。
【図７】サイクルスリップの発生に伴う状態変数ηの初期化処理の説明図である。
【図８】本実施例による測位装置３４により実現される高い測位精度を実証する試験データである。
【符号の説明】
【００６５】
１０ ＧＰＳ衛星
２０ 基準局
２２ 基準局側ＧＰＳ受信機
３０ 移動局
３２ 移動局側ＧＰＳ受信機
３４ 搬送波位相式測位装置
４０ 演算器
４２ 衛星位置算出部
４４ 動的状態量導入部
４８ 整数値バイアス推定部
【特許請求の範囲】
【請求項１】
移動局及び既知点で測定される衛星信号の搬送波位相の積算値又はその位相差を観測量とし、移動局の位置と搬送波位相の積算値に含まれる整数値バイアス又はその位相差とを状態変数とし、複数の観測周期での衛星データに基づいて前記状態変数を推定して移動局の測位を行う搬送波位相式測位装置であって、
位相積算値のサイクルスリップの有無を検出するサイクルスリップ検出手段を備え、
サイクルスリップが検出された場合は、サイクルスリップの検出された衛星に係る状態変数だけを初期化して測位を継続することを特徴とする、搬送波位相式測位装置。
【請求項２】
サイクルスリップが検出された観測周期では、サイクルスリップの検出された衛星に係る状態変数に初期値を用い、他の衛星に係る状態変数に前回周期で導出した値を用いて、前記状態変数の推定を行う、請求項1に記載の搬送波位相式測位装置。
【請求項３】
各衛星に係る共分散を含む共分散行列を衛星データに基づいて観測周期毎に更新して前記状態変数を推定する請求項1に記載の搬送波位相式測位装置において、
サイクルスリップが検出された場合、前記共分散行列において、該サイクルスリップの検出された衛星に係る共分散に初期値が入れられる、搬送波位相式測位装置。
【請求項４】
サイクルスリップ検出手段は、観測値と理論値との関係値の各観測周期での変化態様に基づいてサイクルスリップを検出する、請求項1に記載の搬送波位相式測位装置。
【請求項５】
前記状態変数の推定は、カルマンフィルタ又は最小二乗法若しくはその類の最小二乗原理の適用により実行される、請求項１に記載の搬送波位相式測位装置。
【請求項１】
移動局及び既知点で測定される衛星信号の搬送波位相の積算値又はその位相差を観測量とし、移動局の位置と搬送波位相の積算値に含まれる整数値バイアス又はその位相差とを状態変数とし、複数の観測周期での衛星データに基づいて前記状態変数を推定して移動局の測位を行う搬送波位相式測位装置であって、
位相積算値のサイクルスリップの有無を検出するサイクルスリップ検出手段を備え、
サイクルスリップが検出された場合は、サイクルスリップの検出された衛星に係る状態変数だけを初期化して測位を継続することを特徴とする、搬送波位相式測位装置。
【請求項２】
サイクルスリップが検出された観測周期では、サイクルスリップの検出された衛星に係る状態変数に初期値を用い、他の衛星に係る状態変数に前回周期で導出した値を用いて、前記状態変数の推定を行う、請求項1に記載の搬送波位相式測位装置。
【請求項３】
各衛星に係る共分散を含む共分散行列を衛星データに基づいて観測周期毎に更新して前記状態変数を推定する請求項1に記載の搬送波位相式測位装置において、
サイクルスリップが検出された場合、前記共分散行列において、該サイクルスリップの検出された衛星に係る共分散に初期値が入れられる、搬送波位相式測位装置。
【請求項４】
サイクルスリップ検出手段は、観測値と理論値との関係値の各観測周期での変化態様に基づいてサイクルスリップを検出する、請求項1に記載の搬送波位相式測位装置。
【請求項５】
前記状態変数の推定は、カルマンフィルタ又は最小二乗法若しくはその類の最小二乗原理の適用により実行される、請求項１に記載の搬送波位相式測位装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【公開番号】特開２００８−３９６９１（Ｐ２００８−３９６９１Ａ）
【公開日】平成２０年２月２１日（２００８．２．２１）
【国際特許分類】
【出願番号】特願２００６−２１７３１５（Ｐ２００６−２１７３１５）
【出願日】平成１８年８月９日（２００６．８．９）
【出願人】（０００００３２０７）トヨタ自動車株式会社 (59,920)
【出願人】（５９３００６６３０）学校法人立命館 (359)
【Ｆターム（参考）】
【公開日】平成２０年２月２１日（２００８．２．２１）
【国際特許分類】
【出願日】平成１８年８月９日（２００６．８．９）
【出願人】（０００００３２０７）トヨタ自動車株式会社 (59,920)
【出願人】（５９３００６６３０）学校法人立命館 (359)
【Ｆターム（参考）】
[ Back to top ]