(57)【要約】【課題】厳しい受信環境下での乱れた走行軌跡を滑らかにし、測位精度の劣化を抑える。【解決手段】測位装置１００は、測位信号を送信する人工衛星から前記測位信号を受信して所定の情報を抽出し、前記所定の情報を使用して測位計算を行って緯度経度座標値を出力する測位計算部１０１と、初期値と初期誤差共分散行列と駆動雑音行列と観測雑音行列を読込み、前記初期値と前記初期誤差共分散行列と前記駆動雑音行列と前記観測雑音行列を使用して時間外挿値と伝播誤差共分散行列を計算し、前記測位計算部１０１が出力した緯度経度座標値を読込み、前記緯度経度座標値を地球楕円体面座標値に変換し、前記地球楕円体面座標値を使用して観測更新値と更新誤差共分散を計算し、前記地球楕円体面座標による観測更新値を緯度経度座標値に変換して出力する整形処理部１０２と、前記整形処理部１０２が出力した緯度経度座標値を表示する測位結果表示部１０３とを備える。【選択図】図４
2007263637.tif 000002
【特許請求の範囲】
【請求項１】
測位信号を送信する人工衛星から前記測位信号を受信し、受信した前記測位信号を用いて測位計算を行い、測位対象となる測位対象移動体の所在位置の測位座標値を出力する測位計算部と、
前記測位対象移動体の速さと方向を状態変数とし、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記測位計算部からの測位対象移動体の所在位置の測位座標値を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する整形処理部と、
前記整形処理部が出力した前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示部とを備えたことを特徴とする測位装置。
【請求項２】
前記整形処理部は、
前記測位対象移動体の所在位置の座標値を推定するとともに、単位時間前の前記測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記測定対象移動体の状態変数の新たな推定値を算出することを特徴とする請求項１に記載の測位装置。
【請求項３】
前記整形処理部は、
単位時間前の前記測位対象移動体の状態変数の推定値と所定の状態遷移関数を用いて時間外挿値を算出し、
単位時間前の前記測定対象移動体の状態変数の推定誤差共分散行列と所定の駆動雑音共分散行列を用いて伝播誤差共分散行列を算出し、
前記時間外挿値と、前記伝播誤差共分散行列と、所定の観測雑音共分散行列と、前記測位計算部からの前記測位対象移動体の所在位置の測位座標値とを用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定することを特徴とする請求項１に記載の測位装置。
【請求項４】
前記整形処理部は、
前記測位対象移動体の所在位置の座標値を推定するとともに、前記伝播誤差共分散行列と、前記観測雑音共分散行列とを用いて、前記測定対象移動体の状態変数の新たな推定値を算出することを特徴とする請求項３に記載の測位装置。
【請求項５】
前記測位計算部は、
前記測位対象移動体の所在位置の測位座標値を地球楕円体面座標値以外の座標値で出力し、
前記整形処理部は、
前記測位計算部から出力された前記測位対象移動体の所在位置の測位座標値を地球楕円体面座標値に変換し、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を地球楕円体面座標値で推定して、地球楕円体面座標値での前記測位計算部からの測位対象移動体の所在位置の測位座標値を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の地球楕円体面座標値での推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の地球楕円体面座標値での推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を地球楕円体面座標値で推定し、前記測位対象移動体の所在位置の推定座標値を前記地球楕円体面座標値以外の座標値に変換し、変換後の前記測位対象移動体の所在位置の推定座標値を前記測位結果表示部に出力し、
前記測位結果表示部は、
前記整形処理部が出力した変換後の前記測位対象移動体の所在位置の推定座標値を表示することを特徴とする請求項１に記載の測位装置。
【請求項６】
前記測位計算部は、
前記測位対象移動体の所在位置の測位座標値を緯度経度座標値で出力し、
前記整形処理部は、
前記測位計算部から出力された前記測位対象移動体の所在位置の測位座標値を地球楕円体面座標値に変換し、
前記測位対象移動体の所在位置の推定座標値を緯度経度座標値に変換し、変換後の前記測位対象移動体の所在位置の推定座標値を前記測位結果表示部に出力することを特徴とする請求項５に記載の測位装置。
【請求項７】
測位対象となる測位対象移動体の速さと方向を状態変数とし、人工衛星からの擬似距離信号を読み込み、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記人工衛星からの擬似距離信号を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する測位計算部と、
前記測位計算部が出力した前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示部とを備えたことを特徴とする測位装置。
【請求項８】
前記整形処理部は、
前記測位対象移動体の所在位置の座標値を推定するとともに、単位時間前の前記測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記測定対象移動体の状態変数の新たな推定値を算出することを特徴とする請求項７に記載の測位装置。
【請求項９】
前記測位計算部は、
単位時間前の前記測位対象移動体の状態変数の推定値と所定の状態遷移関数を用いて時間外挿値を算出し、
単位時間前の前記測定対象移動体の状態変数の推定誤差共分散行列と所定の駆動雑音共分散行列を用いて伝播誤差共分散行列を算出し、
前記時間外挿値と、前記伝播誤差共分散行列と、所定の観測雑音共分散行列と、前記人工衛星からの擬似距離信号とを用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定することを特徴とする請求項７に記載の測位装置。
【請求項１０】
前記測位計算部は、
前記測位対象移動体の所在位置の座標値を推定するとともに、前記伝播誤差共分散行列と、前記観測雑音共分散行列とを用いて、前記測定対象移動体の状態変数の新たな推定値を算出することを特徴とする請求項９に記載の測位装置。
【請求項１１】
前記測位計算部は、
前記人工衛星からの擬似距離信号として、地球中心固定座標値での擬似距離信号を読み込み、
単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を地球中心固定座標値で推定して、地球中心固定座標値での前記人口衛星からの擬似距離信号を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の地球中心固定座標値での推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の地球中心固定座標値での推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を地球中心固定座標値で推定し、前記測位対象移動体の所在位置の推定座標値を地球中心固定座標値以外の座標値に変換し、変換後の前記測位対象移動体の所在位置の推定座標値を前記測位結果表示部に出力し、
前記測位結果表示部は、
前記測位計算部が出力した変換後の前記測位対象移動体の所在位置の推定座標値を表示することを特徴とする請求項７に記載の測位装置。
【請求項１２】
前記測位計算部は、
前記測位対象移動体の所在位置の推定座標値を地球中心固定座標値から緯度経度座標値に変換し、変換後の前記測位対象移動体の所在位置の推定座標値を前記測位結果表示部に出力することを特徴とする請求項１１に記載の測位装置。
【請求項１３】
測位信号を送信する人工衛星から前記測位信号を受信し、受信した前記測位信号を用いて測位計算を行い、測位対象となる測位対象移動体の所在位置の測位座標値を出力する測位計算ステップと、
前記測位対象移動体の速さと方向を状態変数とし、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記測位計算ステップからの測位対象移動体の所在位置の測位座標値を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する整形処理ステップと、
前記整形処理ステップから出力された前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示ステップとを備えたことを特徴とする測位方法。
【請求項１４】
測位対象となる測位対象移動体の速さと方向を状態変数とし、人工衛星からの擬似距離信号を読み込み、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記人工衛星からの擬似距離信号を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する測位計算ステップと、
前記測位計算ステップから出力された前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示ステップとを備えたことを特徴とする測位方法。
【請求項１５】
測位信号を送信する人工衛星から前記測位信号を受信し、受信した前記測位信号を用いて測位計算を行い、測位対象となる測位対象移動体の所在位置の測位座標値を出力する測位計算処理と、
前記測位対象移動体の速さと方向を状態変数とし、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記測位計算処理からの測位対象移動体の所在位置の測位座標値を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する整形処理処理と、
前記整形処理処理から出力された前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示処理とをコンピュータに実行させることを特徴とするプログラム。
【請求項１６】
測位対象となる測位対象移動体の速さと方向を状態変数とし、人工衛星からの擬似距離信号を読み込み、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記人工衛星からの擬似距離信号を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する測位計算処理と、
前記測位計算処理から出力された前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示処理とをコンピュータに実行させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、測位衛星から送信される測位信号を用いて車両等の移動体の位置を決定する測位技術に関する。
【背景技術】
【０００２】
従来の測位装置では、車両が都市部の高層建物街を走行中に衛星電波の遮蔽や多重伝搬等が生じ、著しい走行軌跡の乱れ、測位精度の劣化が生じていた。これは、単なるパラメータ調整によって解決するのは困難であり、測位計算の処理手順を見直す必要があった（例えば特許文献１）。
【特許文献１】特開平６−３４７５３１号公報
【発明の開示】
【発明が解決しようとする課題】
【０００３】
上記のように、厳しい受信環境下では、走行軌跡の乱れが発生する。走行軌跡の乱れの問題点として、例えば、以下のものがある。
厳しい受信環境下では、車両等の移動体が徐行や停車中でも、マルチパス（多重伝搬）等によって、現実には有り得ない測位結果を示すことがある。
また、一時的に測位精度が大幅に劣化し、走行軌跡が車線を外れたり、道路の外に飛び出したりして、カーナビゲーション装置等の利用者に大きな混乱を招くことがある。
走行履歴を表示した場合、移動体の運動特性が考慮されていないため、走行軌跡が不自然なものとなることがある。
【０００４】
一方、車両にジャイロや車速センサを搭載し、これらの機器からの出力を用いれば（ＩＮＳ（（Ｉｎｅｒｔｉａｌ  Ｎａｖｉｇａｔｉｏｎ  Ｓｙｓｔｅｍ））複合やＤｅａｄ  Ｒｅｃｋｏｎｉｎｇ）、厳しい受信環境下にあっても、走行軌跡の乱れは低減される。
しかしながら、新たなセンサとそれに伴う周辺装置等が必要となるので端末価格が高価となる。
【０００５】
また、従来、滑らかな走行軌跡を得る方法としては、過去のデータを用いて平均処理する移動平均による手法がある。
移動平均による手法では、以下の問題点があった。
現在から過去の一定期間のデータを用いて平均処理する方法は、必然的に時間遅れを伴うものであり、現時点での測位結果を充分に反映したものではない。
また、この方法は、移動体の運動特性を考慮したものではないので、例えば曲進等の移動体特有の走行軌跡を描くことは困難である。
【０００６】
本発明では、車両にジャイロや車速センサを搭載することなく（ＩＮＳ複合することなく）、厳しい受信環境下でも、走行軌跡の乱れを低減し、滑らかな走行軌跡を得ることを主な目的とする。
【課題を解決するための手段】
【０００７】
本発明に係る測位装置は、
測位信号を送信する人工衛星から前記測位信号を受信し、受信した前記測位信号を用いて測位計算を行い、測位対象となる測位対象移動体の所在位置の測位座標値を出力する測位計算部と、
前記測位対象移動体の速さと方向を状態変数とし、単位時間前の前記測位対象移動体の状態変数の推定値を用いて、当該時刻の前記測位対象移動体の状態変数の推定値を算出して、前記当該時刻の測位対象移動体の状態変数の推定値を用いて、当該時刻の測位対象移動体の所在位置の測位座標値を推定して、前記測位計算部からの測位対象移動体の所在位置の測位座標値を用いて、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を算出して、前記当該時刻の測位対象移動体の状態変数の推定値と、前記当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差を用いて、前記状態変数の新たな推定値を算出して前記測位対象移動体の所在位置の座標値を推定し、前記測位対象移動体の所在位置の推定座標値を出力する整形処理部と、
前記整形処理部が出力した前記測位対象移動体の所在位置の推定座標値を表示する測位結果表示部とを備えたことを特徴とする。
【発明の効果】
【０００８】
本発明によれば、衛星電波の遮蔽や多重伝搬等の厳しい受信環境下でも、車両にジャイロや車速センサを搭載することなく、フィルタ計算を行って走行軌跡を再計算することで、走行軌跡の乱れを低減し、滑らかな走行軌跡を得ることができる。
【発明を実施するための最良の形態】
【０００９】
実施の形態１．
図１〜図７を使用して本実施の形態を説明する。
本実施の形態は、測位衛星から送信される測位信号を用いて車両等の移動体の位置を決定する衛星測位システムに関する。
【００１０】
図１は、本実施の形態１に係る衛星測位システム３００の構成を示す。
衛星測位システム３００は、測位信号１を送信する測位衛星２００ａ，２００ｂ，２００ｃ・・・２００ｎと、測位信号１を受信し移動体（車両）の位置を決定する測位装置１００とを備える。また、測位衛星は、測位衛星２００ａ，２００ｂ，２００ｃ，２００ｎの４機のみ図示しているが、例示であり４機にかぎらない。
【００１１】
本実施の形態に係る測位装置１００の特徴の一つは、従来、車両の運動モデルとしては、直進運動モデルを用いられていたのに対して、本実施の形態では、車両の速さと向きに着目したモデルを用いた点である。
図２は、従来装置の走行軌跡の例である。また図３は、図２の走行軌跡に対する速度ベクトルの時間的変化（時間は（１）から（２）の方向へ進行）である。
横軸は速度ベクトルの東方向成分、縦軸は速度ベクトルの北方向成分である。速度ベクトルとは、車両の速度の大きさと方向を示すものであり、図３は時間の経過に伴う速度ベクトルの変化の様子を描いたものである。
図３に示すように、滑らかでない走行軌跡に対する速度ベクトルの時間的変化を見ると、変化の仕方が急激であり実際の車両の走行に合致していないと考えられる。
そこで、本実施の形態では、速度ベクトルに着目して、速度ベクトルの大きさ（車速）と向き（ヘッディング）を状態変数としたモデル（車速＋ヨーレートモデル）を用いる。車速がアクセル／ブレーキ、ヨーレートがハンドル操作にそれぞれ対応するため、実際の車両走行特性の模擬であり、より良好な走行軌跡が見込まれる。
これは、以下の理由による。
【００１２】
【数１】
2007263637.tif 000003
【００１３】
これより、走行軌跡は速度ベクトルを時間の経過に沿って積分したものであるから、車両の速度ベクトルの変化を緩やかにすれば、走行軌跡を滑らかにすることが出来ることが分る。
【００１４】
図１を参照して、衛星測位システム３００の概要を説明する。
（１）測位衛星２００ａ・・・２００ｎ等は測位信号１を送信する。
（２）測位装置１００は、この測位信号１を受信して所定の情報を抽出する。所定の情報は、測位装置１００と測位衛星２００との間の擬似距離、ドプラ周波数、及び搬送波位相などを含む。
（３）測位装置１００は、抽出した情報に基づいて移動体の測位情報を作成する。測位情報は、移動体の位置、速度などを含む。
【００１５】
図４は、測位装置１００の構成を示す。
図４に示すように、測位装置１００は、測位計算部１０１、整形処理部１０２、測位結果表示部１０３とを備える。
また、整形処理部１０２は、設定データ読込み部１０２１、時間伝播処理部１０２２、入力データ読込み部１０２３、単位・座標変換部１０２４、観測更新処理部１０２５、単位・座標逆変換部１０２６とを備える。
整形処理部１０２は、測位計算部１０１の出力結果（緯度経度）を用いて、走行軌跡の再計算を行う。再計算とは、図２に示されている（滑らかでない）走行軌跡の計算のことである。本実施の形態に係る測位装置の目的は、測位計算部１０１で計算した走行軌跡を改善することである。本実施の形態に係る測位装置の整形処理部１０２は、測位計算部１０１で計算した走行軌跡結果を入力データとして、走行軌跡を計算し直して出力している。このようにして、測位計算部１０１の走行軌跡を滑らかにしている。
【００１６】
ここで、本実施の形態に係る測位装置１００の動作の概要を図５を参照して説明する。
測位計算部１０１（前段）（測位計算ステップ）は、測位衛星２００からの測位信号１を受信する。この測位信号は地球中心固定座標である。次に、測位計算部１０１は、直進運動モデルによる測位計算を行って、移動体の所在位置の座標値を算出する。次に、算出した移動体の座標値を地球中心固定座標から緯度経度座標に座標変換を行い、緯度経度座標の座標値を測位結果として整形処理部１０２に出力する。
整形処理部１０２（整形処理ステップ）は、測位計算部１０１の出力結果（緯度経度）を用いて、走行軌跡の再計算を行う。なお、車両高度は、緯度経度に比較して変化量が小さいので考慮しないものとする。緯度経度座標は、基準座標（初期座標）を原点とする東方向と北方向の地球楕円体面座標に変換する（後述の数１２の計算式を利用）。
この座標において、状態変数を、東方向変位［ｍ］、北方向変位［ｍ］、速さ［ｍ／ｓ］、ヘッディング［ｒａｄ］、アクセル［ｍ／ｓ＾２］、ヨーレート［ｒａｄ／ｓ］の６変数とする。また、観測値は、測位計算部１０１の出力結果である経度［ｄｅｇｍｉｎ］、緯度［ｄｅｇｍｉｎ］とする。状態遷移方程式は非線形系となるので、拡張カルマンフィルタを構成して、新たな走行軌跡を求める。つまり、速度、ヘッディング、アクセル、ヨーレートを算出し、東方向と北方向の位置と、速度、ヘッディング、アクセル、ヨーレートを状態変数としたカルマンフィルタのモデルを作成する。このモデルにより模擬的なＩＮＳ複合計算を行う（後述の数７、８、１３、１５、１７の計算式を利用）。得られた結果は、地球楕円体面座標から緯度経度座標へ変換する（後述の数１８の計算式を利用）。なお、具体的な計算式については、後述する。
測位結果表示部１０３（測位結果表示ステップ）は、整形処理部１０２の測位結果を表示す。
なお、上記において、地球楕円体面座標とは、基準点を適宜設定して（ここでは車両の初期値）、この基準点を原点として、地球上の位置を、東方向と北方向で［ｍ］単位で示したものである。緯度経度座標の単位が［ｄｅｇ］であるのに対して、地球楕円体面座標の単位は［ｍ］である。なお、上述したように、車両高度は、緯度経度に比較して変化量が小さいので考慮していない。
【００１７】
従来の車速センサとジャイロを用いた複合方式では、速度とヨーレートを外部入力変数として用いてカルマンフィルタを構成している。
一方、本実施の形態では、外部入力変数は、緯度／経度情報のみであり、速度及びその時間変化率であるアクセル、ヨーレート及びその時間積分値であるヘッディングを状態変数として含むカルマンフィルタのモデルを作成している。
位置ｘ、ｙと速度ｖとヘッディングφとの関係は、以下の通りである。
位置ｘ：ｘ＝ｖ×ｃｏｓφ    ［ｍ］
位置ｙ：ｙ＝ｖ×ｓｉｎφ    ［ｍ］
また、緯度経度の単位系から、位置、速度、ヘッディング、アクセル、ヨーレートの単位系へ変更するために、座標変換を行っている。
【００１８】
図１２は、本実施の形態及び以下に述べる実施の形態に示す。測位装置１００のハードウェア資源の一例を示す図である。
図１２において、測位装置１００は、プログラムを実行するＣＰＵ９１１（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサともいう）を備えている。ＣＰＵ９１１は、バス９１２を介して、例えば、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）９１３、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）９１４、通信ボード９１５、表示装置９０１、磁気ディスク装置９２０と接続され、これらのハードウェアデバイスを制御する。表示装置９０１は、タッチパネル式であってもよい。更に、ＣＰＵ９１１は、キーボード９０２、マウス９０３、ＦＤＤ９０４（Ｆｌｅｘｉｂｌｅ  Ｄｉｓｋ  Ｄｒｉｖｅ）、コンパクトディスク装置９０５（ＣＤＤ）と接続していてもよい。また、磁気ディスク装置９２０の代わりに、フラッシュメモリ、光ディスク装置、メモリカード読み書き装置などの記憶装置でもよい。
ＲＡＭ９１４は、揮発性メモリの一例である。ＲＯＭ９１３、ＦＤＤ９０４、ＣＤＤ９０５、磁気ディスク装置９２０の記憶媒体は、不揮発性メモリの一例である。これらは、記憶装置あるいは記憶部の一例である。
通信ボード９１５、キーボード９０２、ＦＤＤ９０４などは、入力部、入力装置の一例である。
また、通信ボード９１５、表示装置９０１などは、出力部、出力装置の一例である。
【００１９】
通信ボード９１５は、図１に示すように、測位衛星２００からの測位信号１を受信可能なように構成されている。また、通信ボード９１５は、地上無線通信網と接続することも可能であり、地上無線通信網を通じて、インターネット、ＷＡＮ（ワイドエリアネットワーク）などに接続されていても構わない。
磁気ディスク装置９２０には、オペレーティングシステム９２１（ＯＳ）、ウィンドウシステム９２２、プログラム群９２３、ファイル群９２４が記憶されている。プログラム群９２３のプログラムは、ＣＰＵ９１１、オペレーティングシステム９２１、ウィンドウシステム９２２により実行される。
【００２０】
上記プログラム群９２３には、以下に述べる説明において「〜部」、「〜手段」として説明する機能を実行するプログラムが記憶されている。プログラムは、ＣＰＵ９１１により読み出され実行される。
ファイル群９２４には、以下に述べる説明において、「〜の判定結果」、「〜の計算結果」、「〜の処理結果」、「〜の評価結果」等として説明する情報やデータや信号値や変数値やパラメータが、「〜ファイル」や「〜データベース」の各項目として記憶されている。「〜ファイル」や「〜データベース」は、ディスクやメモリなどの記録媒体に記憶される。ディスクやメモリになどの記憶媒体に記憶された情報やデータや信号値や変数値やパラメータは、読み書き回路を介してＣＰＵ９１１によりメインメモリやキャッシュメモリに読み出され、抽出・検索・参照・比較・演算・計算・処理・編集・出力・印刷・表示などのＣＰＵの動作に用いられる。抽出・検索・参照・比較・演算・計算・処理・編集・出力・印刷・表示のＣＰＵの動作の間、情報やデータや信号値や変数値やパラメータは、メインメモリ、レジスタ、キャッシュメモリ、バッファメモリ等に一時的に記憶される。
また、以下で説明するフローチャートの矢印の部分は主としてデータや信号の入出力を示し、データや信号値は、ＲＡＭ９１４のメモリ、ＦＤＤ９０４のフレキシブルディスク、ＣＤＤ９０５のコンパクトディスク、磁気ディスク装置９２０の磁気ディスク、その他光ディスク、ミニディスク、ＤＶＤ等の記録媒体に記録される。また、データや信号は、バス９１２や信号線やケーブルその他の伝送媒体によりオンライン伝送される。
【００２１】
また、以下に述べる説明において「〜部」、「〜手段」として説明するものは、「〜回路」、「〜装置」、「〜機器」、「手段」であってもよく、また、「〜ステップ」、「〜手順」、「〜処理」であってもよい。すなわち、「〜部」、「〜手段」として説明するものは、ＲＯＭ９１３に記憶されたファームウェアで実現されていても構わない。或いは、ソフトウェアのみ、或いは、素子・デバイス・基板・配線などのハードウェアのみ、或いは、ソフトウェアとハードウェアとの組み合わせ、さらには、ファームウェアとの組み合わせで実施されても構わない。ファームウェアとソフトウェアは、プログラムとして、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ等の記録媒体に記憶される。プログラムはＣＰＵ９１１により読み出され、ＣＰＵ９１１により実行される。すなわち、プログラムは、以下に述べる「〜部」、「〜手段」としてコンピュータを機能させるものである。あるいは、以下に述べる「〜部」、「〜手段」の手順や方法をコンピュータに実行させるものである。
【００２２】
このように、本実施の形態及び以下に述べる実施の形態に示す測位装置１００は、処理装置たるＣＰＵ、記憶装置たるメモリ、磁気ディスク等、入力装置たるタッチパネル、キーボード、マウス等、出力装置たる表示装置、通信ボード等を備えるコンピュータであり、上記したように「〜部」、「〜手段」として示された機能をこれら処理装置、記憶装置、入力装置、出力装置を用いて実現するものである。
【００２３】
図６は、整形処理部１０２と測位結果表示部１０３の動作例を示すフローチャートである。図６を参照して、整形処理部１０２と測位結果表示部１０３の動作の概略を説明する。
【００２４】
Ｓ１１において、設定データ読込み部１０２１は初期値（ｘｓ（０））や初期誤差共分散行列（Ｓ（０））、駆動雑音共分散行列（Ｑ）、観測雑音共分散行列（Ｒ）を読込む。
初期値（ｘｓ（０））や初期誤差共分散行列（Ｓ（０））、駆動雑音共分散行列（Ｑ）、観測雑音共分散行列（Ｒ）のデータは、測位装置１００のメモリに格納されている。設定データ読込み部１０２１は、このメモリからデータを読み込む。以下に、用語を説明する。
初期値（ｘｓ（０））とは、初期時刻（計算開始時刻）における各状態変数の値（設定値）である。
初期誤差共分散行列（Ｓ（０））とは、各状態変数の初期値の誤差に関する共分散行列（自己分散値と相互分散値を要素とする行列）である。
駆動雑音共分散行列（Ｑ）とは、各状態変数に関して時間刻みτの間に変化する程度を示す統計量である。
観測雑音共分散行列（Ｒ）とは、観測値の誤差に関する共分散行列（自己分散値と相互分散値を要素とする行列）である。
ここで、状態変数は、東方向変位［ｍ］、北方向変位［ｍ］、速さ［ｍ／ｓ］、ヘッディング［ｒａｄ］、アクセル［ｍ／ｓ＾２］、ヨーレート［ｒａｄ／ｓ］の６変数である。
【００２５】
Ｓ１２において、状態変数の時間伝播処理部１０２２は初期値（ｘｓ（０））や初期誤差共分散行列（Ｓ（０））、駆動雑音共分散行列（Ｑ）、観測雑音共分散行列（Ｒ）を用いて、時間外挿値（ｘｐ（ｔ））と伝播誤差共分散行列（Ｐ（ｔ））を計算する。
時間外挿値（ｘｐ（ｔ））とは、時間τ経過後の各状態変数の推定値である。
伝播誤差共分散行列（Ｐ（ｔ））とは、時間τ経過後の各状態変数の推定誤差の共分散行列である。
【００２６】
Ｓ１３において、入力データ読込み部１０２３は測位計算部１０１の出力結果である緯度、経度を読込む。
【００２７】
Ｓ１４において、単位・座標変換部１０２４は入力データ読込み部１０２３が読込んだ緯度経度座標を、地球楕円体面座標に変換する。
【００２８】
Ｓ１５において、状態変数の観測更新処理部１０２５は単位・座標変換部１０２４が計算した地球楕円体面座標を用いて、観測更新値（ｘｓ（ｔ）と更新誤差共分散（Ｓ（ｔ））を計算する。
観測更新値（ｘｓ（ｔ））とは、観測値（緯度経度情報）を得ることによって推定した各状態変数の値である。
更新誤差共分散（Ｓ（ｔ））とは、観測値（緯度経度情報）を得ることによる各状態変数の推定誤差の共分散行列である。
【００２９】
Ｓ１６において、単位・座標逆変換部１０２６は観測更新処理部１０２３が計算した地球楕円体面座標を、緯度経度座標へ変換する。
【００３０】
Ｓ１７において、測位結果表示部１０３は単位・座標変換部１０２４が計算した緯度経度座標（測位情報）を表示す。
【００３１】
このように、整形処理部１０２は、状態変数に測位移動体の速さと方向を含み、初期値と初期誤差共分散行列と駆動雑音行列と観測雑音行列を読込み、前記初期値と前記初期誤差共分散行列と前記駆動雑音行列と前記観測雑音行列を使用して前記状態変数の時間外挿値と伝播誤差共分散行列を計算し、前記測位計算部が出力した緯度経度座標値を読込み、前記緯度経度座標値を地球楕円体面座標値に変換し、前記地球楕円体面座標値を使用して前記状態変数の観測更新値と更新誤差共分散を計算し、前記地球楕円体面座標による観測更新値を緯度経度座標値に変換して出力する。
【００３２】
（拡張カルマンフィルタによる実時間推定）
次に、整形処理部１０２による処理内容の詳細について以下に説明する。
地球楕円体面座標系における東方向変位［ｍ］、北方向変位［ｍ］、速さ［ｍ／ｓ］、ヘッディング［ｒａｄ］、アクセル［ｍ／ｓ＾２］、ヨーレート［ｒａｄ／ｓ］の各変数は、以下の拡張カルマンフィルタを用いて実時間推定する。
【００３３】
（１）対象システムの状態方程式の記述
【００３４】
【数２】
2007263637.tif 000004
【００３５】
ここで、ｔは時刻、ｘ（ｔ）は状態変数、ｆ（ｘ（ｔ））は状態遷移関数、ξ（ｔ）はシステム雑音である。また、状態変数と状態遷移関数の内容は以下の通りである。
【００３６】
【数３】
2007263637.tif 000005
【００３７】
（２）対象システムの観測方程式の記述
【００３８】
【数４】
2007263637.tif 000006
【００３９】
ここで、ｍ（ｔ）は観測量、ｇ（ｘ（ｔ））は観測関数、ζ（ｔ）は観測雑音である。また、観測量の内容は本実施の形態の処理を付加する前の測位結果であり、以下の通りである。
【００４０】
【数５】
2007263637.tif 000007
【００４１】
上記の非線形対象系と線形観測系に対する拡張カルマンフィルタと単位・座標変換等の処理を、下記のように構成する。
【００４２】
（３）初期設定データ読込み
【００４３】
下記の初期値、設定値を読み込む。
【００４４】
【数６】
2007263637.tif 000008
【００４５】
（４）時間外挿計算
【００４６】
時間外挿値ｘｐ（ｔ）は以下のように求める。
【００４７】
【数７】
2007263637.tif 000009
【００４８】
ここで、ｘｓ（ｔ）は観測更新値、τは外挿時間である。
初めて時間外挿値ｘｐ（ｔ）を求める場合は、読み込んだ初期平滑値ｘｓ（０）を用いて、ｘｐ（ｔ）＝ｘｓ（０）＋ｆ（ｘｓ（０））・τとして時間外挿値ｘｐ（ｔ）を求める。
既に、現在時刻からτ［ｓ］前に観測更新値ｘｓ（ｔ−τ）を算出している場合は、このｘｓ（ｔ−τ）を数７に適用して時間外挿値ｘｐ（ｔ）を求める。
【００４９】
（５）誤差共分散伝播計算
【００５０】
伝播誤差共分散Ｐ（ｔ）は以下のように求める。
【００５１】
【数８】
2007263637.tif 000010
【００５２】
ここで、Ｓ（ｔ）は更新誤差共分散である。また、Ａ（ｔ）及びＦ（ｔ）の内容は以下の通りである。Ｑは、数６に示すように、駆動雑音共分散行列である。
【００５３】
【数９】
2007263637.tif 000011
【００５４】
（６）入力データ読込み
【００５５】
下記の入力データを読み込む。
【００５６】
【数１０】
2007263637.tif 000012
【００５７】
（７）単位変換
【００５８】
下記の手順で、緯度経度座標［ｄｅｇ］から地球楕円体面座標［ｍ］への変換を行う。
【００５９】
まず、緯度経度初期値（φ０，λ０）［ｄｅｇ］を、基準座標値（原点）とする。
【００６０】
【数１１】
2007263637.tif 000013
【００６１】
次に、緯度経度座標系から地球楕円体面座標系へ変換を行い、単位を［ｄｅｇ］から［ｍ］に直す。
【００６２】
【数１２】
2007263637.tif 000014
【００６３】
（８）観測更新計算
【００６４】
観測更新値ｘｓ（ｔ）は以下のように求める。
【００６５】
【数１３】
2007263637.tif 000015
【００６６】
ここで、ｍ（ｔ）及びｇ（ｘｐ（ｔ））は以下の通りである。
【００６７】
【数１４】
2007263637.tif 000016
【００６８】
また、Ｋ（ｔ）はカルマンゲインである。カルマンゲインＫ（ｔ）は以下のように求める。
【００６９】
【数１５】
2007263637.tif 000017
【００７０】
なお、Ｃ（ｔ）及びＮ（ｔ）は以下の通りである。
【００７１】
【数１６】
2007263637.tif 000018
【００７２】
（９）誤差共分散更新計算
【００７３】
更新誤差共分散Ｓ（ｔ）は以下のように求める。
【００７４】
【数１７】
2007263637.tif 000019
【００７５】
（１０）単位変換
下記の手順で、地球楕円体面座標［ｍ］から緯度経度座標［ｄｅｇ］への変換を行う。
【００７６】
まず、地球楕円体面座標系から緯度経度座標系へ変換を行い、単位を［ｍ］から［ｄｅｇ］に直す。
【００７７】
【数１８】
2007263637.tif 000020
【００７８】
次に、基準座標値（φ０，λ０）［ｄｅｇ］を加えて、元の緯度経度座標値［ｄｅｇ］を求める。
【００７９】
【数１９】
2007263637.tif 000021
【００８０】
（１１）結果表示
下記の計算結果を書き出す。
λ：経度（ｌｏｎｇｉｔｕｄｅ）［ｄｅｇ］
φ：緯度（ｌａｔｉｔｕｄｅ）［ｄｅｇ］
【００８１】
図２に示した従来装置の測位結果に対して本実施の形態に係る整形処理部１０２の整形処理を適用した場合の走行軌跡を図７に示す。
また図８は、図７の走行軌跡に対する速度ベクトルの時間的変化（時間は（１）から（２）の方向へ進行）である。
従来装置の速度ベクトルの時間的変化は、図３に示したように急激であったが、本実施の形態に係る整形処理部１０２の整形処理を適用した結果、速度ベクトルの時間的変化は緩やかに変化し、それに対応して走行軌跡が滑らかになったことが図７よりわかる。
【００８２】
このように、本実施の形態では、測位計算部１０１が、人工衛星から測位信号を受信し、受信した測位信号を用いて測位計算を行い、測位対象となる測位対象移動体の所在位置の測位座標値を出力する。
そして、整形処理部１０２が、測位対象移動体の速さと方向を状態変数とし（ｘ（ｔ）（数３））、単位時間τ前の測位対象移動体の状態変数の推定値（τ時間前の更新観測値：ｘｓ（ｔ−τ）又は初期平滑値：ｘｓ（０））を用いて、当該時刻の測位対象移動体の状態変数の推定値（時間外挿値：ｘｐ（ｔ）（数７））を算出する。
更に、整形処理部１０２は、当該時刻の測位対象移動体の状態変数の推定値（時間外挿値：ｘｐ（ｔ））を用いて、当該時刻の測位対象移動体の所在位置の測位座標値（伝播誤差共分散：Ｐ（ｔ）（数８））を推定して、測位計算部１０１からの測位対象移動体の所在位置の測位座標値（ｍ（ｔ）（数４）、ｍ１：緯度、ｍ２：経度（数１０））を用いて、当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差（更新誤差共分散：Ｓ（ｔ）（数１７））を算出する。
また、整形処理部１０２は、当該時刻の測位対象移動体の状態変数の推定値（時間外挿値：ｘｐ（ｔ））と、当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差（更新誤差共分散：Ｓ（ｔ））を用いて、状態変数の新たな推定値（観測更新値ｘｓ（ｔ）（数１３））を算出して測位対象移動体の所在位置の座標値を推定し、測位対象移動体の所在位置の推定座標値を測位結果表示部１０３に出力する。
測位結果表示部１０３は、整形処理部１０２が出力した測位対象移動体の所在位置の推定座標値を表示す。
【００８３】
以上のように本実施の形態に係る測位装置は、厳しい受信環境下においても測位精度の劣化を抑えることができ、滑らかな走行軌跡が得られる。
特に、ジャイロや車速センサを搭載することなく、模擬的なＩＮＳ複合計算を行う（データ処理のみでＩＮＳ複合を模擬的に実現する）ことにより、安価で、時間遅れなく、曲進等の移動体特有の走行軌跡を得ることができる。
【００８４】
また、以上のように本実施の形態に係る測位装置は、従来の測位装置の後段に整形処理部１０２を付加するのみで、衛星電波の遮蔽や多重伝搬等による精度の劣化を抑えた測位装置を得ることができる。
【００８５】
以上のように本実施の形態に係る測位方法は、衛星電波の遮蔽や多重伝搬等によって乱れた走行軌跡を、フィルタ計算を行って走行軌跡を再計算した結果、軌跡は滑らかになり、測位精度の劣化を抑えることができる。また、走行軌跡の最適化や微調整は、フィルタ処理における初期値や初期誤差共分散行列、駆動雑音行列、観測雑音行列等のパラメータの調整のみで行うことができる。
【００８６】
また、図４の構成に加え、測位計算部１０１の測位のための補正データを受信する補正データ受信部を追加し、測位計算部１０１は、測位衛星２００からの測位信号１とともに補正データ受信部が受信した補正データを使用して測位計算を行って緯度経度座標値を出力するようにしてもよい。
【００８７】
このように、本実施の形態では、測位信号を送信する人工衛星から前記測位信号を受信して所定の情報を抽出し、前記所定の情報を使用して測位計算を行って緯度経度座標値を出力する測位計算部と、状態変数に測位移動体の速さと方向を含み、初期値と初期誤差共分散行列と駆動雑音行列と観測雑音行列を読込み、前記初期値と前記初期誤差共分散行列と前記駆動雑音行列と前記観測雑音行列を使用して前記状態変数の時間外挿値と伝播誤差共分散行列を計算し、前記測位計算部が出力した緯度経度座標値を読込み、前記緯度経度座標値を地球楕円体面座標値に変換し、前記地球楕円体面座標値を使用して前記状態変数の観測更新値と更新誤差共分散を計算し、前記地球楕円体面座標による観測更新値を緯度経度座標値に変換して出力する整形処理部と、前記整形処理部が出力した緯度経度座標値を表示する測位結果表示部とを備えた測位装置について説明した。
【００８８】
また、本実施の形態では、前記測位計算部の測位のための補正データを受信する補正データ受信部を備え、前記測位計算部が、測位信号を送信する人工衛星から前記測位信号を受信して所定の情報を抽出し、前記所定の情報を使用するとともに、前記補正データ受信部で受信した補正データを使用して測位計算を行って緯度経度座標値を出力する測位装置について説明した。
【００８９】
実施の形態２．
図９は、本実施の形態に係る測位装置１０４の構成を示す。
図９に示すように、測位装置１０４は、測位計算部１０５、測位結果表示部１０３とを備える。
測位計算部１０５は、設定データ読込み部１０２１、時間伝播処理部１０２２、入力データ読込み部１０２７、観測更新処理部１０２５、座標変換部１０２８とを備える。
【００９０】
ここで、本実施の形態に係る測位装置１０４の動作の概要を図１０を参照して説明する。
測位計算部１０５（測位計算ステップ）は、走行軌跡の計算を行う。
状態変数は、地球中心固定座標において、ｘ方向［ｍ］、ｙ方向［ｍ］、ｚ方向［ｍ］、受信機クロック［ｍ］、速さ［ｍ／ｓ］、エレベーション［ｒａｄ］、ヘッディング［ｒａｄ］、クロックドリフト［ｍ／ｓ］、アクセル［ｍ／ｓ＾２］、ピッチレート［ｒａｄ／ｓ］、ヨーレート［ｒａｄ／ｓ］、クロックドリフトレート［ｍ／ｓ＾２］の９変数とする。
また、観測値は、測位衛星２００からの測位信号１に含まれる測位衛星２００との間の擬似距離［ｍ］である。測位計算部１０５が受信する測位信号は地球中心固定座標である。状態方程式と観測方程式は非線形系となるので、拡張カルマンフィルタを構成する。
得られた結果は、地球中心固定座標から緯度経度高度座標へ変換する。
測位結果表示部１０３（測位結果表示ステップ）は、測位計算部１０５の測位結果を表示する。
【００９１】
図１１は、測位計算部１０５と測位結果表示部１０３の動作を示すフローチャートである。図１１を参照して、測位計算部１０５と測位結果表示部１０３の動作の概略を説明する。
【００９２】
Ｓ１１において、設定データ読込み部１０２１は初期値（ｘｓ（０））や初期誤差共分散行列（Ｓ（０））、駆動雑音共分散行列（Ｑ）、観測雑音共分散行列（Ｒ）を読込む。
【００９３】
Ｓ１２において、時間伝播処理部１０２２は初期値（ｘｓ（０））や初期誤差共分散行列（Ｓ（０））、駆動雑音共分散行列（Ｑ）、観測雑音共分散行列（Ｒ）を用いて、時間外挿値（ｘｐ（ｔ））と伝播誤差共分散行列（Ｐ（ｔ））を計算する。
【００９４】
Ｓ１８において、入力データ読込み部１０２７は人工衛星からの擬似距離信号を読込む。
【００９５】
Ｓ１５において、観測更新処理部１０２５は入力データ読込み部１０２７が読み込んだ擬似距離信号を用いて、観測更新値（ｘｓ（ｔ））と更新誤差共分散（Ｓ（ｔ））を計算する。
【００９６】
Ｓ１９において、座標変換部１０２８は観測更新処理部１０２３が計算した地球中心固定座標を、緯度経度高度座標へ変換する。
【００９７】
Ｓ１７において、測位結果表示部１０３は座標変換部１０２８が計算した緯度経度高度座標（測位情報）を表示する。
【００９８】
以上の動作により、本実施の形態においても、図７及び図８と同様に、速度ベクトルの時間的変化が緩やかに変化する、滑らかな走行軌跡を表示することができる。
【００９９】
（拡張カルマンフィルタによる実時間推定）
次に測位計算部１０５による処理内容の詳細について以下に説明する。
地球中心固定座標におけるｘ方向［ｍ］、ｙ方向［ｍ］、ｚ方向［ｍ］、受信機クロック［ｍ］、速さ［ｍ／ｓ］、エレベーション［ｒａｄ］、ヘッディング［ｒａｄ］、クロックドリフト［ｍ／ｓ］、アクセル［ｍ／ｓ＾２］、ピッチレート［ｒａｄ／ｓ］、ヨーレート［ｒａｄ／ｓ］、クロックドリフトレート［ｍ／ｓ＾２］の各変数は、以下の拡張カルマンフィルタを用いて実時間推定する。
【０１００】
（１）対象システムの状態方程式の記述
【０１０１】
【数２０】
2007263637.tif 000022
【０１０２】
ここで、ｔは時刻、ｘ（ｔ）は状態変数、ｆ（ｘ（ｔ））は状態遷移関数、ξ（ｔ）は駆動雑音である。また、状態変数と状態遷移関数の内容は以下の通りである。
【０１０３】
【数２１】
2007263637.tif 000023
【０１０４】
（２）対象システムの観測方程式の記述
【０１０５】
【数２２】
2007263637.tif 000024
【０１０６】
ここで、ｍ（ｔ）は観測量、ｇ（ｘ（ｔ））は観測関数、ζ（ｔ）は観測雑音である。観測量の内容は人工衛星からの擬似距離であり、以下の通りである。
【０１０７】
【数２３】
2007263637.tif 000025
【０１０８】
上記の非線形対象系と線形観測系に対する拡張カルマンフィルタと座標変換の処理を、下記のように構成する。
【０１０９】
（３）初期設定データ読込み
【０１１０】
下記の初期値、設定値を読み込む。
【０１１１】
【数２４】
2007263637.tif 000026
【０１１２】
（４）時間外挿計算
【０１１３】
時間外挿値ｘｐ（ｔ）は以下のように求める。
【０１１４】
【数２５】
2007263637.tif 000027
【０１１５】
ここで、ｘｓ（ｔ）は観測更新値、τは外挿時間である。
【０１１６】
（５）誤差共分散伝播計算
【０１１７】
伝播誤差共分散Ｐ（ｔ）は以下のように求める。
【０１１８】
【数２６】
2007263637.tif 000028
【０１１９】
ここで、Ｓ（ｔ）は更新誤差共分散である。また、Ａ（ｔ）及びＦ（ｔ）の内容は以下の通りである。Ｑは、数２４に示すように、駆動雑音共分散行列である。
【０１２０】
【数２７】
2007263637.tif 000029
【０１２１】
（６）入力データ読込み
【０１２２】
下記の入力データを読み込む。
【０１２３】
【数２８】
2007263637.tif 000030
【０１２４】
（８）観測更新計算
【０１２５】
観測更新値ｘｓ（ｔ）は以下のように求める。
【０１２６】
【数２９】
2007263637.tif 000031
【０１２７】
ここで、ｍ（ｔ）及びｇ（ｘｐ（ｔ））は以下の通りである。
【０１２８】
【数３０】
2007263637.tif 000032
【０１２９】
また、Ｋ（ｔ）はカルマンゲインである。カルマンゲインＫ（ｔ）は以下のように求める。
【０１３０】
【数３１】
2007263637.tif 000033
【０１３１】
なお、Ｃ（ｔ）及びＮ（ｔ）は以下の通りである。
【０１３２】
【数３２】
2007263637.tif 000034
【０１３３】
（９）誤差共分散更新計算
【０１３４】
更新誤差共分散Ｓ（ｔ）は以下のように求める。
【０１３５】
【数３３】
2007263637.tif 000035
【０１３６】
（１０）座標変換
地球中心固定座標から緯度経度高度座標への変換を行う。
【０１３７】
（１１）結果表示
下記の計算結果を書き出す。
λ：経度（ｌｏｎｇｉｔｕｄｅ）［ｄｅｇ］
φ：緯度（ｌａｔｉｔｕｄｅ）［ｄｅｇ］
ｈ：高度［ｍ］
【０１３８】
このように、本実施の形態では、測位計算部１０５が、測位対象となる測位対象移動体の速さと方向を状態変数とし（ｘ（ｔ）（数２１））、人工衛星からの擬似距離信号を読み込み（ｍ（ｔ）（数２２）、ｍ１：人工衛星からの擬似距離（数２８））、単位時間τ前の測位対象移動体の状態変数の推定値（τ時間前の更新観測値：ｘｓ（ｔ−τ）又は初期平滑値：ｘｓ（０））を用いて、当該時刻の前記測位対象移動体の状態変数の推定値（時間外挿値：ｘｐ（ｔ）（数２５））を算出する。
更に、測位計算部１０５は、当該時刻の測位対象移動体の状態変数の推定値（時間外挿値：ｘｐ（ｔ））を用いて、当該時刻の測位対象移動体の所在位置の測位座標値（伝播誤差共分散：Ｐ（ｔ）（数２６））を推定して、人工衛星からの擬似距離信号を用いて、当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差（更新誤差共分散：Ｓ（ｔ）（数３３））を算出する。
また、測位計算部１０５は、当該時刻の測位対象移動体の状態変数の推定値（時間外挿値：ｘｐ（ｔ））と、当該時刻の測位対象移動体の所在位置の測位座標値の推定誤差（更新誤差共分散：Ｓ（ｔ））を用いて、状態変数の新たな推定値（観測更新値ｘｓ（ｔ）（数２９））を算出して測位対象移動体の所在位置の座標値を推定し、測位対象移動体の所在位置の推定座標値を測位結果表示部１０３に出力する。
測位結果表示部１０３は、測位計算部１０５が出力した測位対象移動体の所在位置の推定座標値を表示する。
【０１３９】
本実施の形態に係る測位装置によれば、厳しい受信環境下においても測位精度の劣化を抑えることができ、滑らかな走行軌跡が得られる。
特に、ジャイロや車速センサを搭載することなく、模擬的なＩＮＳ複合計算を行う（データ処理のみでＩＮＳ複合を模擬的に実現する）ことにより、安価で、時間遅れなく、曲進等の移動体特有の走行軌跡を得ることができる。
また、実施の形態１では、前段の測位計算部１０１が計算した走行軌跡（緯度経度座標値）を用いて、整形処理部１０２が、再度測位計算を行って滑らかな走行軌跡を求めるものであった。一方、実施の形態２では、人工衛星が送信する測位信号を直接用いて、１回の測位計算で滑らかな走行軌跡を求める。このため、実施の形態１に比べて装置構成をシンプルにすることができ、また、走行軌跡計算に要する時間を短縮することができる。
【０１４０】
また、図９の構成に加え、測位計算部１０５の測位のための補正データを受信する補正データ受信部を追加し、測位計算部１０５は、測位衛星２００からの測位信号１とともに補正データ受信部が受信した補正データを使用して測位計算を行って緯度経度座標値を出力するようにしてもよい。
【０１４１】
このように、本実施の形態では、測位信号を送信する人工衛星から前記測位信号を受信して所定の情報を抽出し、前記所定の情報を地球中心固定座標値として読込み、状態変数に測位移動体の速さと方向を含み、初期値と初期誤差共分散行列と駆動雑音行列と観測雑音行列を読込み、前記初期値と前記初期誤差共分散行列と前記駆動雑音行列と前記観測雑音行列を使用して前記状態変数の時間外挿値と伝播誤差共分散行列を計算し、前記所定の情報を使用して前記状態変数の観測更新値と更新誤差共分散を計算し、前記地球中心固定座標による観測更新値を緯度経度高度座標値に変換して出力する測位計算部と、前記測位計算部が出力した緯度経度座標値を表示する測位結果表示部とを備えた測位装置について説明した。
【０１４２】
また、本実施の形態では、前記測位計算部が測位のための補正データを受信する補正データ受信部を備え、前記測位計算部が、測位信号を送信する人工衛星から前記測位信号を受信して所定の情報を抽出し、前記所定の情報を使用するとともに、前記補正データ受信部で受信した補正データを使用して測位計算を行って緯度経度座標値を出力する測位装置について説明した。
【図面の簡単な説明】
【０１４３】
【図１】実施の形態１及び２における衛星測位システムの構成例を示す図。
【図２】従来装置の走行軌跡の例を示す図。
【図３】従来装置の走行軌跡に対する速度ベクトルの時間的変化を示す図。
【図４】実施の形態１における測位装置の構成例を示す図。
【図５】実施の形態１における測位装置の動作の概略を説明する図。
【図６】実施の形態１における整形処理部と測位結果表示部の動作を示すフローチャートである。
【図７】実施の形態１における測位装置を適用した場合の走行軌跡を示す図。
【図８】実施の形態１における測位装置を適用した場合の走行軌跡に対する速度ベクトルの時間的変化を示す図。
【図９】実施の形態２における測位装置の構成例を示す図。
【図１０】実施の形態２における測位装置の動作の概略を説明する図。
【図１１】実施の形態２における整形処理部と測位結果表示部の動作を示すフローチャートである。
【図１２】実施の形態１及び２における測位装置のハードウェア構成例を示す図。
【符号の説明】
【０１４４】
１  測位信号、１００  測位装置、１０１  測位計算部、１０２  整形処理部、１０３  測位結果表示部、１０４  測位装置、１０５  測位計算部、２００  測位衛星、３００  衛星測位システム。
【図１】
2007263637.tif 000036
【図２】
2007263637.tif 000037
【図３】
2007263637.tif 000038
【図４】
2007263637.tif 000039
【図５】
2007263637.tif 000040
【図６】
2007263637.tif 000041
【図７】
2007263637.tif 000042
【図８】
2007263637.tif 000043
【図９】
2007263637.tif 000044
【図１０】
2007263637.tif 000045
【図１１】
2007263637.tif 000046
【図１２】
2007263637.tif 000047