周波数間バイアス推定装置及び周波数間バイアス推定方法
【課題】衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定する。
【解決手段】１以上の受信機２４の各々から複数の測位衛星１０ａ、１０ｂ、１０ｃの各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する総電子数算出部と、電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出部と、乱数を発生させる乱数発生部と、総電子数と総電子数モデル値に基づき周波数間バイアスを算出し、乱数に基づき周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定部とを備える。
【解決手段】１以上の受信機２４の各々から複数の測位衛星１０ａ、１０ｂ、１０ｃの各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する総電子数算出部と、電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出部と、乱数を発生させる乱数発生部と、総電子数と総電子数モデル値に基づき周波数間バイアスを算出し、乱数に基づき周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定部とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、衛星から放送される複数の異なる周波数信号による通過経路の電離層総電子数を推定するための周波数間バイアス推定装置及び周波数間バイアス推定方法に関する。
【背景技術】
【０００２】
衛星ナビゲーションシステム「Ｇａｌｉｌｅｏ」やＧＰＳシステム、ＱＺＳＳ（準天頂衛星システム）等の衛星から送信される信号は、地球上の観測点の測位を行うのに利用される。その際、位置のわかっている衛星から観測点までの距離を求めるにあたって、電波の伝搬速度及び要した時間は、重要な要素となる。しかしながら、衛星から地上に向けて送信される電波信号は、電離層を通過する際に伝搬遅延量を生じ、測位誤差に大きな影響を与える。したがって、電離層伝搬遅延量を正確に求めることは、高精度な測位を行うにあたって必要不可欠なものである。
【０００３】
従来から、電離層遅延量を正確に求めるために、様々な手法が提案されている。例えば特許文献１には、正確な電離層遅延量を計算するために必要な周波数間バイアスを算出する周波数間バイアス算出装置及び方法が記載されている。周波数間バイアスは、周波数間での電気的な経路差によるものであり、受信機、衛星等のハードウェアごとに異なる。正確な電離層遅延量を計算するためには、周波数間バイアスを除去することが必要である。
【０００４】
特許文献１に記載された周波数間バイアス算出装置は、ＧＰＳ衛星から複数のＧＰＳ受信機までの第一の周波数による各第一の擬似距離と第二の周波数による各第二の擬似距離とを、当該複数のＧＰＳ受信機から入力するデータ収集部と、このデータ収集部から入力した複数の第一及び第二の擬似距離を所定の演算式に代入して第一及び第二の周波数による周波数間バイアスを算出する周波数間バイアス算出部と、この周波数間バイアス算出部で算出された周波数間バイアスをＧＰＳ受信機へ出力するデータ出力部とを備える。この周波数間バイアス算出装置によれば、簡単な構成で周波数間バイアスをユーザに提供することができる。また、大量の演算を周波数間バイアス算出装置側で行うことにより、ＧＰＳ受信機側すなわちユーザ側の負担を軽減できる。
【０００５】
電離層伝搬遅延量は、伝搬する信号の周波数に依存している。そのため、衛星からの複数の周波数信号を受信することにより、信号の通過経路における総電子数（ＴＥＣ：Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ）を得ることができる。
【０００６】
また、非特許文献１や非特許文献２には、電離層における電子密度モデル関数であるＩＲＩ（Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ）モデルが記載されている。さらに非特許文献３には、やはり電離層電子密度モデル関数の一つであるＧａｌｌａｇｈｅｒのモデルが記載されている。これらのモデル関数を用いて計算によりＴＥＣを求めることもできる。
【特許文献１】特開２００６−２３１４４号公報
【非特許文献１】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ，ｅｔ，ａｌ．，：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ １９９０”，Ｎｏｖｅｍｂｅｒ，１９９０．
【非特許文献２】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ ２０００”，Ｒａｄｉｏ．Ｓｃｉｅｎｃｅ，Ｖｏｌ．３６，Ｎｕｍｂｅｒ２，ＰＰ２６１−２７５，Ｍａｒｃｈ／Ａｐｒｉｌ，２００１
【非特許文献３】Ｇａｌｌａｇｈｅｒ，Ｄ．Ｌ．，Ｐ．Ｄ．Ｃｒａｖｅｎ，ａｎｄ Ｒ．Ｈ．Ｃｏｍｆｏｒｔ，Ｇｌｏｂａｌ ｃｏｒｅ ｐｌａｓｍａ ｍｏｄｅｌ，Ｊ．Ｇｅｏｐｈｙｓ．Ｒｅｓ．１０５，Ａ８，１８，８１９−１８，８３３，２０００．
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、衛星のバイアス及び受信機のバイアスを未知数とする方程式は、１観測地点・観測時刻に対して衛星数と同数存在するのに対して、未知数は、衛星数＋１個（各衛星のバイアス及び受信機のバイアス）存在するので、連立方程式により単純に解を求めることはできない。
【０００８】
ところが、衛星からの複数の周波数信号を受信して信号の通過経路におけるＴＥＣを求める際には、衛星・受信機固有のバイアスが追加されてしまうため、衛星及び受信機のバイアスを求めることは必要不可欠である。
【０００９】
また、上述したような電離層電子密度モデル関数は、月平均レベルのものであるため、観測結果との誤差が存在する。
【００１０】
本発明は上述した従来技術の問題点を解決するもので、衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定する周波数間バイアス推定装置及び周波数間バイアス推定方法を提供することを課題とする。
【課題を解決するための手段】
【００１１】
本発明に係る周波数間バイアス推定装置は、上記課題を解決するために、複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する周波数間バイアス推定装置であって、前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出部と、乱数を発生させる乱数発生部と、前記総電子数算出部により算出された総電子数と前記モデル値算出部により算出された総電子数モデル値に基づき周波数間バイアスを算出し、前記乱数発生部により発生した乱数に基づき前記周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定部とを備えることを特徴とする。
【発明の効果】
【００１２】
本発明によれば、衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定するので、電離層総電子数推定精度を向上し、測位精度や電離層モデルの改善に寄与することができる。
【発明を実施するための最良の形態】
【００１３】
以下、本発明の周波数間バイアス推定装置の実施の形態を、図面に基づいて詳細に説明する。
【実施例１】
【００１４】
以下、本発明の実施例について図面を参照しながら説明する。図１は、本発明の実施例１の周波数間バイアス推定装置３２及び測位システムの構成を示すブロック図である。この周波数間バイアス推定装置３２は、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される衛星信号を受信する１以上の衛星信号受信機２４を有し且つ衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する。
【００１５】
まず、本実施の形態の構成を説明する。本実施例の周波数間バイアス推定装置３２及び測位システムは、図１に示すように、測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定処理系３０、及びインターネットデータ処理系４０で構成されている。測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、及び衛星信号処理系２０は、本発明の測位システムに対応する。なお、インターネットデータ処理系４０も含めて測位システムとすることも可能であるが、インターネットデータ処理系４０は、必ずしも必須のものではなく、付加的なものである。また、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定処理系３０、及びインターネットデータ処理系４０は、お互いに通信回線５０で接続されている。
【００１６】
測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃは、ＧＰＳ、Ｇａｌｉｌｅｏ、準天頂衛星（ＱＺＳＳ）等の航法衛星であり、複数の周波数の衛星信号を送信する。
【００１７】
衛星信号処理系２０は、アンテナ２２、衛星信号受信機２４、及び衛星信号処理装置２６から構成される。衛星信号受信機２４は、本発明の受信機に対応し、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される複数の周波数の衛星信号をアンテナ２２を介して受信する。また、衛星信号処理装置２６は、衛星信号に含まれる測位情報を用いて位置情報を得る。
【００１８】
電離層ＴＥＣ・バイアス推定処理系３０は、周波数間バイアス推定装置３２を有し、受信した信号に基づき衛星位置、信号通過経路の電離層総電子数（ＴＥＣ）、周波数間バイアス等を算出する。受信したデータ・処理結果は、データサーバ系へＬＡＮ経由で伝送され保存される。
【００１９】
インターネットデータ処理系４０は、ルータ４２、ＧＥＯＮＥＴ収集データ処理装置４４、及び外部インターネット網４６で構成されている。ルータ４２は、スイッチングハブでもよい。インターネットデータ処理系４０は、公開されている電離層関連の情報や国土地理院が公開しているＧＰＳ観測データ（ＧＥＯＮＥＴデータ）や国際的に観測結果を公開しているＩＧＳ（Ｉｎｔｅｒｎａｔｉｏｎａｌ ＧＰＳ Ｓｅｒｖｉｃｅ ｆｏｒ Ｇｅｏｄｙｎａｍｉｃｓ）データ等をインターネット経由で収集する装置である。インターネットデータ処理系４０により収集された結果は、電離層ＴＥＣ・バイアス推定処理系３０で処理される。ルータ４２は、セキュリティを考慮して設けられ、ファイアウォールとする。
【００２０】
図２は、本発明の実施例１の周波数間バイアス推定装置３２の詳細なブロック図である。図２に示すように、周波数間バイアス推定装置３２は、総電子数算出部３３、モデル値算出部３４、周波数間バイアス推定部３５、及び乱数発生部３６で構成されている。
【００２１】
総電子数算出部３３は、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。具体的には、総電子数算出部３３は、衛星信号受信機２４から測位衛星１０ａ、１０ｂ、１０ｃの各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。また、総電子数算出部３３は、外部インターネット網４６から得た受信機と衛星との擬似距離に関するデータに基づき衛星信号の通過経路の総電子数を算出することもできる。総電子数の算出方法については、後述する。
【００２２】
モデル値算出部３４は、上述したようなＩＲＩモデルやＧａｌｌａｇｈｅｒモデル等の電離層電子密度モデル関数に基づき総電子数モデル値を算出する。
【００２３】
乱数発生部３６は、乱数を発生させる。
【００２４】
周波数間バイアス推定部３５は、総電子数算出部３３により算出された総電子数とモデル値算出部３４により算出された総電子数モデル値に基づき周波数間バイアスを算出し、乱数発生部３６により発生した乱数に基づき周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定する。さらに、周波数間バイアス推定部３６は、予め精密に測量された位置に基づく位置情報と、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとを真の値であると推定する。
【００２５】
次に、上述のように構成された本実施の形態の作用を説明する。なお、電離層がダイナミックに変化するような時間帯であると推定がより複雑になるため、本実施例において、周波数間バイアス推定装置は、ローカル時間の深夜２時頃に測定したデータを用いるものとする。図３は、本実施例の周波数間バイアス推定装置３２の動作を示すフローチャート図である。まず、モデル値算出部３４は、電離層電子密度モデル関数に基づき総電子数モデル値を算出する（ステップＳ１０１）。ここで、電離層電子密度モデル関数をｆとすると、通過経路の総電子数モデル値は、通過経路の電子密度を積分することにより求められる。
【数１】
【００２６】
（１）式は、衛星Ｍｇから受信観測点Ｒまでに渡り、電離層電子密度関数ｆを経路に沿って積分した結果を示す。ここで、ｓは、位置を示すパラメータである。また、ｇは、衛星を識別する識別子である。ＴＥＣｇｔｒｕｅは、通過経路における真の総電子数を示す。Δｍｏｄｅｌは、電離層密度分布モデルに含まれる誤差を示す。
【００２７】
モデル値算出部３４は、算出した総電子数モデル値を周波数間バイアス推定部３５に出力する。
【００２８】
次に、総電子数算出部３３は、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する（ステップＳ１０３）。ここでは、測位衛星１０ａと衛星信号受信機２４との間の擬似距離について考える。衛星信号通過経路のＴＥＣは、測位衛星１０ａから送信される衛星信号（ここではＬ１とＬ２で示す）に基づき求められる。まず、擬似距離（コード距離、シュードレンジ）と位相距離（フェーズ距離）は、以下のように表すことができる。
【数２】
【００２９】
ここで、ρは、擬似距離を示す。また、Φは、位相距離を示す。ｒは、真の距離を示す。ｃは光速である。また、δｔｕは、受信機時刻誤差を示し、δｔｓは、衛星時刻誤差を示す。なお、本実施例において、右下の添え字は基本的に地上（受信機等）に関連する項を示し、右上の添え字は基本的に上空（衛星等）に関連する項を示す。δｔｕ、Ｌ１ｏｒＬ２ ｂｉａｓは、受信機周波数依存ハードウェア依存バイアスを示し、δｔｓＬ１ｏｒＬ２ ｂｉａｓは、衛星周波数依存ハードウェア依存バイアスを示す。さらに、Ｉは、電離層伝搬遅延量を示し、Ｔは、対流圏伝搬遅延量を示す。Ｎａｍｂは、整数不確定値を示し、εは、観測誤差を示す。
【００３０】
総電子数算出部３３は、２つの周波数の観測値の差をとることにより、ＴＥＣを求めることができる。
【数３】
【００３１】
ここで、ｆは周波数である。また、λは波長を示す。Δｂｉａｓは、周波数間バイアスを示す。また、Δ´ｂｉａｓは、位相についてのバイアスである。
【００３２】
（４）式より、ＴＥＣは、次式のように求められる。
【数４】
【００３３】
ここで、ＴＥＣｔｒｕｅは、真の総電子数を示す。（６）式に示すＴＥＣの推定値は、バイアスに関連した項を含む。このままでは誤差ε´の値が大きいため、総電子数算出部３３は、キャリアスムージングを実施して誤差ε´の値を小さくする。キャリアスムージングとは、ノイズを多く含む擬似距離に精密な搬送波差分を適用し、ノイズの低減を図るものである。キャリアスムージング後のＴＥＣを＜ＴＥＣ＞と表すものとすると、＜ＴＥＣ＞は、以下の式により求められる。
【数５】
【００３４】
ここで、ｍは、データ収集の時間順につけた番号である。また、Ｋは、スムージングの定数であり、適宜変更できるものとする。Ｋは、サンプリング時間間隔にも依存する。例えば、時定数を１８０秒にすると、Ｋ＝１８０／ｄｔとなる。ｄｔは、サンプリング時間間隔を表す。＜Δｂｉａｓ＞は、１ヶ月程度の期間一定であると仮定し、データ収集期間中一定であるものとする。したがって、＜Δｂｉａｓ＞は、添え字ｍを必要としない。
【００３５】
総電子数算出部３３は、算出した総電子数を周波数間バイアス推定部３５に出力する。
【００３６】
次に、周波数間バイアス推定部３５は、総電子数算出部３３により算出された総電子数とモデル値算出部３４により算出された総電子数モデル値に基づき周波数間バイアスを算出する。まず周波数間バイアス推定部３５は、衛星観測のデータに基づき（８）式を用いて総電子数算出部３３により求められた総電子数と、（１）式を用いてモデル値算出部３４により求められた総電子数モデル値との差を求める（ステップＳ１０５）。
【数６】
【００３７】
（１）式で用いたような既存の電離層電子密度モデルは、月平均では観測値とよい一致を示す。そのため、（９）式のΔｇ，ｍｏｄｅｌは、１ヶ月間のデータで平均することにより、ゼロに近い値となる。さらに、衛星や受信機の周波数間バイアスは、時間的変動が小さく、時定数を１ヶ月程度と考え、この間に収集したデータについては一定であるとする。１ヶ月間の平均操作を［］で表すとすると、（９）式に示す観測ＴＥＣとモデルＴＥＣとの差は、以下のようになる。
【数７】
【００３８】
ここで、ｄεは、平均後にも残る微小なモデル誤差を示す。本実施例において、ｄεは、ゼロであるものとする。衛星依存周波数間バイアスと受信機依存周波数間バイアスとの間に相関は無いと考えられる。したがって、観測ＴＥＣとモデルＴＥＣとの差は、（１０）式に示すように２つの独立した項にまとめることができる。Ｂｇは、衛星ｇに関するバイアスを示し、ＢＲは、受信機に関するバイアスを示す。受信機依存周波数間バイアスと衛星ｇのバイアスとが、（１０）式により求められる。したがって、受信機のバイアスが決まれば、各衛星のバイアスが決まる。
【００３９】
（１０）式は、１観測地点・観測時刻に対して、衛星数と同数できる。ここで、衛星数をＮｓａｔとすると、未知数は、観測場所受信機のバイアスを含む（Ｎｓａｔ＋１）個ある。したがって、単純に連立方程式により解を求めることはできない。
【００４０】
そこで、周波数間バイアス推定部３５は、（１０）式の受信機に関連した項を、乱数発生部３６により発生した乱数により仮に決定する（ステップＳ１０７）。さらに、周波数間バイアス推定部３５は、受信機依存周波数間バイアスを仮に決定したので、（１０）式より衛星に関するバイアス成分を決定することができる。
【００４１】
次に、決定したバイアス値に基づき測位が行われる（ステップＳ１０９）。測位を行う主体は、測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、及び衛星信号処理系２０からなる測位システムであり、特に衛星信号処理装置２６である。ただし、周波数間バイアス推定装置３２が測位を行う構成でもよい。測位場所の真値をＸ０とする。Ｘ０は、あらかじめ精密に測量されている。この真値Ｘ０に対して、仮に決定したバイアス値に基づいた測位の誤差をΔＸとする。衛星と受信機との距離の誤差であるΔＲは、以下の式のように表せる。
【数８】
【００４２】
ここで、（１３）式の＜ρ＞は、（８）式のように、キャリアスムージングを行った後の擬似距離を示す。（１１）式において、ｒＲｅｐｈは、精密に測量された場所Ｘ０と衛星の位置から計算された幾何学的距離を示す。この衛星の位置は、ＩＧＳより提供されるＧＰＳ位置データや衛星放送エフェメリスデータから計算されるＧＰＳ位置データである。また、Ｔｅｓｔｇは、対流圏遅延量推定値を示す。この対流圏遅延量を推定するに際し、対流圏遅延量モデルが使用されている。＜ＩＬ１，ｅｓｔｇ＞は、電離層遅延量推定値を示す。また、ΔｔＲ，ｃｌｏｃｋは、受信機のクロックオフセットを示す。衛星のクロックオフセット（Δｔｓ，ｃｌｏｃｋ）は、衛星放送エフェメリスデータにより補正されている。この電離層遅延量推定値＜ＩＬ１，ｅｓｔｇ＞は、バイアスを含む。受信機で生じるバイアスは、各衛星共通であるため、受信機のクロックオフセットと同じ影響を持ち、測位位置誤差に影響を与えない。
【００４３】
（１１）式に示す距離誤差ΔＲと測位誤差ΔＸとの関係は、以下のように示される。
【数９】
【００４４】
ここで、Ａは、精密に測量された場所Ｘ０から衛星までの方向ベクトルを要素とする行列である。（ｌ，ｍ，ｎ）が方向ベクトルを示す。添え字は、衛星番号を示す。したがって、（１７）式（１８）式において、衛星は、Ｎ個存在する。測位位置に関して、ΔＸは、最小二乗法等により求められる。（１９）式は、最小二乗法を用いた計算例である。ΔＸの４番目の項は、受信機クロックオフセットに関連した項である。この項は、（１１）式の２番目の括弧内の項に関係するため、（１１）式の２番目の括弧内の項が反映されたものとなっていて、測位位置に影響しない。（１１）式の１番目の項は、測位位置誤差に影響を与える。衛星依存周波数間バイアスの推定値が正しければ、ΔＸの位置に関する１〜３番目の項の大きさ（誤差）は、小さくなる。
【００４５】
（１０）式の受信機に関連した項を、乱数により決定することで、衛星依存周波数間バイアス値が仮に決定される。ここで、乱数により決定した受信機依存周波数間バイアスをＢＲ，ｒａｎｄとすると、仮に決定した衛星依存周波数間バイアス値Ｂｃａｎｄｇは、以下の式により求められる。
【数１０】
【００４６】
したがって、ＢＲ，ｒａｎｄと、（２０）式により求められたＢｃａｎｄｇとを（１１）式に代入することで、擬似距離誤差ΔＲは、以下のようになる。
【数１１】
【００４７】
ここで、（２１）式の１番目の括弧内にある衛星バイアス分が真の値に近い場合に、ΔＸの位置誤差に関する大きさは、最小になる。ΔＸの位置誤差Ｐｅｒｒは、例えば位置誤差のデカルト座標成分から以下のように表すことができる。
【数１２】
【００４８】
周波数間バイアス推定部３５は、乱数発生部３６により発生した乱数に基づき仮の受信機依存周波数間バイアスＢＲ，ｒａｎｄを決定する。このＢＲ，ｒａｎｄが変わる度に（２２）式の位置誤差Ｐｅｒｒは、変化する。Ｐｅｒｒが最小の値をとる場合に、周波数間バイアス推定部３５は、その際の受信機依存周波数間バイアスＢＲ，ｒａｎｄを真の値であると推定する（ステップＳ１１１）。さらに、真の値であると推定されたＢＲ，ｒａｎｄに基づき、衛星依存周波数間バイアスＢｅｓｔｇは、（２０）式を用いて算出される（ステップＳ１１３）。
【００４９】
乱数で推定値を最後まで求めることは、効率が悪い。そこで、周波数間バイアス推定部３５は、予め測量された位置に基づく位置情報と、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき位置を測位して得られた位置情報との誤差の変化率に基づき受信機依存周波数間バイアスと衛星依存周波数間バイアスとを推定する。具体的には、乱数発生部３６により発生した乱数に基づき、周波数間バイアス推定部３５は、受信機依存周波数間バイアスをｋ回生成（ＢＲ，１，・・・，ＢＲ，ｋ）し、測位位置誤差（Ｐｅｒｒ１，・・・，Ｐｅｒｒｋ）を求める。さらに、周波数間バイアス推定部３５は、受信機に関連した項の変化率に対する測位位置誤差Ｐｅｒｒの変化率を計算する。
【数１３】
【００５０】
（２３）式により、ＢＲに対する測位位置誤差Ｐｅｒｒの変化率ηｋが求められる。このηｋに基づき、受信機依存周波数間バイアスＢＲ，ｅｓｔは、以下のように求められる。
【００５１】
まず、ηｋの符号が、（ＢＲ，１，・・・，ＢＲ，ｋ）の範囲で全てプラス又はマイナスである場合、ＢＲ，ｅｓｔは、以下のようになる。
【数１４】
【００５２】
周波数間バイアス推定部３５は、（２４）式により求められたＢＲ，ｅｓｔを真の受信機依存周波数間バイアスとして推定する。さらに詳細に受信機のバイアス値を求める場合には、周波数間バイアス推定部３５は、（２４）式により求められたＢＲ，ｅｓｔの値の近辺で、再度乱数で受信機に関連した項（ＢＲ，１，・・・，ＢＲ，ｋ）を求め、（２３）式を計算する。
【００５３】
次に、ηｋの符号が、（ＢＲ，１，・・・，ＢＲ，ｋ）の範囲内でマイナスからプラスに変化する場合、位置誤差Ｐｅｒｒの最小値は、符号が変わった近辺に存在する。そこで、符号がかわる 前後のＢＲ，ｋを、ＢＲ，ｋ−、ＢＲ，ｋ＋とすると、求める受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔは、以下のように求められる。
【数１５】
【００５４】
ここで、ηｋ＋、ηｋ−は、それぞれ、プラス符号を持つＢＲ，ｋ＋における変化率とマイナス符号を持つＢＲ，ｋ−の変化率を示す。
【００５５】
また、単純に中間値を推定値とする方法もある。その場合、受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔは、以下のように求められる。
【数１６】
【００５６】
上述したような方法により、受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔが求められると、（１０）式の関係により、各衛星の衛星依存周波数間バイアスＢＲ，ｅｓｔｇは、以下のようにして求められる。
【数１７】
【００５７】
以上のように求められた周波数間バイアスは、特定の受信機Ｒにおいて計算した受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔと各衛星の衛星依存周波数間バイアスＢＲ，ｅｓｔｇである。したがって、この推定結果は、観測地点の地形によるマルチパスの影響などの影響を含む可能性がある。
【００５８】
そこで、その影響を小さくするため、周波数間バイアス推定部３５は、測位システムが受信機（図１における衛星信号処理系２０）を複数有する場合に、複数の受信機の各々で得られた位置情報に基づき１つの測位衛星に対して得られた複数の衛星依存周波数間バイアスを平均化することにより衛星依存周波数間バイアスを推定する。具体的には、衛星依存周波数間バイアスの平均値Ｂｅｓｔｇは、以下の式により求められる。
【数１８】
【００５９】
ここで、Wは処理を行った受信機の数を示している。以上のようにして、周波数間バイアス推定装置３２は、周波数間バイアス（受信機依存周波数間バイアス及び衛星依存周波数間バイアス）を求めることができる。また、以上述べた動作は、周波数間バイアス推定方法に適用可能である。
【００６０】
上述のとおり、本発明の実施例１の形態に係る周波数間バイアス推定装置によれば、乱数に基づき受信機依存周波数間バイアス及び衛星依存周波数間バイアスを決定するので、単純に連立方程式により求めることのできない周波数間バイアスを推定することができる。
【００６１】
また、予め精密に測量された位置に基づく位置情報と、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとを真の値であると推定するので、衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定することができる。
【００６２】
さらに、精度良く推定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき、電離層総電子数推定精度を向上し、ＧＰＳ等の測位システムにおける測位精度や電離層モデルの改善に寄与することができる。また、ＨＦ帯等の周波数の伝搬経路を算出する際に使用する電離層モデルの精度を向上させることができ、より有効な通信を可能にすることができる。また、衛星を使った航法の向上が可能である。
【産業上の利用可能性】
【００６３】
本発明に係る周波数間バイアス推定装置は、測位衛星から送信される衛星信号を利用して測位を行うＧＰＳ等の測位システムに利用可能である。
【図面の簡単な説明】
【００６４】
【図１】本発明の実施例１の形態の周波数間バイアス推定装置及び測位システムの構成を示すブロック図である。
【図２】本発明の実施例１の形態の周波数間バイアス推定装置の詳細なブロック図である。
【図３】本発明の実施例１の形態の周波数間バイアス推定装置の動作を示すフローチャート図である。
【符号の説明】
【００６５】
１０ａ、１０ｂ、１０ｃ 測位衛星
２０ 衛星信号処理系
２２ アンテナ
２４ 衛星信号受信機
２６ 衛星信号処理装置
３０ 電離層ＴＥＣ・バイアス推定処理系
３２ 周波数間バイアス推定装置
３３ 総電子数算出部
３４ モデル値算出部
３５ 周波数間バイアス推定部
３６ 乱数発生部
４０ インターネットデータ処理系
４２ ルータ
４４ ＧＥＯＮＥＴ収集データ処理装置
４６ 外部インターネット網
５０ 通信回線
【技術分野】
【０００１】
本発明は、衛星から放送される複数の異なる周波数信号による通過経路の電離層総電子数を推定するための周波数間バイアス推定装置及び周波数間バイアス推定方法に関する。
【背景技術】
【０００２】
衛星ナビゲーションシステム「Ｇａｌｉｌｅｏ」やＧＰＳシステム、ＱＺＳＳ（準天頂衛星システム）等の衛星から送信される信号は、地球上の観測点の測位を行うのに利用される。その際、位置のわかっている衛星から観測点までの距離を求めるにあたって、電波の伝搬速度及び要した時間は、重要な要素となる。しかしながら、衛星から地上に向けて送信される電波信号は、電離層を通過する際に伝搬遅延量を生じ、測位誤差に大きな影響を与える。したがって、電離層伝搬遅延量を正確に求めることは、高精度な測位を行うにあたって必要不可欠なものである。
【０００３】
従来から、電離層遅延量を正確に求めるために、様々な手法が提案されている。例えば特許文献１には、正確な電離層遅延量を計算するために必要な周波数間バイアスを算出する周波数間バイアス算出装置及び方法が記載されている。周波数間バイアスは、周波数間での電気的な経路差によるものであり、受信機、衛星等のハードウェアごとに異なる。正確な電離層遅延量を計算するためには、周波数間バイアスを除去することが必要である。
【０００４】
特許文献１に記載された周波数間バイアス算出装置は、ＧＰＳ衛星から複数のＧＰＳ受信機までの第一の周波数による各第一の擬似距離と第二の周波数による各第二の擬似距離とを、当該複数のＧＰＳ受信機から入力するデータ収集部と、このデータ収集部から入力した複数の第一及び第二の擬似距離を所定の演算式に代入して第一及び第二の周波数による周波数間バイアスを算出する周波数間バイアス算出部と、この周波数間バイアス算出部で算出された周波数間バイアスをＧＰＳ受信機へ出力するデータ出力部とを備える。この周波数間バイアス算出装置によれば、簡単な構成で周波数間バイアスをユーザに提供することができる。また、大量の演算を周波数間バイアス算出装置側で行うことにより、ＧＰＳ受信機側すなわちユーザ側の負担を軽減できる。
【０００５】
電離層伝搬遅延量は、伝搬する信号の周波数に依存している。そのため、衛星からの複数の周波数信号を受信することにより、信号の通過経路における総電子数（ＴＥＣ：Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ）を得ることができる。
【０００６】
また、非特許文献１や非特許文献２には、電離層における電子密度モデル関数であるＩＲＩ（Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ）モデルが記載されている。さらに非特許文献３には、やはり電離層電子密度モデル関数の一つであるＧａｌｌａｇｈｅｒのモデルが記載されている。これらのモデル関数を用いて計算によりＴＥＣを求めることもできる。
【特許文献１】特開２００６−２３１４４号公報
【非特許文献１】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ，ｅｔ，ａｌ．，：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ １９９０”，Ｎｏｖｅｍｂｅｒ，１９９０．
【非特許文献２】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ ２０００”，Ｒａｄｉｏ．Ｓｃｉｅｎｃｅ，Ｖｏｌ．３６，Ｎｕｍｂｅｒ２，ＰＰ２６１−２７５，Ｍａｒｃｈ／Ａｐｒｉｌ，２００１
【非特許文献３】Ｇａｌｌａｇｈｅｒ，Ｄ．Ｌ．，Ｐ．Ｄ．Ｃｒａｖｅｎ，ａｎｄ Ｒ．Ｈ．Ｃｏｍｆｏｒｔ，Ｇｌｏｂａｌ ｃｏｒｅ ｐｌａｓｍａ ｍｏｄｅｌ，Ｊ．Ｇｅｏｐｈｙｓ．Ｒｅｓ．１０５，Ａ８，１８，８１９−１８，８３３，２０００．
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、衛星のバイアス及び受信機のバイアスを未知数とする方程式は、１観測地点・観測時刻に対して衛星数と同数存在するのに対して、未知数は、衛星数＋１個（各衛星のバイアス及び受信機のバイアス）存在するので、連立方程式により単純に解を求めることはできない。
【０００８】
ところが、衛星からの複数の周波数信号を受信して信号の通過経路におけるＴＥＣを求める際には、衛星・受信機固有のバイアスが追加されてしまうため、衛星及び受信機のバイアスを求めることは必要不可欠である。
【０００９】
また、上述したような電離層電子密度モデル関数は、月平均レベルのものであるため、観測結果との誤差が存在する。
【００１０】
本発明は上述した従来技術の問題点を解決するもので、衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定する周波数間バイアス推定装置及び周波数間バイアス推定方法を提供することを課題とする。
【課題を解決するための手段】
【００１１】
本発明に係る周波数間バイアス推定装置は、上記課題を解決するために、複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する周波数間バイアス推定装置であって、前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出部と、乱数を発生させる乱数発生部と、前記総電子数算出部により算出された総電子数と前記モデル値算出部により算出された総電子数モデル値に基づき周波数間バイアスを算出し、前記乱数発生部により発生した乱数に基づき前記周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定部とを備えることを特徴とする。
【発明の効果】
【００１２】
本発明によれば、衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定するので、電離層総電子数推定精度を向上し、測位精度や電離層モデルの改善に寄与することができる。
【発明を実施するための最良の形態】
【００１３】
以下、本発明の周波数間バイアス推定装置の実施の形態を、図面に基づいて詳細に説明する。
【実施例１】
【００１４】
以下、本発明の実施例について図面を参照しながら説明する。図１は、本発明の実施例１の周波数間バイアス推定装置３２及び測位システムの構成を示すブロック図である。この周波数間バイアス推定装置３２は、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される衛星信号を受信する１以上の衛星信号受信機２４を有し且つ衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する。
【００１５】
まず、本実施の形態の構成を説明する。本実施例の周波数間バイアス推定装置３２及び測位システムは、図１に示すように、測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定処理系３０、及びインターネットデータ処理系４０で構成されている。測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、及び衛星信号処理系２０は、本発明の測位システムに対応する。なお、インターネットデータ処理系４０も含めて測位システムとすることも可能であるが、インターネットデータ処理系４０は、必ずしも必須のものではなく、付加的なものである。また、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定処理系３０、及びインターネットデータ処理系４０は、お互いに通信回線５０で接続されている。
【００１６】
測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃは、ＧＰＳ、Ｇａｌｉｌｅｏ、準天頂衛星（ＱＺＳＳ）等の航法衛星であり、複数の周波数の衛星信号を送信する。
【００１７】
衛星信号処理系２０は、アンテナ２２、衛星信号受信機２４、及び衛星信号処理装置２６から構成される。衛星信号受信機２４は、本発明の受信機に対応し、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される複数の周波数の衛星信号をアンテナ２２を介して受信する。また、衛星信号処理装置２６は、衛星信号に含まれる測位情報を用いて位置情報を得る。
【００１８】
電離層ＴＥＣ・バイアス推定処理系３０は、周波数間バイアス推定装置３２を有し、受信した信号に基づき衛星位置、信号通過経路の電離層総電子数（ＴＥＣ）、周波数間バイアス等を算出する。受信したデータ・処理結果は、データサーバ系へＬＡＮ経由で伝送され保存される。
【００１９】
インターネットデータ処理系４０は、ルータ４２、ＧＥＯＮＥＴ収集データ処理装置４４、及び外部インターネット網４６で構成されている。ルータ４２は、スイッチングハブでもよい。インターネットデータ処理系４０は、公開されている電離層関連の情報や国土地理院が公開しているＧＰＳ観測データ（ＧＥＯＮＥＴデータ）や国際的に観測結果を公開しているＩＧＳ（Ｉｎｔｅｒｎａｔｉｏｎａｌ ＧＰＳ Ｓｅｒｖｉｃｅ ｆｏｒ Ｇｅｏｄｙｎａｍｉｃｓ）データ等をインターネット経由で収集する装置である。インターネットデータ処理系４０により収集された結果は、電離層ＴＥＣ・バイアス推定処理系３０で処理される。ルータ４２は、セキュリティを考慮して設けられ、ファイアウォールとする。
【００２０】
図２は、本発明の実施例１の周波数間バイアス推定装置３２の詳細なブロック図である。図２に示すように、周波数間バイアス推定装置３２は、総電子数算出部３３、モデル値算出部３４、周波数間バイアス推定部３５、及び乱数発生部３６で構成されている。
【００２１】
総電子数算出部３３は、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。具体的には、総電子数算出部３３は、衛星信号受信機２４から測位衛星１０ａ、１０ｂ、１０ｃの各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。また、総電子数算出部３３は、外部インターネット網４６から得た受信機と衛星との擬似距離に関するデータに基づき衛星信号の通過経路の総電子数を算出することもできる。総電子数の算出方法については、後述する。
【００２２】
モデル値算出部３４は、上述したようなＩＲＩモデルやＧａｌｌａｇｈｅｒモデル等の電離層電子密度モデル関数に基づき総電子数モデル値を算出する。
【００２３】
乱数発生部３６は、乱数を発生させる。
【００２４】
周波数間バイアス推定部３５は、総電子数算出部３３により算出された総電子数とモデル値算出部３４により算出された総電子数モデル値に基づき周波数間バイアスを算出し、乱数発生部３６により発生した乱数に基づき周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定する。さらに、周波数間バイアス推定部３６は、予め精密に測量された位置に基づく位置情報と、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとを真の値であると推定する。
【００２５】
次に、上述のように構成された本実施の形態の作用を説明する。なお、電離層がダイナミックに変化するような時間帯であると推定がより複雑になるため、本実施例において、周波数間バイアス推定装置は、ローカル時間の深夜２時頃に測定したデータを用いるものとする。図３は、本実施例の周波数間バイアス推定装置３２の動作を示すフローチャート図である。まず、モデル値算出部３４は、電離層電子密度モデル関数に基づき総電子数モデル値を算出する（ステップＳ１０１）。ここで、電離層電子密度モデル関数をｆとすると、通過経路の総電子数モデル値は、通過経路の電子密度を積分することにより求められる。
【数１】
【００２６】
（１）式は、衛星Ｍｇから受信観測点Ｒまでに渡り、電離層電子密度関数ｆを経路に沿って積分した結果を示す。ここで、ｓは、位置を示すパラメータである。また、ｇは、衛星を識別する識別子である。ＴＥＣｇｔｒｕｅは、通過経路における真の総電子数を示す。Δｍｏｄｅｌは、電離層密度分布モデルに含まれる誤差を示す。
【００２７】
モデル値算出部３４は、算出した総電子数モデル値を周波数間バイアス推定部３５に出力する。
【００２８】
次に、総電子数算出部３３は、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する（ステップＳ１０３）。ここでは、測位衛星１０ａと衛星信号受信機２４との間の擬似距離について考える。衛星信号通過経路のＴＥＣは、測位衛星１０ａから送信される衛星信号（ここではＬ１とＬ２で示す）に基づき求められる。まず、擬似距離（コード距離、シュードレンジ）と位相距離（フェーズ距離）は、以下のように表すことができる。
【数２】
【００２９】
ここで、ρは、擬似距離を示す。また、Φは、位相距離を示す。ｒは、真の距離を示す。ｃは光速である。また、δｔｕは、受信機時刻誤差を示し、δｔｓは、衛星時刻誤差を示す。なお、本実施例において、右下の添え字は基本的に地上（受信機等）に関連する項を示し、右上の添え字は基本的に上空（衛星等）に関連する項を示す。δｔｕ、Ｌ１ｏｒＬ２ ｂｉａｓは、受信機周波数依存ハードウェア依存バイアスを示し、δｔｓＬ１ｏｒＬ２ ｂｉａｓは、衛星周波数依存ハードウェア依存バイアスを示す。さらに、Ｉは、電離層伝搬遅延量を示し、Ｔは、対流圏伝搬遅延量を示す。Ｎａｍｂは、整数不確定値を示し、εは、観測誤差を示す。
【００３０】
総電子数算出部３３は、２つの周波数の観測値の差をとることにより、ＴＥＣを求めることができる。
【数３】
【００３１】
ここで、ｆは周波数である。また、λは波長を示す。Δｂｉａｓは、周波数間バイアスを示す。また、Δ´ｂｉａｓは、位相についてのバイアスである。
【００３２】
（４）式より、ＴＥＣは、次式のように求められる。
【数４】
【００３３】
ここで、ＴＥＣｔｒｕｅは、真の総電子数を示す。（６）式に示すＴＥＣの推定値は、バイアスに関連した項を含む。このままでは誤差ε´の値が大きいため、総電子数算出部３３は、キャリアスムージングを実施して誤差ε´の値を小さくする。キャリアスムージングとは、ノイズを多く含む擬似距離に精密な搬送波差分を適用し、ノイズの低減を図るものである。キャリアスムージング後のＴＥＣを＜ＴＥＣ＞と表すものとすると、＜ＴＥＣ＞は、以下の式により求められる。
【数５】
【００３４】
ここで、ｍは、データ収集の時間順につけた番号である。また、Ｋは、スムージングの定数であり、適宜変更できるものとする。Ｋは、サンプリング時間間隔にも依存する。例えば、時定数を１８０秒にすると、Ｋ＝１８０／ｄｔとなる。ｄｔは、サンプリング時間間隔を表す。＜Δｂｉａｓ＞は、１ヶ月程度の期間一定であると仮定し、データ収集期間中一定であるものとする。したがって、＜Δｂｉａｓ＞は、添え字ｍを必要としない。
【００３５】
総電子数算出部３３は、算出した総電子数を周波数間バイアス推定部３５に出力する。
【００３６】
次に、周波数間バイアス推定部３５は、総電子数算出部３３により算出された総電子数とモデル値算出部３４により算出された総電子数モデル値に基づき周波数間バイアスを算出する。まず周波数間バイアス推定部３５は、衛星観測のデータに基づき（８）式を用いて総電子数算出部３３により求められた総電子数と、（１）式を用いてモデル値算出部３４により求められた総電子数モデル値との差を求める（ステップＳ１０５）。
【数６】
【００３７】
（１）式で用いたような既存の電離層電子密度モデルは、月平均では観測値とよい一致を示す。そのため、（９）式のΔｇ，ｍｏｄｅｌは、１ヶ月間のデータで平均することにより、ゼロに近い値となる。さらに、衛星や受信機の周波数間バイアスは、時間的変動が小さく、時定数を１ヶ月程度と考え、この間に収集したデータについては一定であるとする。１ヶ月間の平均操作を［］で表すとすると、（９）式に示す観測ＴＥＣとモデルＴＥＣとの差は、以下のようになる。
【数７】
【００３８】
ここで、ｄεは、平均後にも残る微小なモデル誤差を示す。本実施例において、ｄεは、ゼロであるものとする。衛星依存周波数間バイアスと受信機依存周波数間バイアスとの間に相関は無いと考えられる。したがって、観測ＴＥＣとモデルＴＥＣとの差は、（１０）式に示すように２つの独立した項にまとめることができる。Ｂｇは、衛星ｇに関するバイアスを示し、ＢＲは、受信機に関するバイアスを示す。受信機依存周波数間バイアスと衛星ｇのバイアスとが、（１０）式により求められる。したがって、受信機のバイアスが決まれば、各衛星のバイアスが決まる。
【００３９】
（１０）式は、１観測地点・観測時刻に対して、衛星数と同数できる。ここで、衛星数をＮｓａｔとすると、未知数は、観測場所受信機のバイアスを含む（Ｎｓａｔ＋１）個ある。したがって、単純に連立方程式により解を求めることはできない。
【００４０】
そこで、周波数間バイアス推定部３５は、（１０）式の受信機に関連した項を、乱数発生部３６により発生した乱数により仮に決定する（ステップＳ１０７）。さらに、周波数間バイアス推定部３５は、受信機依存周波数間バイアスを仮に決定したので、（１０）式より衛星に関するバイアス成分を決定することができる。
【００４１】
次に、決定したバイアス値に基づき測位が行われる（ステップＳ１０９）。測位を行う主体は、測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、及び衛星信号処理系２０からなる測位システムであり、特に衛星信号処理装置２６である。ただし、周波数間バイアス推定装置３２が測位を行う構成でもよい。測位場所の真値をＸ０とする。Ｘ０は、あらかじめ精密に測量されている。この真値Ｘ０に対して、仮に決定したバイアス値に基づいた測位の誤差をΔＸとする。衛星と受信機との距離の誤差であるΔＲは、以下の式のように表せる。
【数８】
【００４２】
ここで、（１３）式の＜ρ＞は、（８）式のように、キャリアスムージングを行った後の擬似距離を示す。（１１）式において、ｒＲｅｐｈは、精密に測量された場所Ｘ０と衛星の位置から計算された幾何学的距離を示す。この衛星の位置は、ＩＧＳより提供されるＧＰＳ位置データや衛星放送エフェメリスデータから計算されるＧＰＳ位置データである。また、Ｔｅｓｔｇは、対流圏遅延量推定値を示す。この対流圏遅延量を推定するに際し、対流圏遅延量モデルが使用されている。＜ＩＬ１，ｅｓｔｇ＞は、電離層遅延量推定値を示す。また、ΔｔＲ，ｃｌｏｃｋは、受信機のクロックオフセットを示す。衛星のクロックオフセット（Δｔｓ，ｃｌｏｃｋ）は、衛星放送エフェメリスデータにより補正されている。この電離層遅延量推定値＜ＩＬ１，ｅｓｔｇ＞は、バイアスを含む。受信機で生じるバイアスは、各衛星共通であるため、受信機のクロックオフセットと同じ影響を持ち、測位位置誤差に影響を与えない。
【００４３】
（１１）式に示す距離誤差ΔＲと測位誤差ΔＸとの関係は、以下のように示される。
【数９】
【００４４】
ここで、Ａは、精密に測量された場所Ｘ０から衛星までの方向ベクトルを要素とする行列である。（ｌ，ｍ，ｎ）が方向ベクトルを示す。添え字は、衛星番号を示す。したがって、（１７）式（１８）式において、衛星は、Ｎ個存在する。測位位置に関して、ΔＸは、最小二乗法等により求められる。（１９）式は、最小二乗法を用いた計算例である。ΔＸの４番目の項は、受信機クロックオフセットに関連した項である。この項は、（１１）式の２番目の括弧内の項に関係するため、（１１）式の２番目の括弧内の項が反映されたものとなっていて、測位位置に影響しない。（１１）式の１番目の項は、測位位置誤差に影響を与える。衛星依存周波数間バイアスの推定値が正しければ、ΔＸの位置に関する１〜３番目の項の大きさ（誤差）は、小さくなる。
【００４５】
（１０）式の受信機に関連した項を、乱数により決定することで、衛星依存周波数間バイアス値が仮に決定される。ここで、乱数により決定した受信機依存周波数間バイアスをＢＲ，ｒａｎｄとすると、仮に決定した衛星依存周波数間バイアス値Ｂｃａｎｄｇは、以下の式により求められる。
【数１０】
【００４６】
したがって、ＢＲ，ｒａｎｄと、（２０）式により求められたＢｃａｎｄｇとを（１１）式に代入することで、擬似距離誤差ΔＲは、以下のようになる。
【数１１】
【００４７】
ここで、（２１）式の１番目の括弧内にある衛星バイアス分が真の値に近い場合に、ΔＸの位置誤差に関する大きさは、最小になる。ΔＸの位置誤差Ｐｅｒｒは、例えば位置誤差のデカルト座標成分から以下のように表すことができる。
【数１２】
【００４８】
周波数間バイアス推定部３５は、乱数発生部３６により発生した乱数に基づき仮の受信機依存周波数間バイアスＢＲ，ｒａｎｄを決定する。このＢＲ，ｒａｎｄが変わる度に（２２）式の位置誤差Ｐｅｒｒは、変化する。Ｐｅｒｒが最小の値をとる場合に、周波数間バイアス推定部３５は、その際の受信機依存周波数間バイアスＢＲ，ｒａｎｄを真の値であると推定する（ステップＳ１１１）。さらに、真の値であると推定されたＢＲ，ｒａｎｄに基づき、衛星依存周波数間バイアスＢｅｓｔｇは、（２０）式を用いて算出される（ステップＳ１１３）。
【００４９】
乱数で推定値を最後まで求めることは、効率が悪い。そこで、周波数間バイアス推定部３５は、予め測量された位置に基づく位置情報と、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき位置を測位して得られた位置情報との誤差の変化率に基づき受信機依存周波数間バイアスと衛星依存周波数間バイアスとを推定する。具体的には、乱数発生部３６により発生した乱数に基づき、周波数間バイアス推定部３５は、受信機依存周波数間バイアスをｋ回生成（ＢＲ，１，・・・，ＢＲ，ｋ）し、測位位置誤差（Ｐｅｒｒ１，・・・，Ｐｅｒｒｋ）を求める。さらに、周波数間バイアス推定部３５は、受信機に関連した項の変化率に対する測位位置誤差Ｐｅｒｒの変化率を計算する。
【数１３】
【００５０】
（２３）式により、ＢＲに対する測位位置誤差Ｐｅｒｒの変化率ηｋが求められる。このηｋに基づき、受信機依存周波数間バイアスＢＲ，ｅｓｔは、以下のように求められる。
【００５１】
まず、ηｋの符号が、（ＢＲ，１，・・・，ＢＲ，ｋ）の範囲で全てプラス又はマイナスである場合、ＢＲ，ｅｓｔは、以下のようになる。
【数１４】
【００５２】
周波数間バイアス推定部３５は、（２４）式により求められたＢＲ，ｅｓｔを真の受信機依存周波数間バイアスとして推定する。さらに詳細に受信機のバイアス値を求める場合には、周波数間バイアス推定部３５は、（２４）式により求められたＢＲ，ｅｓｔの値の近辺で、再度乱数で受信機に関連した項（ＢＲ，１，・・・，ＢＲ，ｋ）を求め、（２３）式を計算する。
【００５３】
次に、ηｋの符号が、（ＢＲ，１，・・・，ＢＲ，ｋ）の範囲内でマイナスからプラスに変化する場合、位置誤差Ｐｅｒｒの最小値は、符号が変わった近辺に存在する。そこで、符号がかわる 前後のＢＲ，ｋを、ＢＲ，ｋ−、ＢＲ，ｋ＋とすると、求める受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔは、以下のように求められる。
【数１５】
【００５４】
ここで、ηｋ＋、ηｋ−は、それぞれ、プラス符号を持つＢＲ，ｋ＋における変化率とマイナス符号を持つＢＲ，ｋ−の変化率を示す。
【００５５】
また、単純に中間値を推定値とする方法もある。その場合、受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔは、以下のように求められる。
【数１６】
【００５６】
上述したような方法により、受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔが求められると、（１０）式の関係により、各衛星の衛星依存周波数間バイアスＢＲ，ｅｓｔｇは、以下のようにして求められる。
【数１７】
【００５７】
以上のように求められた周波数間バイアスは、特定の受信機Ｒにおいて計算した受信機依存周波数間バイアスの推定値ＢＲ，ｅｓｔと各衛星の衛星依存周波数間バイアスＢＲ，ｅｓｔｇである。したがって、この推定結果は、観測地点の地形によるマルチパスの影響などの影響を含む可能性がある。
【００５８】
そこで、その影響を小さくするため、周波数間バイアス推定部３５は、測位システムが受信機（図１における衛星信号処理系２０）を複数有する場合に、複数の受信機の各々で得られた位置情報に基づき１つの測位衛星に対して得られた複数の衛星依存周波数間バイアスを平均化することにより衛星依存周波数間バイアスを推定する。具体的には、衛星依存周波数間バイアスの平均値Ｂｅｓｔｇは、以下の式により求められる。
【数１８】
【００５９】
ここで、Wは処理を行った受信機の数を示している。以上のようにして、周波数間バイアス推定装置３２は、周波数間バイアス（受信機依存周波数間バイアス及び衛星依存周波数間バイアス）を求めることができる。また、以上述べた動作は、周波数間バイアス推定方法に適用可能である。
【００６０】
上述のとおり、本発明の実施例１の形態に係る周波数間バイアス推定装置によれば、乱数に基づき受信機依存周波数間バイアス及び衛星依存周波数間バイアスを決定するので、単純に連立方程式により求めることのできない周波数間バイアスを推定することができる。
【００６１】
また、予め精密に測量された位置に基づく位置情報と、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとを真の値であると推定するので、衛星に起因する周波数間バイアス及び受信機に起因する周波数間バイアスを精度よく推定することができる。
【００６２】
さらに、精度良く推定された受信機依存周波数間バイアスと衛星依存周波数間バイアスとに基づき、電離層総電子数推定精度を向上し、ＧＰＳ等の測位システムにおける測位精度や電離層モデルの改善に寄与することができる。また、ＨＦ帯等の周波数の伝搬経路を算出する際に使用する電離層モデルの精度を向上させることができ、より有効な通信を可能にすることができる。また、衛星を使った航法の向上が可能である。
【産業上の利用可能性】
【００６３】
本発明に係る周波数間バイアス推定装置は、測位衛星から送信される衛星信号を利用して測位を行うＧＰＳ等の測位システムに利用可能である。
【図面の簡単な説明】
【００６４】
【図１】本発明の実施例１の形態の周波数間バイアス推定装置及び測位システムの構成を示すブロック図である。
【図２】本発明の実施例１の形態の周波数間バイアス推定装置の詳細なブロック図である。
【図３】本発明の実施例１の形態の周波数間バイアス推定装置の動作を示すフローチャート図である。
【符号の説明】
【００６５】
１０ａ、１０ｂ、１０ｃ 測位衛星
２０ 衛星信号処理系
２２ アンテナ
２４ 衛星信号受信機
２６ 衛星信号処理装置
３０ 電離層ＴＥＣ・バイアス推定処理系
３２ 周波数間バイアス推定装置
３３ 総電子数算出部
３４ モデル値算出部
３５ 周波数間バイアス推定部
３６ 乱数発生部
４０ インターネットデータ処理系
４２ ルータ
４４ ＧＥＯＮＥＴ収集データ処理装置
４６ 外部インターネット網
５０ 通信回線
【特許請求の範囲】
【請求項１】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する周波数間バイアス推定装置であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、
電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出部と、
乱数を発生させる乱数発生部と、
前記総電子数算出部により算出された総電子数と前記モデル値算出部により算出された総電子数モデル値に基づき周波数間バイアスを算出し、前記乱数発生部により発生した乱数に基づき前記周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定部と、
を備えることを特徴とする周波数間バイアス推定装置。
【請求項２】
前記周波数間バイアス推定部は、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差の変化率に基づき前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを推定することを特徴とする請求項１記載の周波数間バイアス推定装置。
【請求項３】
前記周波数間バイアス推定部は、前記測位システムが前記受信機を複数有する場合に、複数の前記受信機の各々で得られた位置情報に基づき１つの前記測位衛星に対して得られた複数の前記衛星依存周波数間バイアスを平均化することにより前記衛星依存周波数間バイアスを推定することを特徴とする請求項１又は請求項２記載の周波数間バイアス推定装置。
【請求項４】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する周波数間バイアス推定方法であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出ステップと、
電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出ステップと、
乱数を発生させる乱数発生ステップと、
前記総電子数算出ステップにより算出された総電子数と前記モデル値算出ステップにより算出された総電子数モデル値に基づき周波数間バイアスを算出し、前記乱数発生ステップにより発生した乱数に基づき前記周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定ステップと、
を備えることを特徴とする周波数間バイアス推定方法。
【請求項１】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する周波数間バイアス推定装置であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、
電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出部と、
乱数を発生させる乱数発生部と、
前記総電子数算出部により算出された総電子数と前記モデル値算出部により算出された総電子数モデル値に基づき周波数間バイアスを算出し、前記乱数発生部により発生した乱数に基づき前記周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定部と、
を備えることを特徴とする周波数間バイアス推定装置。
【請求項２】
前記周波数間バイアス推定部は、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差の変化率に基づき前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを推定することを特徴とする請求項１記載の周波数間バイアス推定装置。
【請求項３】
前記周波数間バイアス推定部は、前記測位システムが前記受信機を複数有する場合に、複数の前記受信機の各々で得られた位置情報に基づき１つの前記測位衛星に対して得られた複数の前記衛星依存周波数間バイアスを平均化することにより前記衛星依存周波数間バイアスを推定することを特徴とする請求項１又は請求項２記載の周波数間バイアス推定装置。
【請求項４】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する周波数間バイアス推定方法であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出ステップと、
電離層電子密度モデル関数に基づき総電子数モデル値を算出するモデル値算出ステップと、
乱数を発生させる乱数発生ステップと、
前記総電子数算出ステップにより算出された総電子数と前記モデル値算出ステップにより算出された総電子数モデル値に基づき周波数間バイアスを算出し、前記乱数発生ステップにより発生した乱数に基づき前記周波数間バイアスが有する受信機依存周波数間バイアスと衛星依存周波数間バイアスとを仮に決定し、予め測量された位置に基づく位置情報と、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとに基づき前記位置を測位して得られた位置情報との誤差が最も小さい場合に、仮に決定された前記受信機依存周波数間バイアスと前記衛星依存周波数間バイアスとを真の値であると推定する周波数間バイアス推定ステップと、
を備えることを特徴とする周波数間バイアス推定方法。
【図１】
【図２】
【図３】
【図２】
【図３】
【公開番号】特開２００８−２４９４０２（Ｐ２００８−２４９４０２Ａ）
【公開日】平成２０年１０月１６日（２００８．１０．１６）
【国際特許分類】
【出願番号】特願２００７−８８９１２（Ｐ２００７−８８９１２）
【出願日】平成１９年３月２９日（２００７．３．２９）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【Ｆターム（参考）】
【公開日】平成２０年１０月１６日（２００８．１０．１６）
【国際特許分類】
【出願日】平成１９年３月２９日（２００７．３．２９）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【Ｆターム（参考）】
[ Back to top ]