(57)【要約】  観測データ選別部（１０３）は、ｎ機（ｎは３以上の整数）の測位衛星からのｎ個の観測データと前記ｎ個の観測データに対応させたｎ個の補正データとから得られるｎ個のＬ１波疑似距離の一重差量を観測量として最小二乗法測位を実行する。また、観測データ選別部（１０３）は、最小二乗法測位の実行により測位衛星ごとに得られたｎ個の残差の二乗和を計算し、ｎ個の残差を正規化する。更に、観測データ選別部（１０３）は、残差の二乗和とｎ個の正規化後の残差とを評価する。
2015145718.tif 000002
【特許請求の範囲】
【請求項１】
ｎ機（ｎは３以上の整数）の測位衛星からのｎ個の観測データと前記ｎ個の観測データに対応させたｎ個の補正データとから得られるｎ個のＬ１波疑似距離の一重差量を観測量として最小二乗法測位を実行する最小二乗法測位実行部と、
前記最小二乗法測位の実行により測位衛星ごとに得られたｎ個の残差の二乗和を計算する二乗和計算部と、
前記ｎ個の残差を正規化する正規化部と、
前記二乗和計算部により計算された二乗和と前記正規化部によるｎ個の正規化後の残差とを評価する評価部とを有することを特徴とする測位装置。
【請求項２】
前記評価部は、
前記二乗和を第１の閾値と比較し、
前記ｎ個の正規化後の残差の各々を第２の閾値と比較し、
前記二乗和が前記第１の閾値を超え、前記ｎ個の正規化後の残差のうちの少なくともいずれかが前記第２の閾値を超える場合に、前記ｎ個の観測データのうちのいずれかの観測データを破棄することを特徴とする請求項１に記載の測位装置。
【請求項３】
前記評価部は、
前記ｎ個の観測データのうち、正規化後の残差の絶対値が最大の測位衛星からの観測データを破棄することを特徴とする請求項２に記載の測位装置。
【請求項４】
観測したｍ機の測位衛星からのｍ個の観測データと前記ｍ個の観測データに対応させたｍ個の補正データとから得られるｍ個のＬ１波疑似距離の一重差量から求めた観測量のうち、ｎ機（３≦ｎ＜ｍ）となる組み合わせを選択し、組み合せごとに、最小二乗法測位を実行する最小二乗法測位実行部と、
前記最小二乗法測位の実行により測位衛星ごとに得られたｎ個の残差の二乗和を計算する二乗和計算部と、
前記ｎ個の残差を正規化する正規化部と、
前記二乗和計算部により計算された二乗和と前記正規化部によるｎ個の正規化後の残差とを評価する評価部とを有することを特徴とする測位装置。
【請求項５】
前記評価部は、
前記二乗和を第１の閾値と比較し、
前記ｎ個の正規化後の残差の各々を第２の閾値と比較し、
前記二乗和が前記第１の閾値を超え、前記ｎ個の正規化後の残差のうちの少なくともいずれかが前記第２の閾値を超える場合に、前記ｎ個の観測データのうちのいずれかの観測データを破棄することを特徴とする請求項４に記載の測位装置。
【請求項６】
前記評価部は、
前記ｎ個の観測データのうち、正規化後の残差の絶対値が最大の測位衛星からの観測データを破棄することを特徴とする請求項５に記載の測位装置。
【請求項７】
前記評価部は、
前記最小二乗法測位実行部により選択された組み合わせのうち、前記二乗和計算部により算出された二乗和が最小となる組み合わせのｎ個の観測データを、測位計算を行う測位計算部に入力することを特徴とする請求項４に記載の測位装置。
【請求項８】
前記評価部は、
測位衛星ごとに、測位信号の信号強度を第３の閾値と比較し、
測位信号の信号強度が前記第３の閾値を下回る測位衛星からの観測データを破棄することを特徴とする請求項１に記載の測位装置。
【請求項９】
前記評価部は、
測位衛星ごとに、衛星仰角を第４の閾値と比較し、
衛星仰角が前記第４の閾値を下回る測位衛星からの観測データを破棄することを特徴とする請求項１に記載の測位装置。
【請求項１０】
前記評価部は、
測位衛星ごとに、疑似距離二重差残差を第５の閾値と比較し、
疑似距離二重差残差が前記第５の閾値を超える測位衛星からの観測データを破棄することを特徴とする請求項１に記載の測位装置。
【請求項１１】
前記測位装置は、更に、
測位衛星ごとにアンビギュイティを計算するアンビギュイティ計算部を有し、
前記評価部は、
測位衛星ごとに、搬送波位相二重差残差を第６の閾値と比較し、
搬送波位相二重差残差が前記第６の閾値を超える測位衛星のアンビギュイティを再計算するよう前記アンビギュイティ計算部に指示し、
前記アンビギュイティ計算部は、
前記評価部からアンビギュイティの再計算を指示された測位衛星のアンビギュイティを再計算することを特徴とする請求項１に記載の測位装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、衛星測位を行う測位装置に関する。
【背景技術】
【０００２】
ＧＰＳ（Ｇｌｏｂａｌ  Ｐｏｓｉｔｉｏｎｉｎｇ  Ｓｙｓｔｅｍ）等の衛星測位システム（ＧＮＳＳ：Ｇｌｏｂａｌ  Ｎａｖｉｇａｔｉｏｎ  Ｓａｔｅｌｌｉｔｅ  Ｓｙｓｔｅｍ）からのＧＮＳＳ観測データには様々な誤差が含まれている。
高精度の測位を実現するためには、ＧＮＳＳ観測データに含まれる誤差を除去する必要がある。
例えば、特許文献１では、擬似距離と搬送波位相について誤差の二重差計算を行って、誤差を除去することが開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００９−２５７８０２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
測位精度を劣化させる原因の１つとしてマルチパスがある。
マルチパスは、観測データが含まれるＧＮＳＳ信号が複数の経路で測位装置のアンテナに到達することによって生じる現象である。
例えば、建物等に反射した後に測位装置のアンテナに到達するＧＮＳＳ信号（マルチパス信号）と、アンテナに直接到達するＧＮＳＳ信号とが干渉してマルチパスが生じる。
建物が多い都市部、森林が多い山間部では、マルチパスが発生しやすいため、これらの地域で高精度測位を達成するためには、マルチパス信号を取り除くことが必要である。
マルチパス信号は測位点の受信環境に強く依存するため、マルチパス信号を除去するための補正データを配信することはできない。
また、ＧＮＳＳ信号ごとにマルチパスの影響を受ける度合いが異なるため、特許文献１に記載の二重差計算でマルチパス信号を除去することもできない。【０００５】
本発明は、このような事情に鑑みたものであり、高精度測位を実現するために、マルチパスの影響がある観測データを排除することを主な目的とする。
【課題を解決するための手段】
【０００６】
本発明に係る測位装置は、
ｎ機（ｎは３以上の整数）の測位衛星からのｎ個の観測データと前記ｎ個の観測データに対応させたｎ個の補正データとから得られるｎ個のＬ１波疑似距離の一重差量を観測量として最小二乗法測位を実行する最小二乗法測位実行部と、
前記最小二乗法測位の実行により測位衛星ごとに得られたｎ個の残差の二乗和を計算する二乗和計算部と、
前記ｎ個の残差を正規化する正規化部と、
前記二乗和計算部により計算された二乗和と前記正規化部によるｎ個の正規化後の残差とを評価する評価部とを有することを特徴とする。
【発明の効果】
【０００７】
マルチパスの影響がある観測データ間では残差のばらつきが大きく、また、個々の残差の量もマルチパスの影響がある観測データでは多い。
本発明では、残差の二乗和を評価することにより残差のばらつきを評価し、また、正規化後の残差を評価することにより個々の残差の量を評価して、マルチパスの影響がある観測データを排除することができる。
【図面の簡単な説明】
【０００８】
【図１】実施の形態１に係る測位システムの構成例を示す図。
【図２】実施の形態１に係る測位装置の入出力を示す図。
【図３】実施の形態１に係る測位装置の入出力データを説明する図。
【図４】実施の形態１に係る観測データを説明する図。
【図５】実施の形態１に係るＧＰＳ信号に含まれる誤差要因を示す図。
【図６】実施の形態１に係るＧＰＳ信号の誤差と補強情報との関係を説明する図。
【図７】実施の形態１に係るグリッド点を示す図。
【図８】実施の形態１に係る測位装置の構成例を示す図。
【図９】実施の形態１に係る測位装置の構成要素を説明する図。
【図１０】実施の形態１に係る測位装置の中間データを説明する図。
【図１１】実施の形態１に係る観測データ選別部の動作例を示すフローチャート図。
【図１２】実施の形態１に係る観測データの内訳を説明する図。
【図１３】実施の形態１に係る二重差計算を説明する図。
【図１４】実施の形態１に係るカルマンフィルタの処理フローを示す図。
【図１５】実施の形態１に係るカルマンフィルタで用いられるベクトル及び行列を説明する図。
【図１６】実施の形態１に係る最小二乗法測位の手順を示すフローチャート図。
【図１７】実施の形態１に係る最小二乗法測位の計算式を示す図。
【図１８】実施の形態１に係る最小二乗法測位の計算式を示す図。
【図１９】実施の形態１及び実施の形態２に係る測位装置のハードウェア構成例を示す図。
【発明を実施するための形態】
【０００９】
実施の形態１．
１．システム構成
図１は、本実施の形態に係る測位システムの構成例を示す図である。
また、以下では、測位衛星としてＧＰＳ衛星を使用した例を示す。
ＧＰＳ衛星の代わりに、ＧＬＯＮＡＳＳ、Ｇａｌｉｌｅｏ、ＢｅｉＤｏｕ等のＧＮＳＳ、準天頂衛星等のＲＮＳＳ（Ｒｅｇｉｏｎａｌ  Ｎａｖｉｇａｔｉｏｎ  Ｓａｔｅｌｌｉｔｅ  Ｓｙｓｔｅｍ）を用いてもよい。
【００１０】
図１において、測位装置１００は、例えば自動車などの移動体に搭載される。
測位装置１００は、ＧＰＳ衛星３００から送信される観測データ５０１と放送暦５０２を受信する。
放送暦５０２は、エフェメリスともいう。
また、測位装置１００は、補強情報配信衛星２００から送信される補強情報４００を受信する。
本実施の形態及び実施の形態２では、主に測位装置１００の動作を説明する。
【００１１】
補強情報配信衛星２００は、図１に図示していない地上局から補強情報を受信し、受信した補強情報を補強情報４００として配信する。
【００１２】
ＧＰＳ衛星３００は、測位衛星であり、観測データ５０１及び放送暦５０２を送信する。
測位装置１００は、測位のためにＧＰＳ衛星３００を４機以上捕捉する必要がある。
【００１３】
２．測位装置１００の動作の概要
ここでは、本実施の形態及び実施の形態２で説明する測位装置１００の動作の概要を説明する。
測位装置１００は、図２に示すように、ＧＰＳ信号として、観測データ５０１と放送暦５０２を受信し、補強情報４００を受信する。
そして、測位装置１００は補強情報４００と観測データ５０１及び放送暦５０２を用いて測位点の位置（測位装置１００の位置）を算出する。
観測データ５０１、放送暦５０２、補強情報４００、位置の詳細は図３に示す通りである。
【００１４】
３．１．観測データ
観測データ５０１からは、測位点とＧＰＳ衛星３００との間の疑似距離と、搬送波位相を導出することができる。
疑似距離と搬送波位相の詳細は、図４に示す通りである。
観測データ５０１から導出される疑似距離及び搬送波位相には、それぞれ誤差が含まれる。
測位装置１００は、補強情報４００を用いて、疑似距離及び搬送波位相に含まれる誤差を除去する。
なお、以降では、ＧＰＳ衛星ｉのＬ１波疑似距離をＰ（ｉ，１）と表記し、ＧＰＳ衛星ｉのＬ２波疑似距離をＰ（ｉ，２）と表記する。
また、ＧＰＳ衛星ｉのＬ１波搬送波位相をΦ（ｉ，１）と表記し、ＧＰＳ衛星ｉのＬ２波搬送波位相をΦ（ｉ，２）と表記する。
【００１５】
３．２．補強情報
観測データ５０１に含まれるバイアス誤差及び雑音要因を図５に示す。
ＧＰＳ衛星３００に起因する誤差として軌道誤差、衛星時計誤差、周波数間バイアスがあり、信号の伝搬経路に起因する誤差として電離層伝搬遅延誤差（電離層遅延誤差又は電離層遅延量ともいう）及び対流圏伝搬遅延誤差（対流圏遅延誤差又は対流圏遅延量ともいう）がある。
さらに、測位装置１００の受信機に起因する誤差として、受信機時計誤差、受信機雑音、さらに建物に反射したＧＰＳ信号とＧＰＳ衛星３００から直接受信したＧＰＳ信号が干渉して生じるマルチパスがある。
これらの誤差のうち、受信機に起因する誤差は、ユーザが使用する測位装置１００の受信機の性能、また、受信環境によって異なるため補正データ及び補強情報には含まれず、測位装置１００の処理において受信機に起因する誤差が取り除かれる。
ＧＰＳ衛星３００に起因する誤差、信号伝搬経路に起因する誤差は圧縮されて補強情報として補強情報配信衛星２００から配信される。
図５に示す誤差以外に、測位点位置によって異なる地球潮汐効果（Ｅａｒｔｈ  Ｔｉｄｅ効果）による誤差、Ｐｈａｓｅ  Ｗｉｎｄ  Ｕｐ効果による誤差は補正データには含まれるが、補強情報には含まれない。
以上を考慮した補強情報の内訳を図６に示す。
【００１６】
観測データ５０１は電子基準点から地上回線を介して地上局に収集され、地上局で補正データが生成される。
従来では、補正データは、補正データの信頼性情報と合わせて、補強情報として、携帯電話あるいは無線ＬＡＮ（Ｌｏｃａｌ  Ａｒｅａ  Ｎｅｔｗｏｒｋ）経由で更新周期１秒で配信されているが、本実施の形態では、伝送容量の限られた人工衛星の通信回線を介してユーザに補強情報を放送するため、本実施の形態に係る補強情報は時間変動の激しさに応じて誤差を高レート、低レートに分類して時間圧縮を行っている。
より具体的には、高レートの誤差は５秒ごとに更新され、低レートの誤差は３０秒ごとに更新される。
また、従来の補正データは約１０〜３０ｋｍ間隔に設定した電子基準点ごとに生成され、配信されているが、本実施の形態に係る補強情報は信号経路に起因する（空間変動の有る）電離層遅延誤差と対流圏遅延誤差について６０ｋｍ四方のグリッド点（図７）ごとのみ配信するようにして空間圧縮した。
さらに、本実施の形態では、補強情報に含まれる誤差を周波数に依存する誤差（周波数依存誤差）と周波数に依存しない誤差（周波数非依存誤差）に分類している。
周波数非依存誤差は高レートの誤差に分類され、周波数依存誤差は低レートの誤差に分類される。
周波数非依存誤差のうち、衛星時計誤差のみ５秒ごとに更新され、配信される。
周波数非依存誤差の他の誤差、すなわち、衛星軌道誤差、周波数間バイアス、対流圏遅延誤差は３０秒ごとに更新され、配信される。
しかしながら、衛星軌道誤差、周波数間バイアス、対流圏遅延誤差については、５秒ごとに、３０秒ごとの誤差からの変化量が補正値として衛星時計誤差（５秒更新、配信）に付加される。
つまり、衛星軌道誤差、周波数間バイアス、対流圏遅延誤差については、３０秒の間に、３０秒／５秒−１＝５個の補正値が５秒ごとの衛星時計誤差に付加される。
このため、測位装置１００は、３０秒ごとに衛星軌道誤差、周波数間バイアス、対流圏遅延誤差の更新値を受信するとともに、５秒ごとに衛星軌道誤差、周波数間バイアス、対流圏遅延誤差の補正値を受信することができる。
そして、測位装置１００は、３０秒ごとの更新値に５秒ごとの補正値を加算することで、５秒ごとに、衛星軌道誤差、周波数間バイアス、対流圏遅延誤差の補正を行うことができる。
衛星時計誤差に付加される衛星軌道誤差、周波数間バイアス、対流圏遅延誤差の５秒ごとの補正値をコンシステンシーともいう。
このように、本実施の形態では、５秒ごとにコンシステンシーを配信することで、補強情報のデータ量を圧縮している。
なお、周波数間バイアスは、Ｌ１信号を基準とし、各信号の衛星搭載ペイロード内の遅延量を表したものであり、周波数間バイアス（Ｌ１バイアス）は、‘０’、であり、周波数間バイアス（Ｌ２バイアス）は、Ｌ１搬送波信号に対するＬ２搬送波信号の遅延量‘Ｌ２−Ｌ１’を表す。
また、周波数間バイアス（Ｌ０バイアス）は、周波数間バイアスの周波数非依存項を示し、図６の式１に示されるように、周波数間バイアス（Ｌ１バイアス）と周波数間バイアス（Ｌ２バイアス）から算出される。
周波数間バイアスはＬ１信号を基準としたが他の信号を基準にしても問題はなく、信号もＬ５を含めても同じ効果を期待できる。
【００１７】
４．測位装置の構成例
図８は、本実施の形態に係る測位装置１００の構成例を示す。
また、図８に示す各構成要素の概略説明を図９に示し、中間データの概略説明を図１０に示す。
【００１８】
概略位置及び衛星位置計算部１０１は、ＧＰＳ衛星３００から、観測データ５０１及び放送暦５０２を受信し、測位点の概略位置と各ＧＰＳ衛星３００の位置を算出する。
概略位置及び衛星位置計算部１０１の算出結果が、概略位置１５１と衛星位置１５２である。
概略位置１５１は、単独測位によって算出されたメートルオーダーの精度の測位点位置である。
衛星位置１５２は、測位装置１００が観測データを受信した各ＧＰＳ衛星３００の位置である。
【００１９】
補正データ作成部１０２は、補強情報配信衛星２００から補強情報４００を受信し、また、概略位置１５１、衛星位置１５２を取得し、補強情報４００、概略位置１５１、衛星位置１５２から補正データ１５３を算出する。
補正データ１５３には、測位点で各ＧＰＳ衛星３００から受信した観測データ５０１に含まれると予想される誤差が示される。
【００２０】
観測データ選別部１０３は、品質が劣化していると予想される観測データ５０１を除去する。
観測データ選別部１０３は、最小二乗法測位実行部、二乗和計算部、正規化部及び評価部の例に相当する。
【００２１】
観測データ誤差補正部１０４は、二重差計算を行って、観測データの二重差データ１５４を出力する。
二重差データ１５４には、従衛星の観測データ（補正データ１５３を使って補正済みの観測データ）から主衛星の観測データ（補正データ１５３を使って補正済みの観測データ）を差し引いた量が示される。
二重差計算及び二重差データ１５４については後述する。
【００２２】
時間外挿計算部１０５は、時間外挿計算を行って、前のエポックの状態量Ｘ＾（ｔ−Δｔ）から現在エポックの状態量Ｘ（ｔ）を推定する。
なお、図８に示している「Ｘ」の真上に「＾」が位置している表記と、「Ｘ」の右上に「＾」が位置している表記（「Ｘ＾」）は同じ意味である。
また、「＾」は、後述する観測更新計算部１０８により更新された後の状態量であることを意味する。
【００２３】
幾何学距離計算部１０６は、衛星位置１５２から、ＧＰＳ衛星３００から測位点までの幾何学距離１５５を計算する。
【００２４】
残差計算部１０７は、二重差データ１５４と幾何学距離１５５から二重差残差１５６を計算する。
【００２５】
観測更新計算部１０８は、状態量Ｘ（ｔ）の推定誤差が最も小さくなるように状態量Ｘ（ｔ）を更新する。
観測更新計算部１０８により更新された後の状態量Ｘ（ｔ）を状態量Ｘ＾（ｔ）と表記する。
【００２６】
アンビギュイティ計算部１０９は、搬送波位相のバイアス量であるアンビギュイティを算出し、算出結果をもとに状態量Ｘ＾（ｔ）を更新する。
アンビギュイティ計算部１０９により更新された状態量Ｘ＾（ｔ）に含まれる位置の値が、測位結果として出力される。
また、アンビギュイティ計算部１０９により更新された状態量Ｘ＾（ｔ）は、１つ前のエポックの状態量Ｘ＾（ｔ−Δｔ）として、時間外挿計算部１０５により時間外挿計算の対象となる。
【００２７】
二重差残差解析部１１０は、マルチパスの疑いのある観測データを特定し、マルチパスの疑いのある観測データを観測データ選別部１０３に通知する。
また、二重差残差解析部１１０は、サイクルスリップの疑いがある場合に、アンビギュイティの再計算をアンビギュイティ計算部１０９に指示する。
二重差残差解析部１１０は、観測データ選別部１０３とともに評価部の例に相当する。
【００２８】
４．１．観測データ選別部１０３
ここでは、観測データ選別部１０３の動作の詳細を説明する。
【００２９】
観測データ選別部１０３は、一重差を用いた残差検出及び除去方式により、複数の観測データの中からマルチパスの影響のある観測データを検出し、検出した観測データを取り除く。
図１１は、観測データ選別部１０３の動作例を示す。
以下では、図１１を参照して観測データ選別部１０３の動作例を説明する。
【００３０】
観測データ選別部１０３は、ＧＰＳ衛星３００が３機以上捕捉できている場合（Ｓ１１０１でＹＥＳ）に、後述の４．２．１項で説明するＬ１波疑似距離の一重差量（以下、単に一重差ともいう）を観測量とする最小二乗法測位を実行する（Ｓ１１０２）。
つまり、観測データ選別部１０３は、ｎ機（ｎは３以上の整数）のＧＰＳ衛星３００からのｎ個の観測データ５０１とｎ個の観測データ５０１に対応させたｎ個の補正データ１５３とから得られるｎ個のＬ１波疑似距離の一重差量を観測量として最小二乗法測位を実行する。
最小二乗法測位の処理手順は図１６に示す通りであり、また、最小二乗法測位の計算式は図１７及び図１８に示す通りである。
図１６の最小二乗法測位の処理手順、図１７及び図１８の最小二乗法測位の計算式は、それぞれ既知のものであるため、説明を省略する。
最小二乗法測位の実行によりＧＰＳ衛星３００ごとに観測量の残差が得られる。
つまり、ｎ機のＧＰＳ衛星３００に対してｎ個の残差が得られる。
【００３１】
最小二乗法測位により得られたｎ個の残差間のばらつきの大きさは、観測データにマルチパスが含まれていない場合に比べて、観測データにマルチパスが含まれている場合の方が大きい。
また、ｎ個の残差の個々の大きさも、観測データにマルチパスが含まれていない場合に比べて、観測データにマルチパスが含まれている場合の方が大きい。
観測データ選別部１０３は、各衛星の観測量の残差の二乗和を計算する（ｎ個の残差の二乗和を計算する）ことで、ｎ個の残差間の大きさのばらつきを評価する（Ｓ１１０３）。
また、観測データ選別部１０３は、ｎ個の残差を正規化する（Ｓ１１０３）。
二乗和の計算式及び正規化の計算式は以下の通りである。
なお、以下において、Ｒは観測ノイズであり、Ｈは観測行列であり、ｎは捕捉した衛星数であり、添字ｉ及びｊは、衛星番号である。
また、行列Ｍは残差の予測される分散値を示す行列である。
Ａの予測値をＥ＜Ａ＞で表わすとするとｉ行ｊ列のＭの要素Ｍ（ｉ，ｊ）は、
Ｍ（ｉ，ｊ）＝Ｅ＜ｄｚ（ｉ）＊ｄｚ（ｊ）＞
で表わされる。
【００３２】
【数１】
2015145718.tif 000003
【００３３】
次に、観測データ選別部１０３は、残差の二乗和ｚｚを閾値（第１の閾値）と比較する（Ｓ１１０４、Ｓ１１０５）。
残差の二乗和ｚｚが閾値を超えない場合（Ｓ１１０５でＮＯ）は、観測データ選別部１０３は、ｎ個の観測データのいずれにもマルチパスが含まれていないと判断して、処理を終了する。
一方、残差の二乗和ｚｚが閾値を超える場合（Ｓ１１０５でＹＥＳ）は、観測データ選別部１０３は、次に、ｎ個の正規化後の残差ｚｎｉのそれぞれを閾値（第２の閾値）と比較する（Ｓ１１０４、Ｓ１１０６）。
ｎ個の正規化後の残差ｚｎｉの全てが閾値を超えない場合（Ｓ１１０６でＮＯ）は、観測データ選別部１０３は、ｎ個の観測データのいずれにもマルチパスが含まれていないと判断して、処理を終了する。
一方、ｎ個の正規化後の残差ｚｎｉのうちの１つでも閾値を超える場合（Ｓ１１０６でＹＥＳ）は、観測データ選別部１０３は、ｎ個の正規化後の残差ｚｎｉのうち絶対値が最大の残差ｚｎｉに対応するＧＰＳ衛星３００からの観測データをマルチパスを含むとして破棄する（Ｓ１１０７）。
【００３４】
また、観測データ選別部１０３は、以下の１）〜３）のうちのいずれかの条件を満たす観測データはマルチパスの疑いがあると判断して取り除く。
１）信号強度が設定した閾値（第３の閾値）より低い
観測データにはＧＰＳ信号の信号強度も含まれる。
マルチパスが含まれる観測データは信号強度が下がる傾向にあるため、観測データ選別部１０３は、信号強度が閾値より低い観測データは取り除く。
２）衛星仰角が設定した閾値（マスク角）（第４の閾値）より小さい
仰角が低いＧＰＳ衛星からの観測データは通過する電離層及び対流圏の距離も長く、建物に反射する可能性が高いため、観測データ選別部１０３は、仰角が低いＧＰＳ衛星からの観測データを取り除く。
３）擬似距離二重差残差（４．２．２項参照）がカルマンフィルタによる推定誤差（第５の誤差）を超える
二重差残差解析部１１０が、カルマンフィルタが二重差観測量に含まれると推定する推定誤差（観測値誤差）と擬似距離二重差残差とを比較し、擬似距離二重差残差が推定誤差を超える場合に、観測データ選別部１０３は該当する観測データを取り除く。
カルマンフィルタによる推定誤差は以下の式で表わされる行列σの対角成分σｉｉで計算される。
σ＝ＨＰＨ＋Ｒ
ＨとＲは観測更新計算部１０８で観測更新に使用する観測行列（図１４）と観測ノイズ（図１４）を用い、Ｐは時間外挿計算部１０５で推定する誤差共分散行列（図１４）を用いる。
二重差残差解析部１１０が、衛星番号ｉのＧＰＳ衛星の疑似距離二重差の観測値誤差と疑似距離二重差の残差を比較し、残差の方が大きい場合は、衛星番号ｉを観測データ選別部１０３に通知し、観測データ選別部１０３が、衛星番号ｉのＧＰＳ衛星からの観測データを取り除く。
【００３５】
このようにして観測データ選別部１０３によりマルチパスが含まれている観測データが取り除かれ、マルチパスが含まれていない観測データのみが観測データ誤差補正部１０４に入力される。
また、観測データ選別部１０３は、観測された衛星の内、一部の複数の衛星を選択し、選択する衛星を変え、上記一重差を用いた残差検出及び除去方式による処理を実施し、マルチパスの検出を実施し、最も残差が少ない衛星の組み合わせの観測データが誤差補正部１０４に入力される。
この衛星選択において、選択する衛星数は、可視衛星数を最大とし、測位計算可能な４機を最小とする。
この選択は可視衛星のすべての組み合わせで行う。
【００３６】
４．２．１．二重差量
次に、観測データ誤差補正部１０４による二重差計算について説明する。
図１２及び図１３は、二重差計算を説明する。
【００３７】
図４に示したように観測データには疑似距離と搬送波位相の２種類がある。
疑似距離と搬送波位相は、それぞれ、前述の３．２項及び図６で示した誤差を含んでおり、誤差には補正データに含まれる量（図１２の符号１２０１）と、補正データに含まれない受信機に起因する受信機時計誤差（図１２の符号１２０３）及び受信機雑音（図１２の符号１２０４）がある（ここではマルチパスは無視する）。
誤差を全て取り除いた部分が測位点とＧＰＳ衛星３００の間の真の距離（幾何学距離、図１２の符号１２０２）になる。
各ＧＰＳ衛星３００の観測データ（疑似距離と搬送波位相）から補正データを差し引くことで補正データに含まれる誤差量（図１２の符号１２０１）が取り除かれる。
観測データ（疑似距離と搬送波位相）から補正データを差し引いた後の誤差量を一重差量と呼ぶ。
一つのＧＰＳ衛星を主衛星と決め（一般的に天頂付近にあるＧＰＳ衛星が選択される）、主衛星以外のＧＰＳ衛星を従衛星とし、従衛星の一重差から主衛星の一重差を差し引く。
図１３に示しているように、従衛星の一重差と主衛星の一重差は共通の受信機の一重差であるため、従衛星の一重差に含まれる受信機時計誤差（図１２の符号１２０３）と主衛星の一重差に含まれる受信機時計誤差は等しい。
このため、従衛星の一重差から主衛星の一重差を差し引くことで、受信機時計誤差（図１２の符号１２０３）はキャンセルされる。
受信機時計誤差が除去された後の誤差量を二重差量と呼ぶ（単に二重差ともいう）。
二重差量を計算することで主な誤差は取り除かれ、幾何学距離（図１２の符号１２０２）と受信機雑音（図１２の符号１２０４）とアンビギュイティ（図１２の符号１２０５、搬送波位相のみ）が残る。
受信機雑音は時間で平均すると０になる雑音であるため、観測を継続する過程で統計処理（後述の４．２．２項のカルマンフィルタ）を行うことで取り除くことができる。
アンビギュイティは、誤差が取り除かれた疑似距離を使用して推定することができ（後述の４．２．３項参照）、この結果、高精度な位置の推定が可能になる。
【００３８】
４．２．２．カルマンフィルタ
図８と図９に示した時間外挿計算部１０５及び観測更新計算部１０８を実現するカルマンフィルタについて説明する。
図１４は、カルマンフィルタの処理フローを示す。
また、図１５は、カルマンフィルタの処理に用いられる変数の説明を示す。
【００３９】
図８の時間外挿計算部１０５では、図１４に示すカルマンフィルタの時間外挿計算が行われる。
また、図８の観測更新計算部１０８では、図１４に示すカルマンフィルタの観測更新計算が行われる。
時間外挿計算と観測更新計算で１つのループを構成し、時間外挿計算と観測更新計算で構成されるループが繰り返し実行される。
カルマンフィルタは、繰り返される各ループにおいて、推定する状態量（状態量Ｘ）の誤差分散（誤差分散行列Ｐｉｊ＝Ｅ〈ｘｉｘｊ〉、Ｅ〈ａ〉はａの分散）の対角成分が最小になるように状態量を推定する。
【００４０】
カルマンフィルタにおける処理を順を追って説明する。
時間外挿計算では、前の時刻の状態量（ｘ＾（−））と誤差分散行列（Ｐ＾（−））から、採用する運動モデルに従って決定される遷移行列Φによって次の時刻の状態量（ｘ（＋））と誤差共分散行列（Ｐ（＋））が推定される。
この時、運動モデルと実際の現象との間の予想される誤差であるプロセスノイズＱを誤差共分散行列（Ｐ＾（−））に加算する。
このプロセスノイズＱも採用する運動モデルと設計によって決定される。
推定した状態量（ｘ（＋））から、状態量と観測量の関係を表わす観測モデル（ｙ―＝ｆ（ｘ））により推定される観測量相当量ｙ―を求める（ｙ―は、「ｙ」の真上に「−」が位置していることを表す、以下も同様）。
観測更新計算では、実際の観測量との差である残差（ｄｚ＝ｙ−ｙ―）を求め、図１４の式で表わされるカルマンゲインＫを使用して状態量の差に変換（ｄｘ＝Ｋ・ｄｚ）し、状態量を更新する。
観測更新計算で使用する観測行列は観測モデルを表現した行列であり、次式で求められる。
【００４１】
【数２】
2015145718.tif 000004
【００４２】
また、カルマンゲインＫの式で分母に含まれるＲは観測量に含まれると予想される観測ノイズを示す。
【００４３】
４．２．３．アンビギュイティ計算
次に、図８及び図９に示すアンビギュイティ計算部１０９で行われるアンビギュイティ計算について説明する。
従来においても、高精度測位を実現するためには受信機雑音がミリメートルオーダーである搬送波位相を使用して測位を行う必要がある（図４参照のこと）。
搬送波位相にはバイアス量であるアンビギュイティが含まれるが、二重差計算（４．２．１項参照）で取り除くことができないので、状態量（位置、速度）にアンビギュイティを加えて測位装置１００でアンビギュイティを推定して取り除く必要がある。
アンビギュイティを取り除く手順を以下に箇条書きで記す。
（１）疑似距離を使用して位置を推定する（４．２．１項、４．２．２項参照のこと）。
（２）上記（１）で求めた位置から幾何学距離が分かるので搬送波位相からアンビギュイティを推定する（４．２．１項参照のこと）。
（３）上記（１）と（２）をしばらく繰り返し統計的に状態量（位置、速度、アンビギュイティ）の分散が小さくなるようにする。
（４）上記（３）の段階では小数値であるアンビギュイティ（Ｆｌｏａｔ解と呼ぶ）をＬＡＭＢＤＡ法と呼ばれる手法を使って整数値（アンビギュイティは整数値、図４参照のこと）に確定する（Ｆｉｘ解と呼ぶ）。
（５）Ｆｉｘ解が正しいか検定を行い、正しければアンビギュイティを固定してＦｌｏａｔ解とＦｉｘ解との差分の位置を更新してセンチメートル精度を実現する。
【００４４】
ＬＡＭＢＤＡ法はアンビギュイティのＦｌｏａｔ解とＦｌｏａｔ解に対応した成分の誤差共分散行列を入力として、Ｆｌｏａｔ解を中心とする誤差共分散の大きさで決定される超次元楕円体球の中から、Ｆｌｏａｔ解とＦｉｘ解との差の二乗和平均が最も小さくなるＦｉｘ解を探索する手法の一つである。
超次元であるのは二重差の個数＝衛星数−１だけＦｌｏａｔ解があるためで、球ではなく楕円体球であるのは誤差共分散行列の非対角成分が０ではないためである。
上記（５）の検定はＦｉｘ解を用いた残差の大きさなどを判断材料にして行う。
【００４５】
実施の形態２．
５．サイクルスリップ対策
本実施の形態では、図８に示した二重差残差解析部１１０によるサイクルスリップ対策を説明する。
【００４６】
測位装置１００がＧＰＳ信号を受信し続けている間（ＧＰＳ衛星を捕捉している間）は、搬送波位相のアンビギュイティは一定値を保っている。
しかし、捕捉が途切れた後に、測位装置１００がＧＰＳ衛星３００を再捕捉した時にアンビギュイティが別の値にオフセットする場合がある（これをサイクルスリップと呼ぶ）。
サイクルスリップが発生すると、再捕捉前に測位装置１００が保持していたアンビギュイティを使用するとオフセットによる齟齬によって測位精度が劣化する。
従って、マルチパスと同様に、特に都市部、山間部など頻繁にＧＰＳ衛星の捕捉が途切れる地域では、高精度測位を達成するためにサイクルスリップ対策が必要である。
【００４７】
二重差残差解析部１１０は、ＧＰＳ衛星３００ごとに、搬送波位相二重差残差（４．２．２項参照）を閾値（第６の閾値）と比較する。
そして、二重差残差解析部１１０は、搬送波位相二重差残差が閾値を超える場合は、サイクルスリップの疑いがあると判断し、搬送波位相二重差残差が閾値を超えるＧＰＳ衛星３００のアンビギュイティを再計算するようアンビギュイティ計算部１０９に指示する。
アンビギュイティ計算部１０９は、二重差残差解析部１１０からアンビギュイティの計算を指示されたＧＰＳ衛星のアンビギュイティを再計算する。
【００４８】
６．変形例
以上、本発明の実施の形態について説明したが、これらの実施の形態のうち、２つ以上を組み合わせて実施しても構わない。
あるいは、これらの実施の形態のうち、１つを部分的に実施しても構わない。
あるいは、これらの実施の形態のうち、２つ以上を部分的に組み合わせて実施しても構わない。
なお、本発明は、これらの実施の形態に限定されるものではなく、必要に応じて種々の変更が可能である。
【００４９】
７．ハードウェア構成例
最後に、実施の形態１及び実施の形態２に示した測位装置１００のハードウェア構成例を図１９を参照して説明する。
測位装置１００はコンピュータであり、測位装置１００の各要素をプログラムで実現することができる。
測位装置１００のハードウェア構成としては、バスに、演算装置９０１、外部記憶装置９０２、主記憶装置９０３、通信装置９０４、入出力装置９０５が接続されている。
【００５０】
演算装置９０１は、プログラムを実行するＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）である。
外部記憶装置９０２は、例えばＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やフラッシュメモリ、ハードディスク装置である。
主記憶装置９０３は、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）である。
通信装置９０４は、ＧＰＳ衛星からの観測データ及び放送暦を受信し、また、補強情報配信衛星からの補強情報を受信する。
通信装置９０４には、ＡＤ（アナログ−ディジタル）変換機能が含まれている。
入出力装置９０５は、例えばタッチパネル式のディスプレイ装置等である。
【００５１】
プログラムは、通常は外部記憶装置９０２に記憶されており、主記憶装置９０３にロードされた状態で、順次演算装置９０１に読み込まれ、実行される。
プログラムは、図８に示す「〜部」として説明している機能を実現するプログラムである。
更に、外部記憶装置９０２にはオペレーティングシステム（ＯＳ）も記憶されており、ＯＳの少なくとも一部が主記憶装置９０３にロードされ、演算装置９０１はＯＳを実行しながら、図８に示す「〜部」の機能を実現するプログラムを実行する。
また、実施の形態１及び実施の形態２の説明において、「〜の補正」、「〜の生成」、「〜の作成」、「〜の計算」、「〜の算出」、「〜の判断」、「〜の判定」、「〜の更新」、「〜の推定」、「〜の抽出」、「〜の選択」、「〜の受信」等として説明している処理の結果を示す情報やデータや信号値や変数値が主記憶装置９０３にファイルとして記憶されている。
【００５２】
なお、図１９の構成は、あくまでも測位装置１００のハードウェア構成の一例を示すものであり、測位装置１００のハードウェア構成は図１９に記載の構成に限らず、他の構成であってもよい。
【符号の説明】
【００５３】
１００  測位装置、１０１  概略位置及び衛星位置計算部、１０２  補正データ作成部、１０３  観測データ選別部、１０４  観測データ誤差補正部、１０５  時間外挿計算部、１０６  幾何学距離計算部、１０７  残差計算部、１０８  観測更新計算部、１０９  アンビギュイティ計算部、１１０  二重差残差解析部、２００  補強情報配信衛星、３００  ＧＰＳ衛星、４００  補強情報、５０１  観測データ、５０２  放送暦。
【図１】
2015145718.tif 000005
【図２】
2015145718.tif 000006
【図３】
2015145718.tif 000007
【図４】
2015145718.tif 000008
【図５】
2015145718.tif 000009
【図６】
2015145718.tif 000010
【図７】
2015145718.tif 000011
【図８】
2015145718.tif 000012
【図９】
2015145718.tif 000013
【図１０】
2015145718.tif 000014
【図１１】
2015145718.tif 000015
【図１２】
2015145718.tif 000016
【図１３】
2015145718.tif 000017
【図１４】
2015145718.tif 000018
【図１５】
2015145718.tif 000019
【図１６】
2015145718.tif 000020
【図１７】
2015145718.tif 000021
【図１８】
2015145718.tif 000022
【図１９】
2015145718.tif 000023
【手続補正書】
【提出日】平成28年7月12日(2016.7.12)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】００３９
【補正方法】変更
【補正の内容】
【００３９】
図８の時間外挿計算部１０５では、図１４に示すカルマンフィルタの時間外挿計算が行われる。
また、図８の観測更新計算部１０８では、図１４に示すカルマンフィルタの観測更新計算が行われる。
時間外挿計算と観測更新計算で１つのループを構成し、時間外挿計算と観測更新計算で構成されるループが繰り返し実行される。
カルマンフィルタは、繰り返される各ループにおいて、推定する状態量（状態量Ｘ）の誤差共分散（誤差共分散行列Ｐｉｊ＝Ｅ〈ｘｉｘｊ〉、Ｅ〈ａ〉はａの分散）の対角成分が最小になるように状態量を推定する。
【手続補正２】
【補正対象書類名】明細書
【補正対象項目名】００４０
【補正方法】変更
【補正の内容】
【００４０】
カルマンフィルタにおける処理を順を追って説明する。
時間外挿計算では、前の時刻の状態量（ｘ＾（−））と誤差共分散行列（Ｐ＾（−））から、採用する運動モデルに従って決定される遷移行列Φによって次の時刻の状態量（ｘ（＋））と誤差共分散行列（Ｐ（＋））が推定される。
この時、運動モデルと実際の現象との間の予想される誤差であるプロセスノイズＱを誤差共分散行列（Ｐ＾（−））に加算する。
このプロセスノイズＱも採用する運動モデルと設計によって決定される。
推定した状態量（ｘ（＋））から、状態量と観測量の関係を表わす観測モデル（ｙ―＝ｆ（ｘ））により推定される観測量相当量ｙ―を求める（ｙ―は、「ｙ」の真上に「−」が位置していることを表す、以下も同様）。
観測更新計算では、実際の観測量との差である残差（ｄｚ＝ｙ−ｙ―）を求め、図１４の式で表わされるカルマンゲインＫを使用して状態量の差に変換（ｄｘ＝Ｋ・ｄｚ）し、状態量を更新する。
観測更新計算で使用する観測行列は観測モデルを表現した行列であり、次式で求められる。
【手続補正３】
【補正対象書類名】特許請求の範囲
【補正対象項目名】請求項５
【補正方法】変更
【補正の内容】
【請求項５】
前記評価部は、
前記二乗和を第１の閾値と比較し、
前記ｎ個の正規化後の残差の各々を第２の閾値と比較し、
前記二乗和が前記第１の閾値を超え、前記ｎ個の正規化後の残差のうちの少なくともいずれかが前記第２の閾値を超える場合に、ｎ個の観測データのうちのいずれかの観測データを破棄することを特徴とする請求項４に記載の測位装置。
【手続補正４】
【補正対象書類名】図面
【補正対象項目名】図１１
【補正方法】変更
【補正の内容】
【図１１】
2015145718.tif 000024
【国際調査報告】
2015145718.tif 000025
2015145718.tif 000026
2015145718.tif 000027
2015145718.tif 000028