(57)【特許請求の範囲】
【請求項１】
測位演算を行なう衛星信号受信装置であって、
自装置と測位衛星の間の擬似距離を観測する擬似距離観測部と、
前記擬似距離観測部が観測した擬似距離に基づいて、測位時刻より過去の時刻において、当該測位時刻における自装置の予測位置を予測した後に、当該測位時刻において、当該測位時刻における自装置の測位位置を更新する位置計算部と、
前記位置計算部が当該測位時刻より過去の時刻において予測した、当該測位時刻における自装置の予測位置と、前記位置計算部が当該測位時刻において更新した、当該測位時刻における自装置の測位位置と、の差分ベクトルを計算する差分ベクトル計算部と、
当該測位時刻における自装置の速度ベクトルを計算する速度計算部と、
前記差分ベクトル計算部が計算した差分ベクトルが長いほど、さらに当該測位時刻における自装置の速度ベクトル誤差量が大きいほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算する誤差楕円計算部と、
を備えることを特徴とする衛星信号受信装置。
【請求項２】
自装置と測位衛星の間のドップラー周波数を観測するドップラー周波数観測部と、
前記差分ベクトル計算部が計算した差分ベクトルと、前記速度計算部が計算した速度ベクトルと、のなすベクトル間角度を計算し、計算したベクトル間角度に基づいて、速度ベクトル誤差量を計算する速度ベクトル誤差量計算部と、
をさらに備え、
前記速度計算部は、前記ドップラー周波数観測部が観測したドップラー周波数に基づいて、当該測位時刻において、当該測位時刻における自装置の速度ベクトルを計算し、
前記誤差楕円計算部は、前記速度ベクトル誤差量計算部が計算したベクトル間角度が９０°又は２７０°に近いほど、当該測位時刻において、当該測位時刻における前記速度ベクトル誤差量を大きく計算することを特徴とする、請求項１に記載の衛星信号受信装置。
【請求項３】
前記誤差楕円計算部は、ＥＮＵ（Ｌｏｃａｌ  Ｅａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系の高さ（Ｕｐ軸）方向における自装置の測位位置を固定する２次元測位演算が行なわれるときに、ＥＮＵ座標系の水平面（Ｅａｓｔ−Ｎｏｒｔｈ平面）での誤差楕円の長半径及び短半径を、前記差分ベクトル計算部が計算した差分ベクトルに基づいて推定したＥＮＵ座標系の高さ方向における誤差楕円の誤差分散で補正することを特徴とする、請求項１又は請求項２に記載の衛星信号受信装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対して反映させる技術に関する。
【背景技術】
【０００２】
ＧＰＳ（米国）、Ｇａｌｉｌｅｏ（欧州）、ＧＬＯＮＡＳＳ（ロシア）及びＢｅｉＤｏｕ（中国）等の衛星航法システムによる衛星信号を用いて測位演算を行なう、衛星信号受信装置（ＧＮＳＳ（Ｇｌｏｂａｌ  Ｎａｖｉｇａｔｉｏｎ  Ｓａｔｅｌｌｉｔｅ  Ｓｙｓｔｅｍ）受信装置）が、従来から広くユーザに利用されている。
【０００３】
衛星信号受信装置が計算する測位位置は、衛星位置誤差、衛星時計誤差及びマルチパス等の影響による測位誤差を含んでいる。衛星信号受信装置のユーザは、測位位置の確度情報を必要とするため、衛星信号受信装置の測位演算部は、測位位置の精度指標を示す誤差楕円を計算する。
【０００４】
測位位置及び誤差楕円の表示方法を図１に示す。ＥＮＵ（Ｌｏｃａｌ  Ｅａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系は、測位位置を座標の原点とする座標系である。誤差楕円は、測位位置に中心を有し、真位置が楕円の内部に存在する確率が所定の確率（例えば、９５％）であるような軸及び径を有する。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開平０５−３３３１３１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
特許文献１では、誤差楕円の計算方法が開示されている。この文献では、複数の測位衛星の空間的なばらつき指標を示すＤＯＰ（Ｄｉｌｕｔｉｏｎ  Ｏｆ  Ｐｒｅｃｉｓｉｏｎ）に基づいて、誤差楕円を計算している。つまり、測位演算に利用可能な測位衛星数が少ないときや、測位演算に利用可能な測位衛星の配置が偏っているときに、誤差楕円の径が大きく計算される。しかし、図２を用いて以下に説明するように、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対して反映させることができなかった。
【０００７】
マルチパスが測位位置に及ぼす影響を図２に示す。衛星信号受信装置Ｒｘは、高架構造Ｅｌの下部に位置する。測位衛星Ｓ１、Ｓ３、Ｓ４は、衛星信号受信装置Ｒｘから見て天頂以外の方向に位置する。測位衛星Ｓ２は、衛星信号受信装置Ｒｘから見て天頂方向に位置する。反射物体Ｒｆは、衛星信号を反射する物体である。
【０００８】
すると、測位衛星Ｓ１、Ｓ３、Ｓ４からの衛星信号は、高架構造Ｅｌに遮られず、衛星信号受信装置Ｒｘにより受信される。よって、衛星信号受信装置Ｒｘと測位衛星Ｓ１、Ｓ３、Ｓ４の間の擬似距離は、誤差をあまり含んでいない。
【０００９】
しかし、測位衛星Ｓ２からの衛星信号は、高架構造Ｅｌに遮られて、反射物体Ｒｆで反射されるマルチパスを経てから、衛星信号受信装置Ｒｘにより受信される。よって、衛星信号受信装置Ｒｘと測位衛星Ｓ２の間の擬似距離は、誤差を大きく含んでおり、擬似距離に誤差が重畳すれば、擬似距離を使用して計算された測位位置の精度も劣化する。
【００１０】
このように、マルチパスの影響があるときには、マルチパスの影響がないときと比べて、測位位置の精度が劣化するため、誤差楕円の径を大きく計算する必要がある。しかし、マルチパスの影響があるときでも、マルチパスの影響がないときと比べて、複数の測位衛星の空間的なばらつき指標を示すＤＯＰが大きくなるわけではない。よって、特許文献１では、マルチパスの影響があるときでも、マルチパスの影響がないときと比べて、ＤＯＰに基づいた誤差楕円の径は大きく計算されない。つまり、この文献での誤差楕円は、本来あるべき誤差楕円から乖離していた。
【００１１】
そこで、前記課題を解決するために、本発明は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることを目的とする。
【課題を解決するための手段】
【００１２】
上記目的を達成するために、測位時刻より過去の時刻において予測した、当該測位時刻における衛星信号受信装置の予測位置と、当該測位時刻において更新した、当該測位時刻における衛星信号受信装置の測位位置と、の差分ベクトルが長いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算することとした。
【００１３】
具体的には、本発明は、測位演算を行なう衛星信号受信装置であって、自装置と測位衛星の間の擬似距離を観測する擬似距離観測部と、前記擬似距離観測部が観測した擬似距離に基づいて、測位時刻より過去の時刻において、当該測位時刻における自装置の予測位置を予測した後に、当該測位時刻において、当該測位時刻における自装置の測位位置を更新する位置計算部と、前記位置計算部が当該測位時刻より過去の時刻において予測した、当該測位時刻における自装置の予測位置と、前記位置計算部が当該測位時刻において更新した、当該測位時刻における自装置の測位位置と、の差分ベクトルを計算する差分ベクトル計算部と、前記差分ベクトル計算部が計算した差分ベクトルが長いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算する誤差楕円計算部と、を備えることを特徴とする衛星信号受信装置である。
【００１４】
この構成によれば、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルが長くなり、差分ベクトルが長いほど、差分ベクトルを用いて算出したプロセス雑音が大きくなり、プロセス雑音を用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。
【００１５】
また、本発明は、自装置と測位衛星の間のドップラー周波数を観測するドップラー周波数観測部と、前記ドップラー周波数観測部が観測したドップラー周波数に基づいて、当該測位時刻において、当該測位時刻における自装置の速度ベクトルを計算する速度計算部と、前記差分ベクトル計算部が計算した差分ベクトルと、前記速度計算部が計算した速度ベクトルと、のなすベクトル間角度を計算し、計算したベクトル間角度に基づいて、速度ベクトル誤差量を計算する速度ベクトル誤差量計算部と、をさらに備え、前記誤差楕円計算部は、前記速度ベクトル誤差量計算部が計算したベクトル間角度が９０°又は２７０°に近いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算することを特徴とする衛星信号受信装置である。
【００１６】
この構成によれば、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルと速度ベクトルのなすベクトル間角度が９０°又は２７０°に近くなり、ベクトル間角度が９０°又は２７０°に近いほど、ベクトル間角度を用いて算出した速度ベクトル誤差量が大きくなり、速度ベクトル誤差量を用いて算出したプロセス雑音が大きくなり、プロセス雑音を用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。
【００１７】
また、本発明の前記誤差楕円計算部は、ＥＮＵ（Ｌｏｃａｌ  Ｅａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系の高さ（Ｕｐ軸）方向における自装置の測位位置を固定する２次元測位演算が行なわれるときに、ＥＮＵ座標系の水平面（Ｅａｓｔ−Ｎｏｒｔｈ平面）での誤差楕円の長半径及び短半径を、ＥＮＵ座標系の高さ方向における誤差楕円の誤差分散で補正することを特徴とする衛星信号受信装置である。
【００１８】
この構成によれば、測位演算に利用可能な測位衛星数が少なく、高さ方向の測位位置を固定する２次元測位演算を行なうときに、誤差楕円を大きくすることができる。逆に、測位演算に利用可能な測位衛星数が多く、高さ方向の測位位置を固定しない３次元測位演算を行なうときは、誤差楕円を大きくし過ぎないようにすることができる。
【発明の効果】
【００１９】
このように、本発明は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。
【図面の簡単な説明】
【００２０】
【図１】測位位置及び誤差楕円の表示方法を示す図である。
【図２】マルチパスが測位位置に及ぼす影響を示す図である。
【図３】本発明の衛星信号受信装置の構成を示す図である。
【図４】本発明の差分ベクトル及びベクトル間角度を示す図である。
【図５】本発明の差分ベクトル及びベクトル間角度を示す図である。
【図６】本発明の差分ベクトルの計算方法を示す図である。
【図７】本発明の速度ベクトル誤差量の計算方法を示す図である。
【図８】本発明の誤差楕円の計算方法を示す図である。
【図９】本発明の誤差楕円の計算方法の全体構成を示す図である。
【図１０】本発明及び比較例の誤差楕円の時系列グラフである。
【発明を実施するための形態】
【００２１】
添付の図面を参照して本発明の実施形態を説明する。以下に説明する実施形態は本発明の実施の例であり、本発明は以下の実施形態に制限されるものではない。なお、本明細書及び図面において符号が同じ構成要素は、相互に同一のものを示すものとする。
【００２２】
本発明の衛星信号受信装置の構成を図３に示す。衛星信号受信装置Ｒｘは、信号受信部１、追尾処理部２、復調処理部３及び測位演算部４から構成される。
【００２３】
信号受信部１は、アンテナを介してＧＮＳＳ信号を受信する。
【００２４】
追尾処理部２は、擬似距離観測部２１、ドップラー周波数観測部２２及び航法データ抽出部２３から構成される。
【００２５】
擬似距離観測部２１は、衛星信号受信装置Ｒｘと測位衛星の間の擬似距離を観測する。ドップラー周波数観測部２２は、衛星信号受信装置Ｒｘと測位衛星の間のドップラー周波数を観測する。航法データ抽出部２３は、ＧＮＳＳ信号から航法データのビット情報を抽出する。
【００２６】
復調処理部３は、航法データのビット情報を復調（又は復号）し、測位演算に必要なエフェメリス及び衛星時計情報を摘出して出力する。
【００２７】
追尾処理部２及び復調処理部３は、最大可視衛星数に応じて各々必要な個数が定まり、１衛星に対して各々１個を用意する必要がある。
【００２８】
測位演算部４は、位置計算部４１、速度計算部４２、差分ベクトル計算部４３、速度ベクトル誤差量計算部４４及び誤差楕円計算部４５から構成される。
【００２９】
位置計算部４１は、擬似距離観測部２１が観測した擬似距離、復調処理部３が出力したエフェメリス及び衛星時計情報、並びに、衛星信号受信装置Ｒｘが有する受信装置時計情報に基づいて、測位時刻より過去の時刻において、当該測位時刻における衛星信号受信装置Ｒｘの予測位置を予測した後に、当該測位時刻において、当該測位時刻における衛星信号受信装置Ｒｘの測位位置を更新する。
【００３０】
速度計算部４２は、ドップラー周波数観測部２２が観測したドップラー周波数、復調処理部３が出力したエフェメリス及び衛星時計情報、並びに、衛星信号受信装置Ｒｘが有する受信装置時計情報に基づいて、当該測位時刻において、当該測位時刻における衛星信号受信装置Ｒｘの速度ベクトルを計算する。
【００３１】
差分ベクトル計算部４３は、位置計算部４１が当該測位時刻より過去の時刻において予測した、当該測位時刻における衛星信号受信装置Ｒｘの予測位置と、位置計算部４１が当該測位時刻において更新した、当該測位時刻における衛星信号受信装置Ｒｘの測位位置と、の差分ベクトルを計算する。
【００３２】
速度ベクトル誤差量計算部４４は、差分ベクトル計算部４３が計算した差分ベクトルと、速度計算部４２が計算した速度ベクトルと、のなすベクトル間角度を計算し、計算したベクトル間角度に基づいて、速度ベクトル誤差量を計算する。
【００３３】
誤差楕円計算部４５は、差分ベクトル計算部４３が計算した差分ベクトルと、速度ベクトル誤差量計算部４４が計算した速度ベクトル誤差量と、を用いてプロセス雑音を計算し、計算したプロセス雑音と、位置計算部４１で使用した情報と、を用いて誤差楕円の径の長さ及び誤差楕円の軸の方向を計算する。
【００３４】
誤差楕円計算部４５は、差分ベクトル計算部４３が計算した差分ベクトルが長いほど、かつ、速度ベクトル誤差量計算部４４が計算したベクトル間角度が９０°又は２７０°に近いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算する。
【００３５】
測位演算部４は、位置計算部４１が計算した測位位置と、誤差楕円計算部４５が計算した誤差楕円の径の長さ（長半径、短半径）及び誤差楕円の軸の方向と、を外部出力する。外部出力された情報は、衛星信号受信装置Ｒｘを使用するユーザが、ナビゲーション等に利用する。
【００３６】
本発明の差分ベクトル及びベクトル間角度を図４、５に示す。ｐ＾ｋ｜ｋ−１は、位置計算部４１が時刻ｋ−１において予測した、時刻ｋにおける衛星信号受信装置Ｒｘの予測位置である。ｐ＾ｋ｜ｋは、位置計算部４１が時刻ｋにおいて更新した、時刻ｋにおける衛星信号受信装置Ｒｘの測位位置である。ｖ＾ｋ｜ｋは、速度計算部４２が時刻ｋにおいて計算した、時刻ｋにおける衛星信号受信装置Ｒｘの速度ベクトルである。Δｔは、時刻ｋ−１と時刻ｋの間の時間であり、ｐ＾ｋ｜ｋ−１は、ｐ＾ｋ−１｜ｋ−１＋ｖ＾ｋ−１｜ｋ−１Δｔと表わせる。
【００３７】
εｋは、ｐ＾ｋ｜ｋ−１とｐｋ｜ｋの差分ベクトルである。θｋは、εｋとｖ＾ｋ｜ｋのなすベクトル間角度である。図４では、εｋは有限の長さを有し、θｋは０°又は１８０°に近く、図５では、εｋは有限の長さを有し、θｋは９０°又は２７０°に近い。差分ベクトルεｋの長さは、予測位置ｐ＾ｋ｜ｋ−１と実計算した位置ｐ＾ｋ｜ｋとの乖離度合いを示しており、マルチパスが及ぼす影響度合いが大きい状況であれば長くなる。ベクトル間角度θｋは、道路の両脇に高層ビルがある場合や高架下走行時のような、マルチパスが及ぼす影響の度合いが非常に強い状況下では、９０°又は２７０°に近くなる。
【００３８】
ここで、一般的に、進行方向には障害物はなく、進行方向と直角（９０°又は２７０°）の方向に障害物があるため、マルチパスが及ぼす影響を受けている状況下では、図４の状態ではなく、図５の状態になる。つまり、図５に示した状況は、図４に示した状況より、マルチパスが及ぼす影響度合いが強い状況である。
【００３９】
測位演算部４を構成する位置計算部４１、速度計算部４２、差分ベクトル計算部４３、速度ベクトル誤差量計算部４４及び誤差楕円計算部４５の処理の流れを図６〜８に示す。
【００４０】
本発明の差分ベクトルの計算方法を図６に示す。
【００４１】
位置計算部４１及び速度計算部４２は、以下に示す数式１〜７を用いて、ｐ＾ｋ｜ｋ及びｖ＾ｋ｜ｋをそれぞれ計算する（ステップＳ１）。位置計算部４１及び速度計算部４２は、状態方程式及び観測方程式をそれぞれ立式し、カルマンフィルタによる手法で、観測量である擬似距離情報及びドップラー周波数情報を用いて、衛星信号受信装置Ｒｘの測位位置ｐ、速度ベクトルｖ、時計誤差及び時計誤差ドリフトを算出する。状態量ｘは、測位位置ｐ、速度ベクトルｖ、時計誤差及び時計誤差ドリフトを構成成分とし、観測量ｚは、擬似距離の観測量及びドップラー周波数の観測量を構成成分とする。
【００４２】
カルマンフィルタに適用する状態方程式は、数式１のように表わされる。
【数１】
0006546730.tif 000002
ｘｋは、時刻ｋにおける状態量である。Ｆｋは、時刻ｋにおけるシステムモデルである。ωｋは、平均が０、共分散行列がＱｋである多変数正規分布に従うプロセス雑音である。
【００４３】
カルマンフィルタに適用する観測方程式は、数式２のように表わされる。
【数２】
0006546730.tif 000003
ｚｋは、時刻ｋにおける観測量である。Ｈｋは、時刻ｋにおける観測モデルである。ｕｋは、平均が０、共分散行列がＲｋである多変数正規分布に従う観測雑音である。
【００４４】
状態方程式と観測方程式を用いてカルマンフィルタの予測過程と更新過程を計算する。
【００４５】
カルマンフィルタの予測過程は、数式３、４のように表わされる。
【数３】
0006546730.tif 000004
【数４】
0006546730.tif 000005
ｘ＾ｋ｜ｋ−１は、時刻ｋ−１において予測された、時刻ｋにおける推定状態量である。ｘ＾ｋ−１｜ｋ−１は、時刻ｋ−１において更新された、時刻ｋ−１における推定状態量である。Ｐｋ｜ｋ−１は、ｃｏｖ（ｘｋ−ｘ＾ｋ｜ｋ−１）から分かるように、時刻ｋ−１において予測された、時刻ｋにおける推定状態量の誤差共分散行列である。Ｐｋ−１｜ｋ−１は、時刻ｋ−１において更新された、時刻ｋ−１における推定状態量の誤差共分散行列である。
【００４６】
カルマンフィルタの更新過程は、数式５〜７のように表わされる。
【数５】
0006546730.tif 000006
【数６】
0006546730.tif 000007
【数７】
0006546730.tif 000008
ｘ＾ｋ｜ｋは、時刻ｋにおいて更新された、時刻ｋにおける推定状態量である。Ｐｋ｜ｋは、ｃｏｖ（ｘｋ−ｘ＾ｋ｜ｋ）から分かるように、時刻ｋにおいて更新された、時刻ｋにおける推定状態量の誤差共分散行列である。Ｋｋは、時刻ｋにおける最適カルマンゲインである。
【００４７】
差分ベクトル計算部４３は、以下に示す数式８又は数式９を用いて、差分ベクトルεｋを計算する（ステップＳ２）。
【数８】
0006546730.tif 000009
【数９】
0006546730.tif 000010
【００４８】
p ＾ｋ｜ｋは、カルマンフィルタの更新過程で計算された推定状態量ｘ＾ｋ｜ｋから取り出した衛星信号受信装置Ｒｘの測位位置である。p ＾ｋ｜ｋ-1は、推定状態量ｘ＾ｋ｜ｋ-1から取り出した衛星信号受信装置Ｒｘの予測位置である。p ＾ｋ-1｜ｋ-1は、推定状態量ｘ＾ｋ-1｜ｋ-1から取り出した衛星信号受信装置Ｒｘの測位位置である。ｖ ＾ｋ-1｜ｋ-1は、推定状態量ｘ＾ｋ-1｜ｋ-1から取り出した衛星信号受信装置Ｒｘの速度ベクトルである。Δｔは、時刻ｋ−１と時刻ｋの間の時間である。数式８及び数式９は等価な式であり、計算処理のし易さに応じて選択する。
【００４９】
本発明の速度ベクトル誤差量の計算方法を図７に示す。
【００５０】
速度計算部４２は、速度ベクトルｖの計算で用いた数式２の時刻ｋにおける観測モデルＨｋを用いて、ＨＤＯＰｖも計算する（ステップＳ３）。ＨＤＯＰｖは、“Ｈｏｌｉｚｏｎｔａｌ  Ｄｉｌｕｔｉｏｎ  Ｏｆ  Ｐｒｅｃｉｓｉｏｎ”の頭文字をつなぎ合わせた言葉であり、測位衛星の水平面上での空間的なばらつき指標であることを意味する。ＤＯＰの計算方法は、ＧＰＳ分野では一般的であり、ＥＣＥＦ（Ｅａｒｔｈ−Ｃｅｎｔｅｒｅｄ、Ｅａｒｔｈ−Ｆｉｘｅｄ）座標系から、測位位置を座標の原点とするＥＮＵ（Ｌｏｃａｌ  Ｅａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系へと、観測モデルＨｋを座標系変換し、変換後の行列と転置後の行列を乗算し、乗算後の行列の逆行列を計算し、対角成分の目的要素を加算して平方根をとる。また、ＨＤＯＰｖの添え字のＶは、速度ベクトルｖのＨＤＯＰという意味である。つまり、ＨＤＯＰｖは、速度ベクトルｖを計算するために用いた測位衛星の水平面上での空間的なばらつき指標であることを意味する。
【００５１】
速度ベクトル誤差量計算部４４は、衛星測位に利用可能な測位衛星数が、所定閾値以下であるかどうか確認する（ステップＳ４）。所定閾値は、例えば、５衛星と定義される。
【００５２】
衛星測位に利用可能な測位衛星数が、所定閾値以下であるときは（ステップＳ４においてＹＥＳ）、速度ベクトル誤差量の計算方法が実行される（ステップＳ５、ステップＳ６）。衛星測位に利用可能な測位衛星数が、所定閾値より多いときは（ステップＳ４においてＮＯ）、速度ベクトル誤差量の計算方法が実行されない。
【００５３】
以下の説明では、衛星測位に利用可能な測位衛星数が、所定閾値以下であるとき（ステップＳ４においてＹＥＳ）を想定する。
【００５４】
速度ベクトル誤差量計算部４４は、以下に示す数式１０を用いて、｜ｓｉｎθｋ｜を計算する（ステップＳ５）。
【数１０】
0006546730.tif 000011
【００５５】
ｖ＾ｋ｜ｋは、推定状態量ｘ＾ｋ｜ｋから取り出した衛星信号受信装置Ｒｘの速度ベクトルである。εｋは、数式８又は数式９で計算された差分ベクトルである。｜・｜は絶対値の操作を表し、｜｜・｜｜は、ノルム計算の操作を表す。
【００５６】
速度ベクトル誤差量計算部４４は、以下に示す数式１１を用いて、速度ベクトル誤差量σｖを計算する（ステップＳ６）。ＨＤＯＰｖは、速度計算部４２がステップＳ３で計算したものである。σｃｖは、速度観測誤差定数（例えば、１．９ｍ／ｓ＝観測誤差１０Ｈｚ×搬送波の波長０．１９ｍ）である。
【数１１】
0006546730.tif 000012
【００５７】
ＨＤＯＰｖは、速度ベクトル誤差量σｖが大きくなり過ぎることを防ぐ目的で、上限値（例えば、１０）を設ける。速度ベクトル誤差量σｖは、図４又は図５の速度ベクトルｖの劣化度合いであり、｜ｓｉｎθｋ｜項の効果により、マルチパスが及ぼす影響が小さい図４の状況からマルチパスが及ぼす影響が大きい図５の状況へと近づくほど、大きくなる。
【００５８】
本発明の誤差楕円の計算方法を図８に示す。
【００５９】
誤差楕円計算部４５は、以下に示す数式１２を用いて、上述のプロセス雑音Ｑｋに代わる新たなプロセス雑音Ｗｋを計算する（ステップＳ７）。
【数１２】
0006546730.tif 000013
【００６０】
数式１２の右辺第１項は、差分ベクトル計算部４３がステップＳ２で計算した、εｋが関わる項である。数式１２の右辺第２項は、速度ベクトル誤差量計算部４４がステップＳ６で計算した、σｖが関わる項である。ｄｉａｇＬは、ＥＮＵ座標系での対角行列である。ＴＬ→Ｇは、ＥＮＵ座標系からＥＣＥＦ座標系への座標系変換行列である。
【００６１】
ここで、速度ベクトル誤差量σｖは、ＥＮＵ座標系の水平面（Ｅａｓｔ−Ｎｏｒｔｈ平面）で表されている。一方で、差分ベクトルεｋは、ＥＣＥＦ座標系で表されている。そこで、数式１２において、座標系を統一するため、速度ベクトル誤差量σｖに対して、ＥＮＵ座標系からＥＣＥＦ座標系への座標系変換を行なうのである。
【００６２】
なお、ステップＳ４においてＮＯとなり、速度ベクトル誤差量σｖの計算方法が実行されないときは、数式１２において速度ベクトル誤差量σｖはゼロになる。
【００６３】
誤差楕円計算部４５は、以下に示す数式１３〜１５を用いて、誤差楕円用の誤差共分散行列Ｐ〜ｋ｜ｋを計算する（ステップＳ８）。数式１３〜１５では、数式４、６、７と異なり、位置ｐを状態量とするが、速度ベクトルｖ、時計誤差及び時計誤差ドリフトを状態量としない。このため、数式１３では、数式４と異なり、数式１２で算出した新たなプロセス雑音Ｗｋを採用しており、数式４で採用したプロセス雑音Ｑｋを採用していない。また、数式１４及び数式１５のＨｋ,pとＲｋ,pは、数式２のＨｋとＲｋから位置ｐに関する成分のみを取り出した行列である。さらに、英文字の上部にチルダ（〜）を付けた行列は、誤差楕円用の行列であることを表す。
【数１３】
0006546730.tif 000014
【数１４】
0006546730.tif 000015
【数１５】
0006546730.tif 000016
【００６４】
誤差楕円計算部４５は、数式１６〜１９を用いて、誤差楕円の軸方向及び径の長さを計算する（ステップＳ９）。
【数１６】
0006546730.tif 000017
【数１７】
0006546730.tif 000018
【数１８】
0006546730.tif 000019
【数１９】
0006546730.tif 000020
【００６５】
σ２Ｅ、σ２ＥＮ、σ２ＥＵ、σ２ＮＥ、σ２Ｎ、σ２ＮＵ、σ２ＵＥ、σ２ＵＮ、σ２Ｕは、ＥＮＵ座標系での誤差楕円用の誤差共分散行列の各種成分である。ＴＧ→Ｌは、ＥＣＥＦ座標系からＥＮＵ座標系への座標系変換行列である。ｘは、測位位置を座標の原点とするＥＮＵ座標系での誤差楕円の円周である。ｒは、真位置が誤差楕円の内部に存在する確率に依存するパラメータであり、例えば、当該確率＝９５％であれば、ｒ＝２．４である。ｕＥ、ｕＮは、ｘの成分であり、数式１７で求められる（ｕＥ、ｕＮ）Ｔをプロットすることで、図１で示したようなＥＮＵ座標系のＥ−Ｎ平面での誤差楕円を描画できる。
【００６６】
λｉは、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円用の誤差共分散行列の固有値である。（Ｘｉ、Ｙｉ）Ｔは、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円用の誤差共分散行列の固有ベクトルである。数式１９によって、前記固有値λｉと前記固有ベクトル（Ｘｉ、Ｙｉ）Ｔから、誤差楕円の長軸及び短軸の軸方向並びに長半径及び短半径の長さを計算する。ここで、λ１≧λ２であるから、ｒ√λ１は、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円の長半径の長さとなり、ｒ√λ２は、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円の短半径の長さになる。
【００６７】
つまり、誤差楕円計算部４５は、数式１６〜１９を用いて、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円用の誤差共分散行列の固有ベクトル及び固有値を計算することにより、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円の軸方向及び径の長さを計算することができる。
【００６８】
誤差楕円計算部４５は、位置計算部４１での衛星信号受信装置Ｒｘの位置算出時にて、２次元測位演算が行なわれているかどうか確認する（ステップＳ１０）。２次元測位演算は、ＥＮＵ座標系での高さ方向（Ｕｐ軸）の位置を固定して位置計算する手法であり、測位衛星数が３衛星程度と少ない場合に実施される。２次元測位演算は、ＧＰＳ分野では一般的な手法であり、測位衛星数が少ない時に位置を求めることができるというメリットがあるが、求めた位置の誤差が固定する高さ方向の位置誤差に比例して大きくなるというデメリットもある。
【００６９】
２次元測位演算が行なわれるときは（ステップＳ１０においてＹＥＳ）、２次元測位演算時における誤差楕円の補正方法が実行される（ステップＳ１１）。２次元測位演算が行なわれないときは（ステップＳ１０においてＮＯ）、２次元測位演算時における誤差楕円の補正方法が実行されない。
【００７０】
以下の説明では、２次元測位演算が行なわれるとき（ステップＳ１０においてＹＥＳ）を想定する。
【００７１】
誤差楕円計算部４５は、数式２０を用いて、誤差楕円を補正する（ステップＳ１１）。
【数２０】
0006546730.tif 000021
数式２０の第１項は、数式１９のｒ√λｉそのものである。数式２０の第２項は、数式１９のｒ√λｉの補正項である。βは、測位衛星数に依存するパラメータである。位置計算部４１にて測位衛星数が少なく、２次元測位演算が行なわれたときは、βは１以上の値となり、位置計算部４１にて測位衛星数が多く、２次元測位演算が行なわれないときは、βは０近傍の値となり、つまり、βは測位衛星数に反比例する補正係数である。ｒは、数式１７に示したｒそのものである。σＵ２は、数式１６に示したσＵ２そのものである。
【００７２】
つまり、誤差楕円計算部４５は、位置計算部４１にて測位衛星数が少なく、２次元測位演算が行なわれたときに、水平面（Ｅ−Ｎ平面）での誤差楕円の長半径及び短半径の長さを、測位衛星数に反比例する補正係数βを乗じた高さ方向（Ｕｐ軸）の誤差分散で補正する。２次元測位演算では、高さ方向の位置を固定するという演算手法の特性上、高さ方向の位置誤差が水平面での位置誤差として現れるという特徴があるため、高さ方向の誤差分散を使用して、誤差楕円の水平面の長半径及び短半径を補正する。
【００７３】
本発明の誤差楕円の計算方法の全体構成を図９に示す。
【００７４】
ステップＳ２１の予測過程について説明する。この予測過程は、前回の更新過程で計算した推定状態量ｘ＾ｋ−１｜ｋ−１及び誤差共分散行列Ｐｋ−１｜ｋ−１を入力し、さらにプロセス雑音Ｑｋを入力する。次に、入力した各種情報をカルマンフィルタの予測過程の数式３、４に適用し、推定状態量ｘ＾ｋ｜ｋ−１及び誤差共分散行列Ｐｋ｜ｋ−１を計算し出力する。
【００７５】
ステップＳ２２の更新過程について説明する。この更新過程は、ステップＳ２１の予測過程で計算した推定状態量ｘ＾ｋ｜ｋ−１及び誤差共分散行列Ｐｋ｜ｋ−１を入力し、さらに観測量ｚｋ及び観測雑音Ｒｋを入力する。次に、入力した各種情報をカルマンフィルタの更新過程の数式５〜７に適用し、推定状態量ｘ＾ｋ｜ｋ及び誤差共分散行列Ｐｋ｜ｋを計算し出力する。
【００７６】
ステップＳ２３の計算過程について説明する。この計算過程は、ステップＳ２１の予測過程で計算した推定状態量ｘ＾ｋ｜ｋ−１を入力し、さらにステップＳ２２の更新過程で計算した推定状態量ｘ＾ｋ｜ｋを入力する。入力した各種推定状態量から衛星信号受信装置Ｒｘの位置p ＾ｋ｜ｋ、 p ＾ｋ｜ｋ-1及び速度ベクトルｖ ＾ｋ｜ｋを取り出す。次に、取り出した各種情報を数式８〜１２に適用し、εｋ、｜ｓｉｎθｋ｜、σｖ及び上述のプロセス雑音Ｑｋに代わる新たなプロセス雑音Ｗｋを計算し、プロセス雑音Ｗｋを出力する。
【００７７】
ステップＳ２４の誤差楕円の予測過程について説明する。この誤差楕円の予測過程は、前回の更新過程で計算した誤差楕円用の誤差共分散行列Ｐ~ｋ-1｜ｋ−１を入力し、さらにステップＳ２３の計算過程で計算した新たなプロセス雑音Ｗｋを入力する。次に、入力した情報を数式１３に適用し、誤差楕円用の誤差共分散行列Ｐ~ｋ｜ｋ−１を計算し出力する。
【００７８】
ステップＳ２５の誤差楕円の更新過程について説明する。この誤差楕円の更新過程は、ステップＳ２４の誤差楕円の予測過程で計算した誤差楕円用の誤差共分散行列Ｐ~ｋ｜ｋ−１を入力し、さらに観測雑音Ｒｋを入力する。入力した観測雑音Ｒｋから位置ｐに関する成分のみを取り出す。次に、入力した情報及び取り出した情報を数式１４、１５に適用し、誤差楕円用の誤差共分散行列Ｐ~ｋ｜ｋを計算し出力する。
【００７９】
ステップＳ２６の誤差楕円の軸方向及び径の長さの計算過程について説明する。ステップＳ２６の計算過程は、ステップＳ２５の誤差楕円の再更新過程で計算した誤差楕円用の誤差共分散行列Ｐ~ｋ｜ｋを入力する。次に、入力した情報を数式１６〜１９に適用し、誤差楕円の軸方向及び径の長さを計算する。また、２次元測位演算が行われていれば、数式２０によって、誤差楕円の径の長さを補正する。計算された誤差楕円の軸方向及び径の長さは、衛星信号受信装置Ｒｘを使用するユーザ側に外部出力される。
【００８０】
本発明及び比較例の誤差楕円の時系列グラフを図１０に示す。図１０に示した時系列グラフは、衛星信号受信装置Ｒｘの搭載車両が図２のような高架下を走行時の、誤差楕円の時系列グラフである。「位置誤差」は、測位位置と真位置間の距離を示す。「従来誤差楕円」は、従来技術の誤差楕円の長半径を示す。「本発明誤差楕円」は、本発明の誤差楕円の長半径を示す。
【００８１】
誤差楕円は、以下の条件を満たすことが理想的である。（１）誤差楕円は位置誤差を上回ること。誤差楕円が位置誤差を下回れば、測位位置の誤差が大きいにも関わらず、ユーザは測位位置の誤差が実際より小さいと認識してしまう。（２）誤差楕円は位置誤差を過剰に上回らず少しだけ上回ること。誤差楕円が位置誤差を過剰に上回れば、測位位置の誤差が小さいにも関わらず、ユーザは測位位置の誤差が実際より大きすぎると認識してしまう。
【００８２】
従来の誤差楕円は、測位時刻２３：１７：３７以降において、位置誤差を下回っているため、上述の条件を満たしておらず理想的ではない。本発明の誤差楕円は、測位時刻２３：１７：３７以降を含めて、時系列の全測位時刻において、位置誤差を少しだけ上回っているため、上述の条件を満たしており理想的である。
【００８３】
以上に説明の事項は、以下のようにまとめられる。
【００８４】
本発明では、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルεｋが長くなり、差分ベクトルεｋが長いほど、差分ベクトルεｋを用いて算出したプロセス雑音Ｗｋが大きくなり、プロセス雑音Ｗｋを用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。
【００８５】
そして、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルεｋと速度ベクトルｖ＾ｋ｜ｋのなすベクトル間角度θｋが９０°又が２７０°に近くなり、ベクトル間角度θｋが９０°又が２７０°に近いほど、ベクトル間角度θｋを用いて算出した速度ベクトル誤差量σｖが大きくなり、速度ベクトル誤差量σｖを用いて算出したプロセス雑音Ｗｋが大きくなり、プロセス雑音Ｗｋを用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。
【００８６】
さらに、測位演算に利用可能な測位衛星数が少なく、ＥＮＵ座標系の高さ方向（Ｕｐ軸）の測位位置を固定する２次元測位演算を行なうときに、ＥＮＵ座標系の水平面（Ｅ−Ｎ平面）の誤差楕円の長半径及び短半径を、ＥＮＵ座標系の高さ方向（Ｕｐ軸）における誤差楕円の誤差分散で補正することで、誤差楕円を大きくすることができる。逆に、測位演算に利用可能な測位衛星数が多く、ＥＮＵ座標系の高さ方向（Ｕｐ軸）の測位位置を固定しない３次元測位演算を行なうときは、ＥＮＵ座標系の水平面（Ｅ−Ｎ平面）の誤差楕円の長半径及び短半径を、ＥＮＵ座標系の高さ方向（Ｕｐ軸）における誤差楕円の誤差分散で補正しないため、誤差楕円を大きくし過ぎないようにすることができる。
【産業上の利用可能性】
【００８７】
本発明の衛星信号受信装置は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。
【符号の説明】
【００８８】
Ｒｘ：衛星信号受信装置
Ｓ１、Ｓ２、Ｓ３、Ｓ４：測位衛星
Ｅｌ：高架構造
Ｒｆ：反射物体
１：信号受信部
２：追尾処理部
３：復調処理部
４：測位演算部
２１：擬似距離観測部
２２：ドップラー周波数観測部
２３：航法データ抽出部
４１：位置計算部
４２：速度計算部
４３：差分ベクトル計算部
４４：速度ベクトル誤差量計算部
４５：誤差楕円計算部
【図１】
0006546730.tif 000022
【図２】
0006546730.tif 000023
【図３】
0006546730.tif 000024
【図４】
0006546730.tif 000025
【図５】
0006546730.tif 000026
【図６】
0006546730.tif 000027
【図７】
0006546730.tif 000028
【図８】
0006546730.tif 000029
【図９】
0006546730.tif 000030
【図１０】
0006546730.tif 000031