JP6546730B2

JP6546730B2 - 衛星信号受信装置

Info

Publication number: JP6546730B2
Application number: JP2014207790A
Authority: JP
Inventors: 真嗣小田
Original assignee: Japan Radio Co Ltd
Current assignee: Japan Radio Co Ltd
Priority date: 2014-10-09
Filing date: 2014-10-09
Publication date: 2019-07-17
Anticipated expiration: 2034-10-09
Also published as: JP2016075646A

Description

本発明は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対して反映させる技術に関する。

ＧＰＳ（米国）、Ｇａｌｉｌｅｏ（欧州）、ＧＬＯＮＡＳＳ（ロシア）及びＢｅｉＤｏｕ（中国）等の衛星航法システムによる衛星信号を用いて測位演算を行なう、衛星信号受信装置（ＧＮＳＳ（ＧｌｏｂａｌＮａｖｉｇａｔｉｏｎＳａｔｅｌｌｉｔｅＳｙｓｔｅｍ）受信装置）が、従来から広くユーザに利用されている。

衛星信号受信装置が計算する測位位置は、衛星位置誤差、衛星時計誤差及びマルチパス等の影響による測位誤差を含んでいる。衛星信号受信装置のユーザは、測位位置の確度情報を必要とするため、衛星信号受信装置の測位演算部は、測位位置の精度指標を示す誤差楕円を計算する。

測位位置及び誤差楕円の表示方法を図１に示す。ＥＮＵ（ＬｏｃａｌＥａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系は、測位位置を座標の原点とする座標系である。誤差楕円は、測位位置に中心を有し、真位置が楕円の内部に存在する確率が所定の確率（例えば、９５％）であるような軸及び径を有する。

特開平０５−３３３１３１号公報

特許文献１では、誤差楕円の計算方法が開示されている。この文献では、複数の測位衛星の空間的なばらつき指標を示すＤＯＰ（ＤｉｌｕｔｉｏｎＯｆＰｒｅｃｉｓｉｏｎ）に基づいて、誤差楕円を計算している。つまり、測位演算に利用可能な測位衛星数が少ないときや、測位演算に利用可能な測位衛星の配置が偏っているときに、誤差楕円の径が大きく計算される。しかし、図２を用いて以下に説明するように、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対して反映させることができなかった。

マルチパスが測位位置に及ぼす影響を図２に示す。衛星信号受信装置Ｒｘは、高架構造Ｅｌの下部に位置する。測位衛星Ｓ１、Ｓ３、Ｓ４は、衛星信号受信装置Ｒｘから見て天頂以外の方向に位置する。測位衛星Ｓ２は、衛星信号受信装置Ｒｘから見て天頂方向に位置する。反射物体Ｒｆは、衛星信号を反射する物体である。

すると、測位衛星Ｓ１、Ｓ３、Ｓ４からの衛星信号は、高架構造Ｅｌに遮られず、衛星信号受信装置Ｒｘにより受信される。よって、衛星信号受信装置Ｒｘと測位衛星Ｓ１、Ｓ３、Ｓ４の間の擬似距離は、誤差をあまり含んでいない。

しかし、測位衛星Ｓ２からの衛星信号は、高架構造Ｅｌに遮られて、反射物体Ｒｆで反射されるマルチパスを経てから、衛星信号受信装置Ｒｘにより受信される。よって、衛星信号受信装置Ｒｘと測位衛星Ｓ２の間の擬似距離は、誤差を大きく含んでおり、擬似距離に誤差が重畳すれば、擬似距離を使用して計算された測位位置の精度も劣化する。

このように、マルチパスの影響があるときには、マルチパスの影響がないときと比べて、測位位置の精度が劣化するため、誤差楕円の径を大きく計算する必要がある。しかし、マルチパスの影響があるときでも、マルチパスの影響がないときと比べて、複数の測位衛星の空間的なばらつき指標を示すＤＯＰが大きくなるわけではない。よって、特許文献１では、マルチパスの影響があるときでも、マルチパスの影響がないときと比べて、ＤＯＰに基づいた誤差楕円の径は大きく計算されない。つまり、この文献での誤差楕円は、本来あるべき誤差楕円から乖離していた。

そこで、前記課題を解決するために、本発明は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることを目的とする。

上記目的を達成するために、測位時刻より過去の時刻において予測した、当該測位時刻における衛星信号受信装置の予測位置と、当該測位時刻において更新した、当該測位時刻における衛星信号受信装置の測位位置と、の差分ベクトルが長いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算することとした。

具体的には、本発明は、測位演算を行なう衛星信号受信装置であって、自装置と測位衛星の間の擬似距離を観測する擬似距離観測部と、前記擬似距離観測部が観測した擬似距離に基づいて、測位時刻より過去の時刻において、当該測位時刻における自装置の予測位置を予測した後に、当該測位時刻において、当該測位時刻における自装置の測位位置を更新する位置計算部と、前記位置計算部が当該測位時刻より過去の時刻において予測した、当該測位時刻における自装置の予測位置と、前記位置計算部が当該測位時刻において更新した、当該測位時刻における自装置の測位位置と、の差分ベクトルを計算する差分ベクトル計算部と、前記差分ベクトル計算部が計算した差分ベクトルが長いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算する誤差楕円計算部と、を備えることを特徴とする衛星信号受信装置である。

この構成によれば、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルが長くなり、差分ベクトルが長いほど、差分ベクトルを用いて算出したプロセス雑音が大きくなり、プロセス雑音を用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。

また、本発明は、自装置と測位衛星の間のドップラー周波数を観測するドップラー周波数観測部と、前記ドップラー周波数観測部が観測したドップラー周波数に基づいて、当該測位時刻において、当該測位時刻における自装置の速度ベクトルを計算する速度計算部と、前記差分ベクトル計算部が計算した差分ベクトルと、前記速度計算部が計算した速度ベクトルと、のなすベクトル間角度を計算し、計算したベクトル間角度に基づいて、速度ベクトル誤差量を計算する速度ベクトル誤差量計算部と、をさらに備え、前記誤差楕円計算部は、前記速度ベクトル誤差量計算部が計算したベクトル間角度が９０°又は２７０°に近いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算することを特徴とする衛星信号受信装置である。

この構成によれば、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルと速度ベクトルのなすベクトル間角度が９０°又は２７０°に近くなり、ベクトル間角度が９０°又は２７０°に近いほど、ベクトル間角度を用いて算出した速度ベクトル誤差量が大きくなり、速度ベクトル誤差量を用いて算出したプロセス雑音が大きくなり、プロセス雑音を用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。

また、本発明の前記誤差楕円計算部は、ＥＮＵ（ＬｏｃａｌＥａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系の高さ（Ｕｐ軸）方向における自装置の測位位置を固定する２次元測位演算が行なわれるときに、ＥＮＵ座標系の水平面（Ｅａｓｔ−Ｎｏｒｔｈ平面）での誤差楕円の長半径及び短半径を、ＥＮＵ座標系の高さ方向における誤差楕円の誤差分散で補正することを特徴とする衛星信号受信装置である。

この構成によれば、測位演算に利用可能な測位衛星数が少なく、高さ方向の測位位置を固定する２次元測位演算を行なうときに、誤差楕円を大きくすることができる。逆に、測位演算に利用可能な測位衛星数が多く、高さ方向の測位位置を固定しない３次元測位演算を行なうときは、誤差楕円を大きくし過ぎないようにすることができる。

このように、本発明は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。

測位位置及び誤差楕円の表示方法を示す図である。マルチパスが測位位置に及ぼす影響を示す図である。本発明の衛星信号受信装置の構成を示す図である。本発明の差分ベクトル及びベクトル間角度を示す図である。本発明の差分ベクトル及びベクトル間角度を示す図である。本発明の差分ベクトルの計算方法を示す図である。本発明の速度ベクトル誤差量の計算方法を示す図である。本発明の誤差楕円の計算方法を示す図である。本発明の誤差楕円の計算方法の全体構成を示す図である。本発明及び比較例の誤差楕円の時系列グラフである。

添付の図面を参照して本発明の実施形態を説明する。以下に説明する実施形態は本発明の実施の例であり、本発明は以下の実施形態に制限されるものではない。なお、本明細書及び図面において符号が同じ構成要素は、相互に同一のものを示すものとする。

本発明の衛星信号受信装置の構成を図３に示す。衛星信号受信装置Ｒｘは、信号受信部１、追尾処理部２、復調処理部３及び測位演算部４から構成される。

信号受信部１は、アンテナを介してＧＮＳＳ信号を受信する。

追尾処理部２は、擬似距離観測部２１、ドップラー周波数観測部２２及び航法データ抽出部２３から構成される。

擬似距離観測部２１は、衛星信号受信装置Ｒｘと測位衛星の間の擬似距離を観測する。ドップラー周波数観測部２２は、衛星信号受信装置Ｒｘと測位衛星の間のドップラー周波数を観測する。航法データ抽出部２３は、ＧＮＳＳ信号から航法データのビット情報を抽出する。

復調処理部３は、航法データのビット情報を復調（又は復号）し、測位演算に必要なエフェメリス及び衛星時計情報を摘出して出力する。

追尾処理部２及び復調処理部３は、最大可視衛星数に応じて各々必要な個数が定まり、１衛星に対して各々１個を用意する必要がある。

測位演算部４は、位置計算部４１、速度計算部４２、差分ベクトル計算部４３、速度ベクトル誤差量計算部４４及び誤差楕円計算部４５から構成される。

位置計算部４１は、擬似距離観測部２１が観測した擬似距離、復調処理部３が出力したエフェメリス及び衛星時計情報、並びに、衛星信号受信装置Ｒｘが有する受信装置時計情報に基づいて、測位時刻より過去の時刻において、当該測位時刻における衛星信号受信装置Ｒｘの予測位置を予測した後に、当該測位時刻において、当該測位時刻における衛星信号受信装置Ｒｘの測位位置を更新する。

速度計算部４２は、ドップラー周波数観測部２２が観測したドップラー周波数、復調処理部３が出力したエフェメリス及び衛星時計情報、並びに、衛星信号受信装置Ｒｘが有する受信装置時計情報に基づいて、当該測位時刻において、当該測位時刻における衛星信号受信装置Ｒｘの速度ベクトルを計算する。

差分ベクトル計算部４３は、位置計算部４１が当該測位時刻より過去の時刻において予測した、当該測位時刻における衛星信号受信装置Ｒｘの予測位置と、位置計算部４１が当該測位時刻において更新した、当該測位時刻における衛星信号受信装置Ｒｘの測位位置と、の差分ベクトルを計算する。

速度ベクトル誤差量計算部４４は、差分ベクトル計算部４３が計算した差分ベクトルと、速度計算部４２が計算した速度ベクトルと、のなすベクトル間角度を計算し、計算したベクトル間角度に基づいて、速度ベクトル誤差量を計算する。

誤差楕円計算部４５は、差分ベクトル計算部４３が計算した差分ベクトルと、速度ベクトル誤差量計算部４４が計算した速度ベクトル誤差量と、を用いてプロセス雑音を計算し、計算したプロセス雑音と、位置計算部４１で使用した情報と、を用いて誤差楕円の径の長さ及び誤差楕円の軸の方向を計算する。

誤差楕円計算部４５は、差分ベクトル計算部４３が計算した差分ベクトルが長いほど、かつ、速度ベクトル誤差量計算部４４が計算したベクトル間角度が９０°又は２７０°に近いほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算する。

測位演算部４は、位置計算部４１が計算した測位位置と、誤差楕円計算部４５が計算した誤差楕円の径の長さ（長半径、短半径）及び誤差楕円の軸の方向と、を外部出力する。外部出力された情報は、衛星信号受信装置Ｒｘを使用するユーザが、ナビゲーション等に利用する。

本発明の差分ベクトル及びベクトル間角度を図４、５に示す。ｐ^＾ _{ｋ｜ｋ−１}は、位置計算部４１が時刻ｋ−１において予測した、時刻ｋにおける衛星信号受信装置Ｒｘの予測位置である。ｐ^＾ _ｋ｜ｋは、位置計算部４１が時刻ｋにおいて更新した、時刻ｋにおける衛星信号受信装置Ｒｘの測位位置である。ｖ^＾ _ｋ｜ｋは、速度計算部４２が時刻ｋにおいて計算した、時刻ｋにおける衛星信号受信装置Ｒｘの速度ベクトルである。Δｔは、時刻ｋ−１と時刻ｋの間の時間であり、ｐ^＾ _{ｋ｜ｋ−１}は、ｐ^＾ _{ｋ−１｜ｋ−１}＋ｖ^＾ _{ｋ−１｜ｋ−１}Δｔと表わせる。

ε_ｋは、ｐ^＾ _{ｋ｜ｋ−１}とｐ_ｋ｜ｋの差分ベクトルである。θ_ｋは、ε_ｋとｖ^＾ _ｋ｜ｋのなすベクトル間角度である。図４では、ε_ｋは有限の長さを有し、θ_ｋは０°又は１８０°に近く、図５では、ε_ｋは有限の長さを有し、θ_ｋは９０°又は２７０°に近い。差分ベクトルε_ｋの長さは、予測位置ｐ^＾ _{ｋ｜ｋ−１}と実計算した位置ｐ^＾ _ｋ｜ｋとの乖離度合いを示しており、マルチパスが及ぼす影響度合いが大きい状況であれば長くなる。ベクトル間角度θ_ｋは、道路の両脇に高層ビルがある場合や高架下走行時のような、マルチパスが及ぼす影響の度合いが非常に強い状況下では、９０°又は２７０°に近くなる。

ここで、一般的に、進行方向には障害物はなく、進行方向と直角（９０°又は２７０°）の方向に障害物があるため、マルチパスが及ぼす影響を受けている状況下では、図４の状態ではなく、図５の状態になる。つまり、図５に示した状況は、図４に示した状況より、マルチパスが及ぼす影響度合いが強い状況である。

測位演算部４を構成する位置計算部４１、速度計算部４２、差分ベクトル計算部４３、速度ベクトル誤差量計算部４４及び誤差楕円計算部４５の処理の流れを図６〜８に示す。

本発明の差分ベクトルの計算方法を図６に示す。

位置計算部４１及び速度計算部４２は、以下に示す数式１〜７を用いて、ｐ^＾ _ｋ｜ｋ及びｖ^＾ _ｋ｜ｋをそれぞれ計算する（ステップＳ１）。位置計算部４１及び速度計算部４２は、状態方程式及び観測方程式をそれぞれ立式し、カルマンフィルタによる手法で、観測量である擬似距離情報及びドップラー周波数情報を用いて、衛星信号受信装置Ｒｘの測位位置ｐ、速度ベクトルｖ、時計誤差及び時計誤差ドリフトを算出する。状態量ｘは、測位位置ｐ、速度ベクトルｖ、時計誤差及び時計誤差ドリフトを構成成分とし、観測量ｚは、擬似距離の観測量及びドップラー周波数の観測量を構成成分とする。

カルマンフィルタに適用する状態方程式は、数式１のように表わされる。

ｘ_ｋは、時刻ｋにおける状態量である。Ｆ_ｋは、時刻ｋにおけるシステムモデルである。ω_ｋは、平均が０、共分散行列がＱ_ｋである多変数正規分布に従うプロセス雑音である。

カルマンフィルタに適用する観測方程式は、数式２のように表わされる。

ｚ_ｋは、時刻ｋにおける観測量である。Ｈ_ｋは、時刻ｋにおける観測モデルである。ｕ_ｋは、平均が０、共分散行列がＲ_ｋである多変数正規分布に従う観測雑音である。

状態方程式と観測方程式を用いてカルマンフィルタの予測過程と更新過程を計算する。

カルマンフィルタの予測過程は、数式３、４のように表わされる。

ｘ^＾ _{ｋ｜ｋ−１}は、時刻ｋ−１において予測された、時刻ｋにおける推定状態量である。ｘ^＾ _{ｋ−１｜ｋ−１}は、時刻ｋ−１において更新された、時刻ｋ−１における推定状態量である。Ｐ_{ｋ｜ｋ−１}は、ｃｏｖ（ｘ_ｋ−ｘ^＾ _{ｋ｜ｋ−１}）から分かるように、時刻ｋ−１において予測された、時刻ｋにおける推定状態量の誤差共分散行列である。Ｐ_{ｋ−１｜ｋ−１}は、時刻ｋ−１において更新された、時刻ｋ−１における推定状態量の誤差共分散行列である。

カルマンフィルタの更新過程は、数式５〜７のように表わされる。

ｘ^＾ _ｋ｜ｋは、時刻ｋにおいて更新された、時刻ｋにおける推定状態量である。Ｐ_ｋ｜ｋは、ｃｏｖ（ｘ_ｋ−ｘ^＾ _ｋ｜ｋ）から分かるように、時刻ｋにおいて更新された、時刻ｋにおける推定状態量の誤差共分散行列である。Ｋ_ｋは、時刻ｋにおける最適カルマンゲインである。

差分ベクトル計算部４３は、以下に示す数式８又は数式９を用いて、差分ベクトルε_ｋを計算する（ステップＳ２）。

p^＾ _ｋ｜ｋは、カルマンフィルタの更新過程で計算された推定状態量ｘ^＾ _ｋ｜ｋから取り出した衛星信号受信装置Ｒｘの測位位置である。p^＾ _ｋ｜ｋ-1は、推定状態量ｘ^＾ _ｋ｜ｋ-1から取り出した衛星信号受信装置Ｒｘの予測位置である。p^＾ _{ｋ-1｜ｋ-1}は、推定状態量ｘ^＾ _{ｋ-1｜ｋ-1}から取り出した衛星信号受信装置Ｒｘの測位位置である。ｖ^＾ _{ｋ-1｜ｋ-1}は、推定状態量ｘ^＾ _{ｋ-1｜ｋ-1}から取り出した衛星信号受信装置Ｒｘの速度ベクトルである。Δｔは、時刻ｋ−１と時刻ｋの間の時間である。数式８及び数式９は等価な式であり、計算処理のし易さに応じて選択する。

本発明の速度ベクトル誤差量の計算方法を図７に示す。

速度計算部４２は、速度ベクトルｖの計算で用いた数式２の時刻ｋにおける観測モデルＨ_ｋを用いて、ＨＤＯＰ_ｖも計算する（ステップＳ３）。ＨＤＯＰ_ｖは、“ＨｏｌｉｚｏｎｔａｌＤｉｌｕｔｉｏｎＯｆＰｒｅｃｉｓｉｏｎ”の頭文字をつなぎ合わせた言葉であり、測位衛星の水平面上での空間的なばらつき指標であることを意味する。ＤＯＰの計算方法は、ＧＰＳ分野では一般的であり、ＥＣＥＦ（Ｅａｒｔｈ−Ｃｅｎｔｅｒｅｄ、Ｅａｒｔｈ−Ｆｉｘｅｄ）座標系から、測位位置を座標の原点とするＥＮＵ（ＬｏｃａｌＥａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系へと、観測モデルＨ_ｋを座標系変換し、変換後の行列と転置後の行列を乗算し、乗算後の行列の逆行列を計算し、対角成分の目的要素を加算して平方根をとる。また、ＨＤＯＰ_ｖの添え字のＶは、速度ベクトルｖのＨＤＯＰという意味である。つまり、ＨＤＯＰ_ｖは、速度ベクトルｖを計算するために用いた測位衛星の水平面上での空間的なばらつき指標であることを意味する。

速度ベクトル誤差量計算部４４は、衛星測位に利用可能な測位衛星数が、所定閾値以下であるかどうか確認する（ステップＳ４）。所定閾値は、例えば、５衛星と定義される。

衛星測位に利用可能な測位衛星数が、所定閾値以下であるときは（ステップＳ４においてＹＥＳ）、速度ベクトル誤差量の計算方法が実行される（ステップＳ５、ステップＳ６）。衛星測位に利用可能な測位衛星数が、所定閾値より多いときは（ステップＳ４においてＮＯ）、速度ベクトル誤差量の計算方法が実行されない。

以下の説明では、衛星測位に利用可能な測位衛星数が、所定閾値以下であるとき（ステップＳ４においてＹＥＳ）を想定する。

速度ベクトル誤差量計算部４４は、以下に示す数式１０を用いて、｜ｓｉｎθ_ｋ｜を計算する（ステップＳ５）。

ｖ^＾ _ｋ｜ｋは、推定状態量ｘ^＾ _ｋ｜ｋから取り出した衛星信号受信装置Ｒｘの速度ベクトルである。ε_ｋは、数式８又は数式９で計算された差分ベクトルである。｜・｜は絶対値の操作を表し、｜｜・｜｜は、ノルム計算の操作を表す。

速度ベクトル誤差量計算部４４は、以下に示す数式１１を用いて、速度ベクトル誤差量σ_ｖを計算する（ステップＳ６）。ＨＤＯＰ_ｖは、速度計算部４２がステップＳ３で計算したものである。σ_ｃｖは、速度観測誤差定数（例えば、１．９ｍ／ｓ＝観測誤差１０Ｈｚ×搬送波の波長０．１９ｍ）である。

ＨＤＯＰ_ｖは、速度ベクトル誤差量σ_ｖが大きくなり過ぎることを防ぐ目的で、上限値（例えば、１０）を設ける。速度ベクトル誤差量σ_ｖは、図４又は図５の速度ベクトルｖの劣化度合いであり、｜ｓｉｎθ_ｋ｜項の効果により、マルチパスが及ぼす影響が小さい図４の状況からマルチパスが及ぼす影響が大きい図５の状況へと近づくほど、大きくなる。

本発明の誤差楕円の計算方法を図８に示す。

誤差楕円計算部４５は、以下に示す数式１２を用いて、上述のプロセス雑音Ｑ_ｋに代わる新たなプロセス雑音Ｗ_ｋを計算する（ステップＳ７）。

数式１２の右辺第１項は、差分ベクトル計算部４３がステップＳ２で計算した、ε_ｋが関わる項である。数式１２の右辺第２項は、速度ベクトル誤差量計算部４４がステップＳ６で計算した、σ_ｖが関わる項である。ｄｉａｇ_Ｌは、ＥＮＵ座標系での対角行列である。Ｔ_Ｌ→Ｇは、ＥＮＵ座標系からＥＣＥＦ座標系への座標系変換行列である。

ここで、速度ベクトル誤差量σ_ｖは、ＥＮＵ座標系の水平面（Ｅａｓｔ−Ｎｏｒｔｈ平面）で表されている。一方で、差分ベクトルε_ｋは、ＥＣＥＦ座標系で表されている。そこで、数式１２において、座標系を統一するため、速度ベクトル誤差量σ_ｖに対して、ＥＮＵ座標系からＥＣＥＦ座標系への座標系変換を行なうのである。

なお、ステップＳ４においてＮＯとなり、速度ベクトル誤差量σ_ｖの計算方法が実行されないときは、数式１２において速度ベクトル誤差量σ_ｖはゼロになる。

誤差楕円計算部４５は、以下に示す数式１３〜１５を用いて、誤差楕円用の誤差共分散行列Ｐ^〜 _ｋ｜ｋを計算する（ステップＳ８）。数式１３〜１５では、数式４、６、７と異なり、位置ｐを状態量とするが、速度ベクトルｖ、時計誤差及び時計誤差ドリフトを状態量としない。このため、数式１３では、数式４と異なり、数式１２で算出した新たなプロセス雑音Ｗ_ｋを採用しており、数式４で採用したプロセス雑音Ｑ_ｋを採用していない。また、数式１４及び数式１５のＨ_ｋ,pとＲ_ｋ,pは、数式２のＨ_ｋとＲ_ｋから位置ｐに関する成分のみを取り出した行列である。さらに、英文字の上部にチルダ（〜）を付けた行列は、誤差楕円用の行列であることを表す。

誤差楕円計算部４５は、数式１６〜１９を用いて、誤差楕円の軸方向及び径の長さを計算する（ステップＳ９）。

σ^２ _Ｅ、σ^２ _ＥＮ、σ^２ _ＥＵ、σ^２ _ＮＥ、σ^２ _Ｎ、σ^２ _ＮＵ、σ^２ _ＵＥ、σ^２ _ＵＮ、σ^２ _Ｕは、ＥＮＵ座標系での誤差楕円用の誤差共分散行列の各種成分である。Ｔ_Ｇ→Ｌは、ＥＣＥＦ座標系からＥＮＵ座標系への座標系変換行列である。ｘは、測位位置を座標の原点とするＥＮＵ座標系での誤差楕円の円周である。ｒは、真位置が誤差楕円の内部に存在する確率に依存するパラメータであり、例えば、当該確率＝９５％であれば、ｒ＝２．４である。ｕ_Ｅ、ｕ_Ｎは、ｘの成分であり、数式１７で求められる（ｕ_Ｅ、ｕ_Ｎ）^Ｔをプロットすることで、図１で示したようなＥＮＵ座標系のＥ−Ｎ平面での誤差楕円を描画できる。

λ_ｉは、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円用の誤差共分散行列の固有値である。（Ｘ_ｉ、Ｙ_ｉ）^Ｔは、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円用の誤差共分散行列の固有ベクトルである。数式１９によって、前記固有値λ_ｉと前記固有ベクトル（Ｘ_ｉ、Ｙ_ｉ）^Ｔから、誤差楕円の長軸及び短軸の軸方向並びに長半径及び短半径の長さを計算する。ここで、λ_１≧λ_２であるから、ｒ√λ_１は、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円の長半径の長さとなり、ｒ√λ_２は、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円の短半径の長さになる。

つまり、誤差楕円計算部４５は、数式１６〜１９を用いて、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円用の誤差共分散行列の固有ベクトル及び固有値を計算することにより、ＥＮＵ座標系のＥ−Ｎ平面での誤差楕円の軸方向及び径の長さを計算することができる。

誤差楕円計算部４５は、位置計算部４１での衛星信号受信装置Ｒｘの位置算出時にて、２次元測位演算が行なわれているかどうか確認する（ステップＳ１０）。２次元測位演算は、ＥＮＵ座標系での高さ方向（Ｕｐ軸）の位置を固定して位置計算する手法であり、測位衛星数が３衛星程度と少ない場合に実施される。２次元測位演算は、ＧＰＳ分野では一般的な手法であり、測位衛星数が少ない時に位置を求めることができるというメリットがあるが、求めた位置の誤差が固定する高さ方向の位置誤差に比例して大きくなるというデメリットもある。

２次元測位演算が行なわれるときは（ステップＳ１０においてＹＥＳ）、２次元測位演算時における誤差楕円の補正方法が実行される（ステップＳ１１）。２次元測位演算が行なわれないときは（ステップＳ１０においてＮＯ）、２次元測位演算時における誤差楕円の補正方法が実行されない。

以下の説明では、２次元測位演算が行なわれるとき（ステップＳ１０においてＹＥＳ）を想定する。

誤差楕円計算部４５は、数式２０を用いて、誤差楕円を補正する（ステップＳ１１）。

数式２０の第１項は、数式１９のｒ√λ_ｉそのものである。数式２０の第２項は、数式１９のｒ√λ_ｉの補正項である。βは、測位衛星数に依存するパラメータである。位置計算部４１にて測位衛星数が少なく、２次元測位演算が行なわれたときは、βは１以上の値となり、位置計算部４１にて測位衛星数が多く、２次元測位演算が行なわれないときは、βは０近傍の値となり、つまり、βは測位衛星数に反比例する補正係数である。ｒは、数式１７に示したｒそのものである。σ_Ｕ ^２は、数式１６に示したσ_Ｕ ^２そのものである。

つまり、誤差楕円計算部４５は、位置計算部４１にて測位衛星数が少なく、２次元測位演算が行なわれたときに、水平面（Ｅ−Ｎ平面）での誤差楕円の長半径及び短半径の長さを、測位衛星数に反比例する補正係数βを乗じた高さ方向（Ｕｐ軸）の誤差分散で補正する。２次元測位演算では、高さ方向の位置を固定するという演算手法の特性上、高さ方向の位置誤差が水平面での位置誤差として現れるという特徴があるため、高さ方向の誤差分散を使用して、誤差楕円の水平面の長半径及び短半径を補正する。

本発明の誤差楕円の計算方法の全体構成を図９に示す。

ステップＳ２１の予測過程について説明する。この予測過程は、前回の更新過程で計算した推定状態量ｘ^＾ _{ｋ−１｜ｋ−１}及び誤差共分散行列Ｐ_{ｋ−１｜ｋ−１}を入力し、さらにプロセス雑音Ｑ_ｋを入力する。次に、入力した各種情報をカルマンフィルタの予測過程の数式３、４に適用し、推定状態量ｘ^＾ _{ｋ｜ｋ−１}及び誤差共分散行列Ｐ_{ｋ｜ｋ−１}を計算し出力する。

ステップＳ２２の更新過程について説明する。この更新過程は、ステップＳ２１の予測過程で計算した推定状態量ｘ^＾ _{ｋ｜ｋ−１}及び誤差共分散行列Ｐ_{ｋ｜ｋ−１}を入力し、さらに観測量ｚ_ｋ及び観測雑音Ｒ_ｋを入力する。次に、入力した各種情報をカルマンフィルタの更新過程の数式５〜７に適用し、推定状態量ｘ^＾ _ｋ｜ｋ及び誤差共分散行列Ｐ_ｋ｜ｋを計算し出力する。

ステップＳ２３の計算過程について説明する。この計算過程は、ステップＳ２１の予測過程で計算した推定状態量ｘ^＾ _{ｋ｜ｋ−１}を入力し、さらにステップＳ２２の更新過程で計算した推定状態量ｘ^＾ _ｋ｜ｋを入力する。入力した各種推定状態量から衛星信号受信装置Ｒｘの位置p^＾ _ｋ｜ｋ、 p^＾ _ｋ｜ｋ-1及び速度ベクトルｖ^＾ _ｋ｜ｋを取り出す。次に、取り出した各種情報を数式８〜１２に適用し、ε_ｋ、｜ｓｉｎθ_ｋ｜、σ_ｖ及び上述のプロセス雑音Ｑ_ｋに代わる新たなプロセス雑音Ｗ_ｋを計算し、プロセス雑音Ｗ_ｋを出力する。

ステップＳ２４の誤差楕円の予測過程について説明する。この誤差楕円の予測過程は、前回の更新過程で計算した誤差楕円用の誤差共分散行列Ｐ~_{ｋ-1｜ｋ−１}を入力し、さらにステップＳ２３の計算過程で計算した新たなプロセス雑音Ｗ_ｋを入力する。次に、入力した情報を数式１３に適用し、誤差楕円用の誤差共分散行列Ｐ~_{ｋ｜ｋ−１}を計算し出力する。

ステップＳ２５の誤差楕円の更新過程について説明する。この誤差楕円の更新過程は、ステップＳ２４の誤差楕円の予測過程で計算した誤差楕円用の誤差共分散行列Ｐ~_{ｋ｜ｋ−１}を入力し、さらに観測雑音Ｒ_ｋを入力する。入力した観測雑音Ｒ_ｋから位置ｐに関する成分のみを取り出す。次に、入力した情報及び取り出した情報を数式１４、１５に適用し、誤差楕円用の誤差共分散行列Ｐ~_ｋ｜ｋを計算し出力する。

ステップＳ２６の誤差楕円の軸方向及び径の長さの計算過程について説明する。ステップＳ２６の計算過程は、ステップＳ２５の誤差楕円の再更新過程で計算した誤差楕円用の誤差共分散行列Ｐ~_ｋ｜ｋを入力する。次に、入力した情報を数式１６〜１９に適用し、誤差楕円の軸方向及び径の長さを計算する。また、２次元測位演算が行われていれば、数式２０によって、誤差楕円の径の長さを補正する。計算された誤差楕円の軸方向及び径の長さは、衛星信号受信装置Ｒｘを使用するユーザ側に外部出力される。

本発明及び比較例の誤差楕円の時系列グラフを図１０に示す。図１０に示した時系列グラフは、衛星信号受信装置Ｒｘの搭載車両が図２のような高架下を走行時の、誤差楕円の時系列グラフである。「位置誤差」は、測位位置と真位置間の距離を示す。「従来誤差楕円」は、従来技術の誤差楕円の長半径を示す。「本発明誤差楕円」は、本発明の誤差楕円の長半径を示す。

誤差楕円は、以下の条件を満たすことが理想的である。（１）誤差楕円は位置誤差を上回ること。誤差楕円が位置誤差を下回れば、測位位置の誤差が大きいにも関わらず、ユーザは測位位置の誤差が実際より小さいと認識してしまう。（２）誤差楕円は位置誤差を過剰に上回らず少しだけ上回ること。誤差楕円が位置誤差を過剰に上回れば、測位位置の誤差が小さいにも関わらず、ユーザは測位位置の誤差が実際より大きすぎると認識してしまう。

従来の誤差楕円は、測位時刻２３：１７：３７以降において、位置誤差を下回っているため、上述の条件を満たしておらず理想的ではない。本発明の誤差楕円は、測位時刻２３：１７：３７以降を含めて、時系列の全測位時刻において、位置誤差を少しだけ上回っているため、上述の条件を満たしており理想的である。

以上に説明の事項は、以下のようにまとめられる。

本発明では、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルε_ｋが長くなり、差分ベクトルε_ｋが長いほど、差分ベクトルε_ｋを用いて算出したプロセス雑音Ｗ_ｋが大きくなり、プロセス雑音Ｗ_ｋを用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。

そして、マルチパスが測位位置に及ぼす影響が大きいほど、差分ベクトルε_ｋと速度ベクトルｖ^＾ _ｋ｜ｋのなすベクトル間角度θ_ｋが９０°又が２７０°に近くなり、ベクトル間角度θ_ｋが９０°又が２７０°に近いほど、ベクトル間角度θ_ｋを用いて算出した速度ベクトル誤差量σ_ｖが大きくなり、速度ベクトル誤差量σ_ｖを用いて算出したプロセス雑音Ｗ_ｋが大きくなり、プロセス雑音Ｗ_ｋを用いて算出した誤差楕円が大きくなるため、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。

さらに、測位演算に利用可能な測位衛星数が少なく、ＥＮＵ座標系の高さ方向（Ｕｐ軸）の測位位置を固定する２次元測位演算を行なうときに、ＥＮＵ座標系の水平面（Ｅ−Ｎ平面）の誤差楕円の長半径及び短半径を、ＥＮＵ座標系の高さ方向（Ｕｐ軸）における誤差楕円の誤差分散で補正することで、誤差楕円を大きくすることができる。逆に、測位演算に利用可能な測位衛星数が多く、ＥＮＵ座標系の高さ方向（Ｕｐ軸）の測位位置を固定しない３次元測位演算を行なうときは、ＥＮＵ座標系の水平面（Ｅ−Ｎ平面）の誤差楕円の長半径及び短半径を、ＥＮＵ座標系の高さ方向（Ｕｐ軸）における誤差楕円の誤差分散で補正しないため、誤差楕円を大きくし過ぎないようにすることができる。

本発明の衛星信号受信装置は、マルチパスが測位位置に及ぼす影響を測位位置の誤差楕円に対してよりよく反映させることができる。

Ｒｘ：衛星信号受信装置
Ｓ１、Ｓ２、Ｓ３、Ｓ４：測位衛星
Ｅｌ：高架構造
Ｒｆ：反射物体
１：信号受信部
２：追尾処理部
３：復調処理部
４：測位演算部
２１：擬似距離観測部
２２：ドップラー周波数観測部
２３：航法データ抽出部
４１：位置計算部
４２：速度計算部
４３：差分ベクトル計算部
４４：速度ベクトル誤差量計算部
４５：誤差楕円計算部

Claims

測位演算を行なう衛星信号受信装置であって、
自装置と測位衛星の間の擬似距離を観測する擬似距離観測部と、
前記擬似距離観測部が観測した擬似距離に基づいて、測位時刻より過去の時刻において、当該測位時刻における自装置の予測位置を予測した後に、当該測位時刻において、当該測位時刻における自装置の測位位置を更新する位置計算部と、
前記位置計算部が当該測位時刻より過去の時刻において予測した、当該測位時刻における自装置の予測位置と、前記位置計算部が当該測位時刻において更新した、当該測位時刻における自装置の測位位置と、の差分ベクトルを計算する差分ベクトル計算部と、
当該測位時刻における自装置の速度ベクトルを計算する速度計算部と、
前記差分ベクトル計算部が計算した差分ベクトルが長いほど、さらに当該測位時刻における自装置の速度ベクトル誤差量が大きいほど、当該測位時刻において、当該測位時刻における誤差楕円を大きく計算する誤差楕円計算部と、
を備えることを特徴とする衛星信号受信装置。
自装置と測位衛星の間のドップラー周波数を観測するドップラー周波数観測部と、
前記差分ベクトル計算部が計算した差分ベクトルと、前記速度計算部が計算した速度ベクトルと、のなすベクトル間角度を計算し、計算したベクトル間角度に基づいて、速度ベクトル誤差量を計算する速度ベクトル誤差量計算部と、
をさらに備え、
前記速度計算部は、前記ドップラー周波数観測部が観測したドップラー周波数に基づいて、当該測位時刻において、当該測位時刻における自装置の速度ベクトルを計算し、
前記誤差楕円計算部は、前記速度ベクトル誤差量計算部が計算したベクトル間角度が９０°又は２７０°に近いほど、当該測位時刻において、当該測位時刻における前記速度ベクトル誤差量を大きく計算することを特徴とする、請求項１に記載の衛星信号受信装置。
前記誤差楕円計算部は、ＥＮＵ（ＬｏｃａｌＥａｓｔ、Ｎｏｒｔｈ、Ｕｐ）座標系の高さ（Ｕｐ軸）方向における自装置の測位位置を固定する２次元測位演算が行なわれるときに、ＥＮＵ座標系の水平面（Ｅａｓｔ−Ｎｏｒｔｈ平面）での誤差楕円の長半径及び短半径を、前記差分ベクトル計算部が計算した差分ベクトルに基づいて推定したＥＮＵ座標系の高さ方向における誤差楕円の誤差分散で補正することを特徴とする、請求項１又は請求項２に記載の衛星信号受信装置。