JP5305416B2 - 衛星航法システムにおける電離圏異常を検出する方法及びその装置。 - Google Patents

Description

この発明は、衛星航法システムにおける電離圏異常を検出する方法及びその装置に関し、特に、１周波数の観測データのみで電離圏遅延を推定し、電離圏異常を検出する方法及びその装置に関するものである。

一般に、人工衛星や天体から発せられる電磁波信号が、地上に到達するまでの間に、電離圏及び対流圏を通過するが、それぞれの領域を電磁波信号が通過する際に遅延が生じる。これらの遅延は、それぞれ電離圏遅延及び対流圏遅延と呼ばれている。従って、この電磁波信号を測位信号として用いる場合、これらの電離圏遅延及び対流圏遅延が測位誤差源のひとつとなっている。

測位誤差の原因のひとつである電離圏遅延量の補正は勿論のこと、その他の補正情報を提供するために、世界的に実用化が進められているシステムとしては、航法衛星を使って補正情報を送信する地上型航法衛星補強システム（Ｇｒｏｕｎｄ−Ｂａｓｅｄ Ａｕｇｍｅｎｔａｔｉｏｎ Ｓｙｓｔｅｍ：以下、ＧＢＡＳと記す。）がある。

このＧＢＡＳは、航空機の着陸のために地上局によって補強されたＧＮＳＳ（Ｇｌｏｂａｌ Ｎａｖｉｇａｔｉｏｎ Ｓａｔｅｌｌｉｔｅ Ｓｙｓｔｅｍ）を用いた衛星航法システムで、空港周辺の約２０ＮＭ（約３７ｋｍ）以内の空域において、衛星信号（ＧＮＳＳ信号）を地上から補強することにより、航空機の進入着陸誘導に利用しようとする次世代の航空機着陸誘導システムである。なお、静止衛星を使った補正情報を送信する衛星航法補強システム（Ｓａｔｅｌｌｉｔｅ Ｂａｓｅｄ Ａｕｇｍｅｎｔａｔｉｏｎ Ｓｙｓｔｅｍ：以下、ＳＢＡＳと記す。）があり、このＳＢＡＳは日本が受け持つ管制空域内全域を飛行する航空機を誘導するために使用されるシステムである。

従って、ＧＢＡＳでは、非常に高い精度の補強情報が要求される。そのため、基準局とモニタ局間の基線長が長い場合（長基線）でも短い場合（短基線）でも、航法衛星から受信した衛星信号の電離圏遅延や対流圏遅延等の影響を除去しようとの要望がある。

そこで、従来から電離圏遅延や対流圏遅延等の影響を除去するための各種の方法があるが、その１つの方法として、測位用衛星から送信されるＬ１波とＬ２波との２つの周波数の衛星信号を用いて、基準周波数信号とする整数値バイアスを、基線長の長短に係わらず高精度に決定することが出来る相対測位装置がある。

これは、下記特許文献１に記載のように、測位用衛星から送信されるそれぞれ異なる２つの周波数の衛星信号を受信し、推定手段で、１つの基準周波数の衛星信号に対して異なる周波数の衛星信号を差分合成してなるワイドレーン結合の整数値バイアスを推定演算するとともに、基準周波数の衛星信号の電離圏遅延バイアスを遅延バイアス平均成分と遅延バイアス変動成分とに分解して推定演算し、ワイドレーン結合の観測方程式と電離圏フリー結合の観測方程式とにより推定演算するようにしたものである。

特開２００５−３２１３６２号公報

図４に示すように、航法衛星５０からの衛星信号は、それぞれ基準局５１とモニタ局５２、５３、５４で受信されるが、この場合、基準局（受信機ａ）５１とそれぞれモニタ局（受信機ｂ）５２、５３、５４との間隔（以下、基線長と記す）が、長い長基線Ａ_３あるいは中基線Ａ_２の場合は勿論のこと、短い短基線Ａ_１の場合でも電離圏異常による電離圏遅延差が大きな測位誤差となる。

例えば、長基線Ａ_３の場合の相対測位では、電離圏遅延や対流圏遅延が測位誤差として大きく現れるため、上記特許文献１に記載のように、２つの周波数の衛星信号（２周波）の観測データを利用してモデル補正や推定手法を用いることにより測位精度の信頼性を高めている。一方、短基線Ａ_１の場合の相対測位では、それらは相殺されるものと仮定される。その上、２周波観測データを用いているので、例えば、電離圏異常時にはシンチレーションという現象により受信が途切れやすく、何らかの理由によりいずれか一方の周波数の観測データが途切れた場合には、２周波を結合して同時に完全に利用する方法の場合、モニタリングができなくなるという問題がある。

又、長基線Ａ_３は勿論のこと、短基線Ａ_１の場合であっても、プラズマバブルやＳＥＤ（Ｓｔｏｒｍ Ｅｎｈａｎｃｅｄ Ｄｅｎｓｉｔｙ）等により、電離圏に異常が発生した場合、例えば、図４に示すように、電離圏の厚い箇所ｔ_１や薄い箇所ｔ_２が発生し、さらに厚い箇所ｔ_１と薄い箇所ｔ_２との間に勾配があるような場合には、電離圏異常をモニタリングすることが出来なくなり、測位誤差が増大するという問題がある。

そこで、発明者は、本願発明のように、１つの周波数の衛星信号（１周波）の観測データを利用するとともに、基線長を種々変えた場合について実験し、検討した。

例えば、本願発明のように、１周波の搬送波位相を利用するとともに、基線長を数百ｍオーダーの短基線Ａ_１とした場合、搬送波位相には搬送波位相バイアスが含まれるため、通常は搬送波位相バイアスを決定する必要があるが、基線長を短くすれば搬送波位相バイアスを決定する必要は無くなる。しかしながら、基線長が短い場合、電離圏遅延差も小さくなるため、搬送波位相の観測雑音やマルチパスを適切に調整する必要が生じる。そのため、サイティング要件が厳しくなる。さらに、基線長に完全に依存した方法となり、利用出来る範囲が非常に限定され、応用範囲が狭くなるという問題がある。

次いで、１周波のコード疑似距離を利用するとともに、基線長を１０−３０ｋｍ程度の中基線Ａ_２とした場合、コード疑似距離は、搬送波位相に比べると、観測雑音やマルチパスは１００倍程度ともなりその影響が大きい。そこで、距離（基線長）を短くしていくと、電離圏遅延差が小さくなり、電離圏異常の検出が困難となる。従って、ある程度基線長を長くした上で、観測雑音やマルチパス等の雑音と区別出来る程度に電離圏遅延差が大きくならないと、即ち、変動の大きい電離圏遅延差しか検出出来ない。従って、基線長を長くしているため、その間の電離圏異常をモニタリングすることが出来ないという問題があり、ＧＢＡＳの電離圏フィールドモニタ（Ｉｏｎｏｓｐｈｅｒｉｃ Ｆｉｅｌｄ Ｍｏｎｉｔｏｒ：以下、ＩＦＭと記す。)としては、利用出来ない可能性がある。

基線長を５０ｋｍ以上の長基線Ａ_３とした場合、対流圏遅延、航法衛星の軌道誤差等の影響も考慮する必要がある。従って、パラメータを増やす必要があり、１周波のみでは電離圏異常を検出することが困難となる。そこで、このような場合には、２周波を利用するのが一般的である。なお、２周波を使う方法は、基線長の長短に関係なく、如何なる距離でも使用されているが、上記のような問題がある。

本願発明は、上記のような問題点を鑑みなされたもので、１周波数のみの観測データを用い、基準局とモニタ局間の基線長に影響されることのない衛星航法システムにおける電離圏異常を検出する方法及びその装置を提供することを目的としている。

請求項１に係る発明は、複数の周波数で送信される航法衛星からの衛星信号を、基準局とモニタ局で受信し、この受信した衛星信号から取得した衛星データから基準局とモニタ局間の電離圏遅延の異常を検出する方法において、衛星データの１周波数のみのコード疑似距離と搬送波位相の観測データを、それぞれ基準局及びモニタ局について求め、基準局とモニタ局における１周波数のみの搬送波位相及び電離圏フリー線形結合（ＧＲＡＰＨＩＣ：ＧＲｏｕｐ Ａｎｄ Ｐｈａｓｅ Ｉｏｎｏｓｐｈｅｒｉｃ Ｃｏｒｒｅｃｔｉｏｎ）の受信機間一重差（ＳＤ：Ｓｉｎｇｌｅ Ｄｉｆｆｅｒｅｎｃｅ）を求め、同じエポックで基準局とモニタ局とで共通に観測された全航法衛星について、ＳＤ電離圏遅延、ＳＤ搬送波位相バイアス、ＳＤ受信機時計誤差を状態変数とする観測方程式を構成し、この観測方程式の状態変数を、最小二乗原則でエポック毎に推定して、観測方程式のフロート解を求め、次いで、最小二乗原則で求めたＳＤ搬送波位相バイアスを、エポック毎に他の航法衛星から求めたＳＤ搬送波位相バイアスとの間の二重差（ＤＤ：Ｄｏｕｂｌｅ Ｄｉｆｆｅｒｅｎｃｅ）であるＤＤ搬送波位相バイアスを求め、整数不定性決定手法によりエポック毎に求めたＤＤ搬送波位相バイアスを整数化した整数解の候補を求め、この求めた整数解の候補に対し検定を実行して観測方程式のフィックス解を求め、基準局とモニタ局間の推定したＳＤ電離圏遅延、ＳＤ受信機時計誤差を修正するとともに、最終的に得られたＳＤ電離圏遅延を、閾値と比較することにより電離圏異常を検出することを特徴とする衛星航法システムにおける電離圏異常を検出する方法である。

請求項２に係る発明は、請求項１に記載の発明において、基準局とモニタ局間の基線長を短基線としたものである。

請求項３に係る発明は、請求項１〜請求項２の何れかに記載の発明において、最小二乗原則は、最小二乗法を用いることとしたものである。

請求項４に係る発明は、請求項１〜請求項２の何れかに記載の発明において、最小二乗原則は、カルマンフィルタを用いることとしたものである。

請求項５に係る発明は、請求項１〜請求項４の何れかに記載の発明において、整数不定性決定手法は、ＬＡＭＢＤＡ（Ｌｅａｓｔ−ｓｑｕａｒｅｓ ＡＭＢｉｇｕｉｔｙ Ｄｅｃｏｒｒｅｌａｔｉｏｎ Ａｄｊｕｓｔｍｅｎｔ）法を適用することとしたものである。

請求項６に係る発明は、請求項１〜請求項５の何れかに記載の発明において、ＤＤ搬送波位相バイアスのフィックス解を拘束条件として、次のエポックの観測方程式に適用することとしたものである。

請求項７に係る発明は、複数の周波数で送信される航法衛星からの衛星信号を、ＧＰＳ受信機を備えた基準局とＧＰＳ受信機を備えたモニタ局で受信し、この受信した衛星信号から取得した衛星データから電離圏遅延の異常を検出する装置において、下記（１）〜（１０）の機能を備えたモニタリング装置と、（１）基準局及びモニタ局でそれぞれ受信した衛星信号から取得した衛星データの１周波数のみのコード疑似距離と搬送波位相の観測データを作成する機能と、（２）基準局とモニタ局における１周波数のみの搬送波位相及び電離圏フリー線形結合の受信機間一重差を求める機能と、（３）同じエポックで基準局とモニタ局とで共通に観測された全航法衛星について、ＳＤ電離圏遅延、ＳＤ搬送波位相バイアス、ＳＤ受信機時計誤差を状態変数とする観測方程式を構成する機能と、（４）この観測方程式の状態変数を、最小二乗原則でエポック毎に演算し、推定する機能と、（５）推定した状態変数を有する観測方程式のフロート解を求める機能と、（６）最小二乗原則で求めたＳＤ搬送波位相バイアスを、エポック毎に他の航法衛星から求めたＳＤ搬送波位相バイアスとの間の二重差であるＤＤ搬送波位相バイアスを求める機能と、（７）整数不定性決定手法によりエポック毎に求めたＤＤ搬送波位相バイアスを整数化した整数解の候補を求める機能と、（８）この整数解の候補に対し検定を実行して観測方程式のフィックス解を求める機能と、（９）基準局とモニタ局間の推定したＳＤ電離圏遅延、ＳＤ受信機時計誤差を修正する機能と、（１０）最終的に得られたＳＤ電離圏遅延と電離圏異常を判定するための閾値とを比較することにより電離圏異常を検出する機能とを備えていることを特徴とする衛星航法システムにおける電離圏異常を検出する装置である。

請求項８に係る発明は、請求項７に記載の発明において、基準局とモニタ局間の基線長を短基線としたものである。

請求項９に係る発明は、請求項７〜請求項８の何れかに記載の発明において、最小二乗原則は、最小二乗法を用いることとしたものである。

請求項１０に係る発明は、請求項７〜請求項８の何れかに記載の発明において、最小二乗原則は、カルマンフィルタを用いることとしたものである。

請求項１１に係る発明は、請求項７〜請求項１０の何れかに記載の発明において、整数不定性決定手法は、ＬＡＭＢＤＡ法を適用することとしたものである。

請求項１２に係る発明は、請求項７〜請求項１１の何れかに記載の発明において、モニタリング装置は、ＤＤ搬送波位相バイアスのフィックス解を拘束条件として、次のエポックの観測方程式を作成し、適用する機能とを備えたものである。

請求項１及び請求項７に係る発明は、上記のように構成したので、ＧＢＡＳにおいて問題とされている電離圏異常を検出することが出来る。その上、急激な測位誤差の増大を回避することが出来るので、ロバストな測位システムを実現することが出来る。又、航法衛星からの衛星信号のうち１周波数のみの衛星データを利用しているので、従来の２周波数を利用するシステムのように、電離圏異常時に発生するシンチレーション現象により衛星信号の受信が途切れるという問題もなく、電離圏異常を安定的にモニタリングすることが出来る。その上、安価なシステムを構成することが出来る。

ＳＤ電離圏フリー線形結合を用いることにより、電離圏異常時でもそれらの影響を受けず、搬送波位相バイアスを高速かつ正確に推定することが可能である。又、従来のように、単純にコード疑似距離を利用する方法に比べて、電離圏フリー線形結合は観測ノイズが低減されており、より信頼性の高い推定が出来る。搬送波位相バイアスを高速かつ正確に推定出来る。

受信機間一重差（ＳＤ）を行う事により、航法衛星側に起因する誤差要因（主に衛星時計誤差，衛星側の初期位相成分）を相殺することが出来る。又、その他の誤差要因である衛星の軌道誤差もまたほぼ相殺出来る。

さらに、他の航法衛星から求めたＳＤ搬送波位相バイアスとの間の二重差（ＤＤ）であるＤＤ搬送波位相バイアスを求めているので、受信機側の初期位相成分が相殺され、観測方程式を整数として扱うことが出来、整数不安定性決定手法を用いることが出来、パラメータを修正することが出来る。

請求項２及び請求項８に係る発明は、上記のように構成したので、請求項１及び請求項７に係る発明と同様の効果がある。さらに、基準局とモニタ局間の基線長は、短基線であるから、対流圏遅延は相殺されると仮定することが出来る。又、航法衛星の軌道誤差は、ほぼ０とすることが出来る。

請求項３及び請求項９に係る発明は、上記のように構成したので、請求項１〜請求項２及び請求項７〜請求項８に係る発明と同様の効果がある。

請求項４及び請求項１０に係る発明は、上記のように構成したので、請求項１〜請求項２及び請求項７〜請求項８に係る発明と同様の効果がある。さらに、推定値とその推定誤差共分散をエポック間で引き継ぎが可能となるので、適切な時間変化のモデルを与えれば、遂次的に最適な最小二乗原則に基づいて推定値が得られる。

請求項５及び請求項１１に係る発明は、上記のように構成したので、請求項１〜請求項４及び請求項７〜請求項１０に係る発明と同様な効果がある。さらに、ＬＡＭＢＤＡ法を適用することにより、整数解の候補が得られる。

請求項６及び請求項１２に係る発明は、上記のように構成したので、請求項１〜請求項５及び請求項７〜請求項１１に係る発明と同様な効果がある。さらに、測定精度の信頼性が向上する。

この発明の実施例を示すもので、本願発明の衛星航法システムにおける電離圏異常を検出する方法を説明するための模式図である。 この発明の実施例を示すもので、本願発明の衛星航法システムにおける電離圏異常を検出する方法を説明するためのブロック図である。 この発明の実施例を示すもので、本願発明の衛星航法システムにおける電離圏異常を検出する方法における電離圏遅延差の推定を行う際の処理フローである。 従来例を示す模式図である。

１周波数のみの観測データを用い、基準局とモニタ局間の基線長に影響されることのない衛星航法システムにおける電離圏異常を検出する方法及びその装置を提供とすることを目的として、複数の周波数で送信される航法衛星からの衛星信号を、基準局とモニタ局で受信し、この受信した衛星信号から取得した衛星データの１周波数のみのコード疑似距離と搬送波位相の観測データを、それぞれ基準局及びモニタ局について求め、基準局とモニタ局における１周波数のみの搬送波位相及び電離圏フリー線形結合の受信機間一重差を求め、同じエポックで基準局とモニタ局とで共通に観測された全航法衛星について、ＳＤ電離圏遅延、ＳＤ搬送波位相バイアス、ＳＤ受信機時計誤差を状態変数とする観測方程式を構成し、この観測方程式の状態変数を、カルマンフィルタでエポック毎に推定して、観測方程式のフロート解を求め、次いで、カルマンフィルタで求めたＳＤ搬送波位相バイアスを、エポック毎に他の航法衛星から求めたＳＤ搬送波位相バイアスとの間の二重差であるＤＤ搬送波位相バイアスを求め、ＬＡＭＢＤＡ法によりエポック毎に求めたＤＤ搬送波位相バイアスを整数化した整数解の候補を求め、この求めた整数解の候補に対し検定を実行して観測方程式のフィックス解を求め、基準局とモニタ局間の推定したＳＤ電離圏遅延、ＳＤ受信機時計誤差を修正するとともに、最終的に得られたＳＤ電離圏遅延を、閾値と比較することにより電離圏異常を検出し、ＤＤ搬送波位相バイアスのフィックス解を拘束条件として、次のエポックの観測方程式に適用することで実現する。

この発明の実施例を、図１〜図３に基づいて詳細に説明する。
図１〜図３は、この発明の実施例を示すもので、図１及び図２は本願発明の衛星航法システムにおける電離圏異常を検出する方法を説明するための模式図及びブロック図で、図３は本願発明の衛星航法システムにおける電離圏異常を検出する方法における電離圏遅延差の推定を行う際の処理フローである。

図１及び図２において、モニタリング装置１は、基準局３及びモニタ局４から送信された航法衛星２（２ａ、２ｂ・・・）の観測データをもとに電離圏異常を検出する機能を有する装置である。

航法衛星２（２ａ、２ｂ・・・）は、この実施例ではＧＰＳ衛星であり、それぞれ衛星信号（ＧＰＳ衛星自身の詳細な軌道情報等の衛星データを含む）を、複数の周波数（Ｌ１：１．５ＧＨｚ帯／Ｌ２：１．２ＧＨｚ帯等）で送信している。

基準局３及びモニタ局４は、最終的には航法衛星２（２ａ、２ｂ・・・）から送信される衛星信号（以下、単に衛星信号と記す。）による測位精度を向上させる目的で、場所が既知である複数の既知点にそれぞれ設置される。また、基準局３及びモニタ局４は、必要に応じてそれぞれ１又は複数局設置され、これら基準局３及びモニタ局４とモニタリング装置１とは、ネットワーク（図示せず）を介して接続されている。

基準局３及びモニタ局４は、ＧＰＳ衛星から送信されている複数の周波数の衛星信号の内、１つの周波数（Ｌ１）のみを受信可能な、いわゆる１周波型のＧＰＳ受信機１１及びＧＰＳ受信機１２をそれぞれ有している。基準局３及びモニタ局４は、このＧＰＳ受信機１１及びＧＰＳ受信機１２で衛星信号を受信することにより得られた観測データ１３及び観測データ１４を、ネットワークを介してモニタリング装置１へ送信する。
また、この実施例では、説明を簡単にするために、基準局３とモニタ局４との基線長は短基線としている。

ＧＰＳ受信機１１及びＧＰＳ受信機１２は、本願発明では前述の通り、いわゆる１周波型のＧＰＳ受信機を使用しているが、当然ながら２つの周波数（Ｌ１、Ｌ２）とも受信可能な２周波型ＧＰＳ受信機を有するものでも良い。また、ＧＰＳ受信機１１及びＧＰＳ受信機１２は、Ｌ１以外の周波数、例えば、Ｌ２やＬ５等のいずれか１つの周波数を受信可能なＧＰＳ受信機を有するものでも良い。

観測データ１３及び観測データ１４は、それぞれ、ＧＰＳ受信機１１及びＧＰＳ受信機１２において受信した衛星信号から取得した衛星データをもとに算出されたデータであって、航法衛星２（２ａ、２ｂ・・・）と基準局３又はモニタ局４とのコード疑似距離（以下、単にコード疑似距離と記す。）や、基準局３又はモニタ局４における衛星信号の搬送波の位相（以下、単に搬送波位相と記す。）等の情報を含むデータである。また、この実施例の場合、ＧＰＳ受信機１１及びＧＰＳ受信機１２は、前述の通り、１周波型のＧＰＳ受信機を使用しているため、観測データ１３及び観測データ１４は、１周波数のみのデータとなっている。

なお、モニタリング装置１の代わりに、基準局２やモニタ局３の一つにモニタリング装置１の機能を持たせても良い。また、基準局３とモニタ局４との基線長は中長基線であっても、本願発明の手法は適用可能である。

次に、作用動作について、図１〜図２に基づいて説明する。

まず、基準局３及びモニタ局４は、航法衛星２（２ａ、２ｂ・・・）からの衛星信号をＧＰＳ受信機１１及びＧＰＳ受信機１２でそれぞれ受信し、この受信した衛星信号に含まれる衛星データをもとに算出された観測データ１３及び観測データ１４を、ネットワークを介してモニタリング装置１へ送信する。

次に、モニタリング装置１では、このネットワークを介して基準局３及びモニタ局４から得られた観測データ１３及び観測データ１４に対して、適切な前処理を行う（図２ １５）。この前処理１５は、複数の受信機を使う場合の定常手段であって、時刻合わせ等の処理を行っている。前処理１５は、複数の受信機を使う場合、モニタリング装置１にデータを送って処理する際に、受信機自体も離れているために時間の遅延があり、また、受信機間の時間のズレが大きすぎると、航法衛星も電離圏自体も動いているので、時刻合わせ等の処理をする必要があるために行っている。

そして、この前処理１５の後に、電離圏遅延差推定部において電離圏の遅延差の推定を行い（図２ １６）、最終的にこの推定された電離圏の遅延をもとに電離圏異常を検知する（図２ １７）。

次に、電離圏遅延差推定部１６において行われる電離圏遅延差の推定の処理フローにおけるアルゴリズムを、図３に基づいて説明する。最初に、このアルゴリズムの全体の流れと概要について説明する。

まず、基準局３とモニタ局４において取得した１周波数のみの搬送波位相及びコード疑似距離を利用して（ステップ２１）、そのときの搬送波位相および電離圏フリー線形結合の受信機間一重差（ＳＤ：Ｓｉｎｇｌｅ Ｄｉｆｆｅｒｅｎｃｅ）を観測量とした観測方程式を構成する（ステップ２５）。次いで、この観測方程式に含まれる状態変数であって、推定すべきパラメータであるＳＤ電離圏遅延、ＳＤ搬送波位相バイアス及びＳＤ受信機時計誤差をエポック毎に推定する（ステップ２２）。

なお、エポックとは、ＧＰＳ受信機で用いられる時間の単位である。また、上記の電離圏遅延や搬送波位相バイアス等の各パラメータの前に付しているＳＤは、当該各パラメータの受信機間のデータの差をとった受信機間一重差を意味する。例えば、上記のＳＤ受信機時計誤差とは、受信機時計誤差の受信機間一重差という意味である。以後、特に記載のない限り、この様に記載する。

さらに、後述するＤＤ（ＤＤ：Ｄｏｕｂｌｅ Ｄｉｆｆｅｒｅｎｃｅ）は、各パラメータについて受信機間で一重差をとったもの（ＳＤ）と、さらに別の航法衛星において各パラメータについて受信機間で一重差をとったもの（ＳＤ）との差をとる航法衛星間の二重差を意味する。このＤＤについてもＳＤと同様に、各パラメータの前にＤＤを付した場合には、特に記載のない限り、当該パラメータの二重差を意味するものとして記載する。

次に、上記状態変数は、最小二乗法またはカルマンフィルタ等の最小二乗原則を用いてエポック毎に推定する（ステップ２４，２６）。このとき得られた推定値は、フロート解（Ｆｌｏａｔ解）と呼ばれる（ステップ２７）。カルマンフィルタを使うことで、推定値と、その推定誤差共分散をエポック間で引き継ぐことが可能となるので、適切な時間変化のモデルを与えてやれば、逐次的に最適な最小二乗原則での推定値が得られることとなる。

このカルマンフィルタによって得られたフロート解において、ＳＤ搬送波位相バイアスは、受信機の初期位相成分が含まれるので、基本的に整数としては扱うことが出来ない。しかしながら、基準となる航法衛星を選択した上で（ステップ２８）、エポック毎にＳＤ搬送波位相バイアスをＤＤ搬送波位相バイアスに変換することにより、整数解として扱うことが出来るようになる（ステップ２９）。

なお、本願発明のように、この二重差に変換することにより、搬送波位相バイアスを整数解として扱うことが出来る理由は、以下の通りである。まず、受信機間一重差を求める際には、各受信機で共通にある特定の航法衛星を対象としてデータを取得しているので、航法衛星側の誤差成分（衛星時計誤差、衛星軌道誤差等）を相殺することができる。この受信機間の一重差を求めた状態で、さらにデータ取得する対象の航法衛星をもう一つ増やして同様に受信機間の一重差を求める。この受信機間の一重差同士の差をさらに求めることにより、今度は受信機側を共通とすることができるので、受信機側の誤差成分を相殺することが出来る。この受信機側の誤差成分として、ＳＤ搬送波位相バイアスに含まれている受信機側の初期位相を相殺することが出来るので、整数解として扱うことが出来る。

このようにして、搬送波位相バイアスが整数解として扱うことが出来るようになれば、整数不定性決定手法を適用することが出来る。この整数不定性決定手法を適用する事により、さらにパラメータをより最適な値へ修正をすることが出来る。この実施例では、この整数不定性決定手法は、現時点において、最も実用化が進んでおり、かつ、この分野で最も利用されているＬＡＭＢＤＡ（Ｌｅａｓｔ−ｓｑｕａｒｅｓ ＡＭＢｉｇｕｉｔｙ Ｄｅｃｏｒｒｅｌａｔｉｏｎ Ａｄｊｕｓｔｍｅｎｔ）法を利用している。なお、本願発明は、ＬＡＭＢＤＡ法に限定されるものではなく、同様な手法であれば、他の方法であっても良い。

この実施例では、このＬＡＭＢＤＡ法を利用するために、カルマンフィルタによって得られた推定値及び推定誤差共分散の搬送波位相バイアスの部分をＤＤに変換して（ステップ２９）、ＬＡＭＢＤＡ法を適用する（ステップ３０）。ＬＡＭＢＤＡ法を適用することにより、整数解の候補が得られる。この整数解の候補に対して、一定の基準をもとに順列をつけた状態で算出する。

次いで、整数解の候補に対して適切な検定（Ｖａｌｉｄａｔｉｏｎ）を行い（ステップ３１）、最終的な整数解であるフィックス解（Ｆｉｘ解）を求める（ステップ３２）。このフィックス解をもとにして、最終的に得られたＳＤ電離圏遅延を逐次的にモニタリングすることにより、電離圏異常を検出する（ステップ３３）。

なお、最終的に整数解として算出されたＤＤ搬送波位相バイアスの信頼性が十分高ければ、ＤＤ搬送波位相バイアスのフィックス解を拘束条件として、次のエポックで観測方程式に組み込むことも出来る。

次に、図３に基づいて、電離圏遅延差推定部１６において行われる電離圏遅延差の推定の処理フローにおけるアルゴリズムの詳細について説明する。

基準局３の位置を（χ_ａ,ｙ_ａ,ｚ_ａ）、モニタ局４の位置を（χ_ｂ,ｙ_ｂ,ｚ_ｂ）としたとき、基準局３及びモニタ局４において観測される航法衛星２（２ａ、２ｂ・・・）の位置を（χ^ｐ,ｙ^ｐ,ｚ^ｐ）とする。この場合、基準局３及びモニタ局４における搬送波位相（距離）及びコード疑似距離は、下記の数式１で表現できる。なお、座標系はＷＧＳ８４系の地心地球固定座標系（ＥＣＥＦ：Ｅａｒｔｈ Ｃｅｎｔｅｒｅｄ Ｅａｒｔｈ Ｆｉｘｅｄ）としている。また、以後の数式等において用いられているパラメータの添字については、ａは基準局３を、ｂはモニタ局４を、ｐは航法衛星２（２ａ、２ｂ・・・）をそれぞれ表している。

上記数式１において、搬送波位相（距離）の式とコード疑似距離の式を比較すると、電離圏遅延であるδIが、互いに±となっている、即ち、極性が反対になっているのが特徴である。また、数式１における搬送波位相（距離）の式には、Ｎというバイアス成分が含まれているのも特徴である。なお、搬送波位相やコード疑似距離の観測ノイズであるε_Φ及びε_ρは、平均0の正規性白色雑音として、そのときの分散の値に関しては仰角に応じた重み付けを行う。この重み付けの手法としては、仰角が低ければノイズが上がるような係数を掛け算している。また、搬送波位相バイアスＮには、それぞれ初期位相の受信機側・航法衛星側の初期位相成分が含まれているため、整数でなく、実数となる。

ここで、基準局３と航法衛星２（２ａ、２ｂ・・・）との間の幾何学的距離は、下記の数式２で表現できる。

モニタ局４と航法衛星２（２ａ、２ｂ・・・）との間の幾何学的距離も、同様に下記の数式３で表現できる。

さらに、搬送波位相（距離）およびコード疑似距離に基づいた電離圏フリー線形結合（ＧＲＡＰＨＩＣ：ＧＲｏｕｐ Ａｎｄ ＰＨａｓｅ Ｉｏｎｏｓｐｈｅｒｉｃ Ｃｏｒｒｅｃｔｉｏｎ）は、下記の数式４で表現できる。この数式４は、搬送波位相（距離）とコード疑似距離を加算して２分割しているのが特徴的である。上記したように、数式１における搬送波位相（距離）の式とコード疑似距離の式を比較すると、電離圏遅延であるδIが互いに±となっている、即ち、極性が反対になっている。このように、極性が反対となっているので、加算することにより、この電離圏遅延であるδIを相殺出来る。しかしながら、加算することによって係数がつくことになるが、２分割することにより、通常のコード疑似距離等の成分と同様な成分の式となる。

基準局３とモニタ局４との間で、受信機間一重差（ＳＤ）を搬送波位相および電離圏フリー線形結合に対して行うと下記の数式５となる。なお、受信機間一重差を行うことにより、航法衛星側に起因する誤差要因（主に衛星時計誤差、航法衛星側の初期位相成分、衛星軌道誤差）は相殺される。

この数式５に含まれているパラメータのうち、基準局３とモニタ局４の位置は既知であるから、その距離の成分である幾何学的距離は既知として扱うことが出来る。この実施例では、基準局３とモニタ局４の基線長は短基線という条件であるので、対流圏遅延は相殺されると仮定できる。従って、δＴ_ＳＤは殆ど０として扱うことが出来る。それ故、数式５は簡単な式となり、ＳＤ搬送波位相およびＳＤ電離圏フリー線形結合は以下の数式６で表現できる。

数式６に示すように、ＳＤ電離圏フリー線形結合を用いることにより、電離圏異常時でもそれらの影響を受けず、搬送波位相バイアスの推定を容易にすることが出来る。また、単純にコード疑似距離を利用する場合に比べて、電離圏フリー線形結合は観測ノイズが低減されており、より信頼性の高い推定が可能となる。

次いで、数式６において、同じエポックで基準局３とモニタ局４とで共通に観測された全航法衛星（２ａ、２ｂ・・・）に対する式を作成すると、最終的に下記の数式７のようなベクトルと行列式で表現された式となる。なお、このとき、基準局３とモニタ局４とで共通に観測された航法衛星（２ａ、２ｂ・・・）の数は、全部でｎ個とする。

ここで、状態空間モデルを以下の数式８に示す。この状態空間モデルは、二つの式から構成されており、数式８の上の式（ｘ_ｋ）は状態方程式と呼ばれ、時間的な予測をするモデルである。数式８の下の式（ｙ_ｋ）は、観測方程式であり、上記の数式７を簡潔に表現したものである。なお、ｋはエポックを表す。

この数式８において、ｗ_ｋ及びｖ_ｋは、平均０の正規性白色雑音とする。また、数式８の観測方程式におけるＹ、Ｈ、χ、ｖは、下記の数式９で表される。

従って、カルマンフィルタのアルゴリズムは、下記の表１の様に表現出来る。

カルマンフィルタを適用するためには、数式８の上の式（ｘ_ｋ）に示す状態方程式における状態変数のダイナミクスモデルが必要である。この状態方程式におけるχの中のパラメータである電離圏遅延、受信機時計誤差及び搬送波位相バイアス等が、状態変数である。その状態変数のダイナミクスモデルについては、それぞれ複数のダイナミクスモデルが知られている。カルマンフィルタには、この状態変数それぞれに対して予測というプロセスがあり、この予測に関して、どのダイナミクスモデルを使用して予測を行うのかを決定することが重要である。従って、測定を行う地域における電離圏の状態や、使用する受信機の特性等に応じ、状態変数それぞれについて適切なダイナミクスモデルを選択する必要がある。

電離圏遅延のダイナミクスモデルには、下記の表２に示すように、ランダムウォークや1次ガウス・マルコフモデル等の状態変数のモデルがある。なお、この実施例では、1次ガウス・マルコフを適用している。

受信機時計誤差のダイナミクスモデルには、下記の表３に示すように、ホワイトノイズやランダムウォーク、1次ガウス・マルコフモデル等の状態変数のモデルがある。なお、この実施例では、ランダムウォークを適用している。

搬送波位相バイアスのダイナミクスモデルには、下記の表４に示すように、コンスタントやランダムウォーク等の状態変数のモデルがある。なお、この実施例では、ランダムウォークを適用している。

以上の表２〜表４に示した状態変数のダイナミクスモデルに基づいて、下記の数式１０に示す状態遷移行列及びシステム雑音行列が得られる。

また、カルマンフィルタは、最初に初期値を与える必要があり、その初期値については下記の数式１１に示す通りとする。

数式１１において、電離圏遅延のＩ_ＳＤ，０に関しては、その初期値は基本的には０であるが、数式１１における括弧内の式でも初期値を設定可能である。疑似距離ρ_ＳＤのパラメータの中には、幾何学的な距離ｒ_ＳＤと受信機時計誤差ｂ_ＳＤ，０と電離圏遅延Ｉ_ＳＤ，０と、その観測ノイズが含まれている。この観測ノイズには、コードマルチパスを含んでいると考えられる。幾何学的な距離ｒ_ＳＤは既知であり、受信機時計誤差も数式１１に示すように、最小の値が最小二乗法等で求められ、これらの値を疑似距離ρ_ＳＤから減算することにより、電離圏遅延Ｉ_ＳＤ，０と観測ノイズの成分だけとなる。観測ノイズの成分については、観測ノイズに含まれるコードマルチパスが大きければ、電離圏遅延の初期値Ｉ_ＳＤ，０が不適切な値に設定される可能性を含んでいるが、コードマルチパスが小さければ、数式１１における括弧内の式でも初期値を設定可能である。

搬送波位相バイアスＮ_ＳＤ，０に関しては、数式１１に示すように、搬送波位相（距離）Φ_ＳＤから幾何学的距離ｒ_ＳＤと受信機時計誤差ｂ_ＳＤ，０を減算することにより、初期値として求めることが出来る。搬送波位相（距離）Φ_ＳＤのパラメータの中には、幾何学的な距離ｒ_ＳＤと受信機時計誤差ｂ_ＳＤ，０と搬送波位相バイアスＮ_ＳＤ，０と電離圏遅延Ｉ_ＳＤ，０と、その観測ノイズが含まれている。

数式６からも明らかであるように、受信機時計誤差ｂ_ＳＤ，０は、大きな成分であるが、電離圏遅延の初期値Ｉ_ＳＤ，０を求めるときと同様に、最小の値が最小二乗法等で求められる。幾何学的な距離ｒ_ＳＤも前述の通り既知であることから、受信機時計誤差ｂ_ＳＤ，０と幾何学的な距離ｒ_ＳＤを搬送波位相（距離）Φ_ＳＤから減算することにより、搬送波位相バイアスＮ_ＳＤ，０と電離圏遅延Ｉ_ＳＤ，０と観測ノイズの成分だけとなる。

さらに、この残った成分の中で、観測ノイズは殆どないので無視することが出来る。電離圏遅延に関しても、搬送波位相バイアスのＮ_ＳＤ，０の初期値は、大体の値を設定出来れば良いので、電離圏遅延もまた殆どないと仮定できる。以上のことから、数式１１に示すように、搬送波位相バイアスのＮ_ＳＤ，０の初期値は、搬送波位相（距離）Φ_ＳＤから幾何学的距離ｒ_ＳＤと受信機時計誤差ｂ_ＳＤ，０を減算することにより求めることが出来る。

受信機時計誤差の初期値ｂ_ＳＤ，０は、数式１１や前述の通り、それぞれ基準局３とモニタ局４とにおいて、一般的な最小二乗法等で受信機の時計を求めた上でその差を取ることにより求める。この手法は、単独測位等で一般的な手法である。

なお、最初の推定の初期値以外に、例えば新しく見え始めた航法衛星、即ち、新たな可視航法衛星に関しても初期値を設定する必要があるので、それについても同様に設定する。ただし、一旦受信し、何かしらの原因で受信が途切れてしまい、また受信されるという状態もある。このような受信欠落後に、ある一定期間（数エポック〜数十エポック程度）であれば、電離圏遅延に関しては、急に極端な変動は起こり得ないために大きな差はないとして、欠落直前の推定値を保持する。その上で、受信が途切れてから一定期間後に、また受信した場合には、その値を再受信後の初期値として設定する。このようにすることで、推定時の収束を速めることが可能である。

上述のカルマンフィルタにより得られた解

は、フロート解と呼ばれる。このフロート解のＳＤ搬送波位相バイアスの推定値

は、受信機側の初期位相の成分が含まれるため、一般的に整数として扱うことが出来ないので、下記の数式１２に示すように、変換行列を用いてＤＤに変換することで、受信機側の初期位相分を相殺し、整数として扱うことが可能となる。なお、この変換行列が数式１２のＤとＴの式である。また、同様に、その時の推定誤差の共分散の値に関しても、同じように変換する。このように変換した後に、搬送波のバイアス成分とそれ以外（受信機時計誤差及び電離圏遅延に関する共分散）という分け方をして、さらにそれぞれの相関をとる形で、推定誤差の共分散を４分割にしている。

なお、ＳＤからＤＤへの変換に際しての基準となる航法衛星の選択方法については、最高仰角の航法衛星を選択するのが一般的であるが、必ずしもその航法衛星が電離圏異常の影響を受けていないとは限らないため、残差、共分散、航法衛星の配置等を利用して総合的に判断される。

上記の数式１２により、搬送波位相バイアスをＤＤに変換した上で、整数不定性決定手法であるＬＡＭＢＤＡ法を適用する。このＬＡＭＢＤＡ法を適用すると、搬送波位相バイアスの整数解

の候補が得られる。この整数解の候補には順列が付けられており、最も誤差（残差）の少ない候補を最善解

とし、その次を次善解

として、この２つの解を整数解の候補とするのが一般的である。

この得られた２つの候補に対して、それぞれの誤差を評価することにより適切に検定を行う。この検定方法には色々な方法があり、この実施例では下記の数式１３に示すレシオテストを用いている。

上記数式１３に示すような検定を行う必要がある理由は、ＬＡＭＢＤＡ法によって得られた２つの候補を比較した場合に、似たような値であれば、どちらの解を採用すれば良いのかを判断できないからである。しかし、数式１３において、その比がかけ離れていれば（閾値以上であるならば）、圧倒的に最善解の誤差が小さいこととなり、その最善解を採用することができる。

最終的に検定を通過した場合、下記の数式１４によって修正を行うことによって、フィックス解が求められる。

なお、数式１３に示す検定を通過しなかった場合は、上記数式１４による修正は行わず、フィックス解ではなく、フロート解を利用して、後述する電離圏異常を検出する。

以上の処理により、最終的に得られたフィックス解若しくはフロート解におけるＳＤ電離圏遅延を、閾値と比較する等のモニタリングを逐次的に行うことにより、電離圏異常を検出する。

さらにフィックス解として推定されて決定された搬送波位相バイアスの整数解の信頼性が、十分高いと判断出来れば、それを拘束条件として、次のエポックにおいて、数式７に示す観測方程式の一番下の式として組み込んで、推定することが可能である。この場合、次のエポックの観測方程式は、下記の数式１５となる。この数式１５では、ε_Ｎは、極めて小さな雑音を設定している。

なお、信頼性が十分高いという判断としては、一定回数連続で同一の解が得られたとき等が考えられる。このようにして、信頼性が高いかどうかを判断して、信頼性が十分高いと判断されれば、次のエポックからは、フロート解を推定する段階から、拘束条件を付けて推定を行っていくこととなる。

この発明による衛星航法システムにおける電離圏異常を検出する方法及び装置は、移動体の測位システム、誘導システム等に利用可能であり、広くＧＮＳＳ衛星を利用した測位分野に応用することが出来る。

１ モニタリング装置
２（２ａ、２ｂ・・・） 航法衛星
３ 基準局
４ モニタ局

Claims (12)

1. 複数の周波数で送信される航法衛星からの衛星信号を、基準局とモニタ局で受信し、この受信した衛星信号から取得した衛星データから前記基準局と前記モニタ局間の電離圏遅延の異常を検出する方法において、
   前記衛星データの１周波数のみのコード疑似距離と搬送波位相の観測データを、それぞれ前記基準局及び前記モニタ局について求め、
   前記基準局と前記モニタ局における前記１周波数のみの搬送波位相及び電離圏フリー線形結合（ＧＲＡＰＨＩＣ：ＧＲｏｕｐ Ａｎｄ Ｐｈａｓｅ Ｉｏｎｏｓｐｈｅｒｉｃ Ｃｏｒｒｅｃｔｉｏｎ）の受信機間一重差（ＳＤ：Ｓｉｎｇｌｅ Ｄｉｆｆｅｒｅｎｃｅ）を求め、
   同じエポックで前記基準局と前記モニタ局とで共通に観測された全航法衛星について、ＳＤ電離圏遅延、ＳＤ搬送波位相バイアス、ＳＤ受信機時計誤差を状態変数とする観測方程式を構成し、
   この観測方程式の前記状態変数を、最小二乗原則でエポック毎に推定して、前記観測方程式のフロート解を求め、
   次いで、最小二乗原則で求めた前記ＳＤ搬送波位相バイアスを、エポック毎に他の航法衛星から求めたＳＤ搬送波位相バイアスとの間の二重差（ＤＤ：Ｄｏｕｂｌｅ Ｄｉｆｆｅｒｅｎｃｅ）であるＤＤ搬送波位相バイアスを求め、
   整数不定性決定手法により前記エポック毎に求めたＤＤ搬送波位相バイアスを整数化した整数解の候補を求め、
   この求めた整数解の候補に対し検定を実行して前記観測方程式のフィックス解を求め、
   前記基準局と前記モニタ局間の推定したＳＤ電離圏遅延、ＳＤ受信機時計誤差を修正するとともに、最終的に得られたＳＤ電離圏遅延を、閾値と比較することにより電離圏異常を検出すること
   を特徴とする衛星航法システムにおける電離圏異常を検出する方法。
2. 前記基準局と前記モニタ局間の基線長を短基線としたこと
   を特徴とする請求項１に記載の衛星航法システムにおける電離圏異常を検出する方法。
3. 前記最小二乗原則は、最小二乗法を用いること
   を特徴とする請求項１〜請求項２の何れかに記載の衛星航法システムにおける電離圏異常を検出する方法。
4. 前記最小二乗原則は、カルマンフィルタを用いること
   を特徴とする請求項１〜請求項２の何れかに記載の衛星航法システムにおける電離圏異常を検出する方法。
5. 前記整数不定性決定手法は、ＬＡＭＢＤＡ（Ｌｅａｓｔ−ｓｑｕａｒｅｓ ＡＭＢｉｇｕｉｔｙ Ｄｅｃｏｒｒｅｌａｔｉｏｎ Ａｄｊｕｓｔｍｅｎｔ）法を適用すること
   を特徴とする請求項１〜請求項４の何れかに記載の衛星航法システムにおける電離圏異常を検出する方法。
6. 前記ＤＤ搬送波位相バイアスのフィックス解を拘束条件として、次のエポックの観測方程式に適用すること
   特徴とする請求項１〜請求項５の何れかに記載の衛星航法システムにおける電離圏異常を検出する方法。
7. 複数の周波数で送信される航法衛星からの衛星信号を、ＧＰＳ受信機を備えた基準局とＧＰＳ受信機を備えたモニタ局で受信し、この受信した衛星信号から取得した衛星データから電離圏遅延の異常を検出する装置において、
   下記（１）〜（１０）の機能を備えたモニタリング装置と、
   （１）前記基準局及び前記モニタ局でそれぞれ受信した衛星信号から取得した衛星データの１周波数のみのコード疑似距離と搬送波位相の観測データを作成する機能と、
   （２）前記基準局と前記モニタ局における前記１周波数のみの搬送波位相及び電離圏フリー線形結合の受信機間一重差を求める機能と、
   （３）同じエポックで前記基準局と前記モニタ局とで共通に観測された全航法衛星について、ＳＤ電離圏遅延、ＳＤ搬送波位相バイアス、ＳＤ受信機時計誤差を状態変数とする観測方程式を構成する機能と、
   （４）この観測方程式の前記状態変数を、最小二乗原則でエポック毎に演算し、推定する機能と、
   （５）前記推定した状態変数を有する前記観測方程式のフロート解を求める機能と、
   （６）最小二乗原則で求めた前記ＳＤ搬送波位相バイアスを、エポック毎に他の航法衛星から求めたＳＤ搬送波位相バイアスとの間の二重差であるＤＤ搬送波位相バイアスを求める機能と、
   （７）整数不定性決定手法により前記エポック毎に求めたＤＤ搬送波位相バイアスを整数化した整数解の候補を求める機能と、
   （８）この整数解の候補に対し検定を実行して前記観測方程式のフィックス解を求める機能と、
   （９）前記基準局と前記モニタ局間の推定したＳＤ電離圏遅延、ＳＤ受信機時計誤差を修正する機能と、
   （１０）最終的に得られたＳＤ電離圏遅延と電離圏異常を判定するための閾値とを比較することにより電離圏異常を検出する機能と、
   を備えていることを特徴とする衛星航法システムにおける電離圏異常を検出する装置。
8. 前記基準局と前記モニタ局間の基線長を短基線としたこと
   を特徴とする請求項７に記載の衛星航法システムにおける電離圏異常を検出する装置。
9. 前記最小二乗原則は、最小二乗法を用いること
   を特徴とする請求項７〜請求項８の何れかに記載の衛星航法システムにおける電離圏異常を検出する装置。
10. 前記最小二乗原則は、カルマンフィルタを用いること
    を特徴とする請求項７〜請求項８の何れかに記載の衛星航法システムにおける電離圏異常を検出する装置。
11. 前記整数不定性決定手法は、ＬＡＭＢＤＡ法を適用すること
    を特徴とする請求項７〜請求項１０の何れかに記載の衛星航法システムにおける電離圏異常を検出する装置。
12. 前記モニタリング装置は、前記ＤＤ搬送波位相バイアスのフィックス解を拘束条件として、次のエポックの観測方程式を作成し、適用する機能と
    を備えていること特徴とする請求項７〜請求項１１の何れかに記載の衛星航法システムにおける電離圏異常を検出する装置。

Images