電子密度推定装置及び電子密度推定方法
【課題】３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も可能とする電子密度推定装置を提供する。
【解決手段】複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定装置であって、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき衛星信号の通過経路の総電子数を算出する総電子数算出部と、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する３次元電子密度推定空間設定部３６と、総電子数算出部により算出された総電子数に基づいて３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する暫定電子密度推定部３８と、推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する３次元電子密度推定部３９とを備える。
【解決手段】複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定装置であって、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき衛星信号の通過経路の総電子数を算出する総電子数算出部と、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する３次元電子密度推定空間設定部３６と、総電子数算出部により算出された総電子数に基づいて３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する暫定電子密度推定部３８と、推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する３次元電子密度推定部３９とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、衛星から送信される複数の異なる周波数信号による通過経路の電離層電子密度分布を推定する電子密度推定装置及び電子密度推定方法に関する。
【背景技術】
【０００２】
衛星ナビゲーションシステム「Ｇａｌｉｌｅｏ」やＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）、ＱＺＳＳ（Ｑｕａｓｉ Ｚｅｎｉｔｈ Ｓａｔｅｌｌｉｔｅｓ Ｓｙｓｔｅｍ；準天頂衛星システム）等の衛星から送信される信号は、地球上の観測点の測位を行うのに利用される。これらの衛星により放送される複数周波数信号の擬似距離差は、電波の伝搬経路上の電子数に関係していることが知られている。
【０００３】
衛星から地上に向けて送信される信号は、電離層を通過する際に伝搬遅延量を生じる。この電離層伝搬遅延量は、伝搬する信号の周波数に依存している。そのため、衛星からの複数の周波数信号を受信することにより、信号の通過経路における総電子数（ＴＥＣ：Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ）を得ることができる。電波の伝搬経路上の総電子数を推定する場合には、複数周波数信号の送信・受信時刻による擬似距離（コード距離）差を使用する方法や複数周波数のキャリア位相擬似距離（フェーズ距離）の差を利用する方法が知られている。
【０００４】
非特許文献１や非特許文献２には、電離層における電子密度モデル関数であるＩＲＩ（Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ）モデルが記載されている。さらに、非特許文献４には、やはり電離層電子密度モデル関数の一つであるＧａｌｌａｇｈｅｒのモデルが記載されている。これらのモデル関数を用いて計算によりＴＥＣを求めることもできる。
【０００５】
非特許文献５や非特許文献６には、トモグラフィ手法による衛星ＴＥＣを用いた３次元電子密度推定が記載されており、ＭＡＲＴ（Ｍｕｌｔｉｐｌｉｃａｔｉｖｅ Ａｌｇｅｂｒａｉｃ Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ Ｔｅｃｈｎｉｑｕｅ）と呼ばれる代数的再構成法が適用されている。一方、非特許文献７には、ＧＳＶＤ（Ｇｅｎｅｒａｌｉｚｅｄ Ｓｉｎｇｕｌａｒ Ｖａｌｕｅ Ｄｅｃｏｍｐｏｓｉｔｉｏｎ）を用いたトモグラフィック再構成による電子密度推定が記載されている。
【０００６】
また、非特許文献８及び非特許文献９には、電離層を球核の薄い層とみなした電子密度計算について記載されている。さらに、非特許文献１０には、３次元方向の電子密度分布を高さ方向に計算した場合について記載されている。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ ２０００”，Ｒａｄｉｏ．Ｓｃｉｅｎｃｅ，Ｖｏｌ．３６，Ｎｕｍｂｅｒ２，ＰＰ２６１−２７５，Ｍａｒｃｈ／Ａｐｒｉｌ，２００１
【非特許文献２】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ，ｅｔ，ａｌ．，：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ １９９０”，Ｎｏｖｅｍｂｅｒ，１９９０．
【非特許文献３】Ａ．Ｋｏｍｊａｔｈｙ：‘Ｇｌｏｂａｌ Ｉｏｎｏｓｐｈｅｒｉｃ Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ Ｍａｐｐｉｎｇ Ｕｓｉｎｇ ｔｈｅ Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ’ ，ＵＮＢ，Ｔｅｃｈｎｉｃａｌ Ｒｅｐｏｒｔ Ｎｏ．１８８，Ｓｅｐ．１９９７．
【非特許文献４】Ｇａｌｌａｇｈｅｒ，Ｄ．Ｌ．，Ｐ．Ｄ．Ｃｒａｖｅｎ，ａｎｄ Ｒ．Ｈ．Ｃｏｍｆｏｒｔ，Ｇｌｏｂａｌ ｃｏｒｅ ｐｌａｓｍａ ｍｏｄｅｌ，Ｊ．Ｇｅｏｐｈｙｓ．Ｒｅｓ．１０５，Ａ８，１８，８１９−１８，８３３，２０００．
【非特許文献５】Ｒａｙｍｕｎｄ，Ｔ．Ｄ．，Ａｕｓｔｉｎ，Ｊ．Ｒ．，Ｆｒａｎｋｅ，Ｓ．Ｊ．，Ｌｉｕ，Ｃ．Ｈ．，Ｋｌｏｂｕｃｈａｒ，Ｊ．Ａ．，ａｎｄ Ｓｔａｌｋｅｒ，Ｊ．，：“Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｃｏｍｕｔｅｒｉｚｅｄ Ｔｏｍｏｇｒａｐｈｙ ｔｏ ｔｈｅ Ｉｎｖｅｓｔｉｇａｔｉｏｎ ｏｆ Ｉｏｎｏｓｐｈｅｒｉｃ Ｓｔｒｕｃｔｕｒｅｓ”，Ｒａｄｉｏ Ｓｃｉ．，２５（５），７７１−７８９，１９９０．
【非特許文献６】Ｓｔｅｆａｎ Ｓｃｈｌｕｅｒ，ｅｔ ａｌ．：“Ｍｏｎｉｔｏｒｉｎｇ ｔｈｅ ３ Ｄｉｍｅｎｓｉｏｎａｌ Ｉｏｎｏｓｐｈｅｒｉｃ Ｅｌｅｃｔｒｏｎ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｂａｓｅｄ ｏｎ ＧＰＳ Ｍｅａｓｕｒｅｍｅｎｔｓ”
【非特許文献７】Ｋ．Ｂｈｕｙａｎ，ｅｔ ａｌ：“Ｔｏｍｏｇｒａｐｈｉｃｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｏｆ ｔｈｅ Ｉｏｎｏｓｐｈｅｒｅ ｕｓｉｎｇ ｇｅｎｅｒａｌｉｚｅｄ ｓｉｎｇｕｌａｒ ｖａｌｕｅ ｄｅｃｏｍｐｏｓｉｔｉｏｎ” ，Ｃｕｒｒｅｎｔ Ｓｃｉ．，Ｖｏｌ１．８３，Ｎｏ．９，１０ Ｎｏｖ．，２００２．
【非特許文献８】Ｏｔｕｋａ，Ｙ．，ｅｔ．ａｌ．：“Ａ ｎｅｗ ｔｅｃｈｎｉｑｕｅ ｆｏｒ ｍａｐｐｉｎｇ ｏｆ ｔｏｔａｌ ｅｌｅｃｔｒｏｎ ｃｏｎｔｅｎｔ ｕｓｉｎｇ ＧＰＳ ｎｅｔｗｏｒｋ ｉｎ Ｊａｐａｎ，”，Ｅａｒｔｈ Ｐｌａｎｅｔｓ Ｓｐａｃｅ，１１１−１２０，２００１．
【非特許文献９】Ｍａ，Ｇ．，ａｎｄ Ｔ．Ｍａｒｕｙａｍａ（２００３），Ｄｅｒｉｖａｔｉｏｎ ｏｆ ＴＥＣ ａｎｄ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｉｎｓｔｒｕｍｅｎｔａｌ ｂｉａｓｅｓ ｆｒｏｍ ＧＥＯＮＥＴ ｉｎ Ｊａｐａｎ，Ａｎｎ．Ｇｅｏｐｈｙｓ．，２１，２０８３−２０９３．
【非特許文献１０】五十嵐来良，斉藤昭則，大塚雄一：“標準電離層モデル（ＩＲＩ）で導出したＴＥＣとＧＰＳ受信機網で観測したＴＥＣの比較解析”，第１回電離圏の利用と影響に関するシンポジウム講和集，ｐｐ２４−１〜８，２００３．
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、非特許文献５や非特許文献６に記載された手法は、観測した衛星の数が限定されるため、限られた空間の領域にしか適用できないという短所がある。特に、非特許文献５に記載された手法は、空間平均を使用しており、電離層が有する特異な高さ方向の電子密度の谷（Ｅ層とＦ層の間にある密度が小さくなる部分）を消してしまうため、電離層の高さ方向変化を捉えることができない。また、非特許文献６には、空間平均についての記載がされていない。
【０００９】
非特許文献８及び非特許文献９に記載された手法は、上述したように電離層を球核の薄い層とみなして電子密度を計算しているため、３次元高さ方向の分布を知ることができない。また、非特許文献１０は、３次元方向の電子密度分布を高さ方向に計算しているが、水平方向の広がりについての言及が無く、緯度・経度方向の広がりがある空間の電子密度を求めていない。これらにより、従来の文献においては、ＧＰＳ等の衛星で観測できない空間部分の取り扱い方法が不十分である。
【００１０】
すなわち、従来の技術は、３次元の電離層電子密度分布を十分推定しきれていない。従来技術における手法は、空間平均を実施して電子密度を推定するため、本来なら通過する電波であっても反射したと判断してしまう場合がある。また、推定された空間とそれ以外の空間との差が大きい場合に、空間的な電子密度分布変化に不連続が生じてしまう。このように、従来手法は、電子密度変化に不連続が生じた結果、電子密度変化率を使用して電波伝搬経路を推定する際に大きな誤差を生じてしまうという問題がある。
【００１１】
本発明は上述した従来技術の問題点を解決するもので、３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も可能とする電子密度推定装置及び電子密度推定方法を提供することを課題とする。
【課題を解決するための手段】
【００１２】
本発明に係る電子密度推定装置は、上記課題を解決するために、複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定装置であって、前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定部と、前記総電子数算出部により算出された総電子数に基づいて前記空間設定部により分割された領域毎の電子密度を推定する第１電子密度推定部と、前記第１電子密度推定部により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定部とを備えることを特徴とする。
【００１３】
本発明に係る電子密度推定方法は、上記課題を解決するために、複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定方法であって、前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出ステップと、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定ステップと、前記総電子数算出ステップにより算出された総電子数に基づいて前記空間設定ステップにより分割された領域毎の電子密度を推定する第１電子密度推定ステップと、前記第１電子密度推定ステップにより推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定ステップとを備えることを特徴とする。
【発明の効果】
【００１４】
本発明によれば、３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も可能とする電子密度推定装置及び電子密度推定方法を提供することができる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施例１の形態の電子密度推定装置及び測位システムの構成を示すブロック図である。
【図２】本発明の実施例１の形態の電子密度推定装置の詳細なブロック図である。
【図３】本発明の実施例１の形態の電子密度推定装置の動作を示すフローチャート図である。
【図４】本発明の実施例１の形態の電子密度推定装置内の３次元電子密度推定空間設定部による３次元空間設定を説明する図である。
【図５】本発明の実施例１の形態の電子密度推定装置内の３次元電子密度推定空間設定部による補正係数設定を説明する図である。
【発明を実施するための形態】
【００１６】
以下、本発明の電子密度推定装置及び電子密度推定方法の実施の形態について図面を参照して詳細に説明する。
【実施例１】
【００１７】
以下、本発明の実施例について図面を参照しながら説明する。図１は、本発明の実施例１の電子密度推定装置３２及び測位システムの構成を示すブロック図である。この電子密度推定装置３２は、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される衛星信号を受信する１以上の衛星信号受信機２４を有し且つ衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する。
【００１８】
まず、本実施の形態の構成を説明する。本実施例の電子密度推定装置３２及び測位システムは、図１に示すように、測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系３０、及びインターネットデータ処理系４０で構成されている。測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、及び衛星信号処理系２０は、本発明の測位システムに対応する。なお、インターネットデータ処理系４０も含めて測位システムとすることも可能であるが、インターネットデータ処理系４０は、必ずしも必須のものではなく、付加的なものである。また、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定処理系３０、及びインターネットデータ処理系４０は、お互いに通信回線５０で接続されている。
【００１９】
測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃは、ＧＰＳ、Ｇａｌｉｌｅｏ、準天頂衛星（ＱＺＳＳ）等の航法衛星であり、複数の周波数の衛星信号を送信する。
【００２０】
衛星信号処理系２０は、アンテナ２２、衛星信号受信機２４、及び衛星信号処理装置２６から構成される。衛星信号受信機２４は、本発明の受信機に対応し、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される複数の周波数の衛星信号をアンテナ２２を介して受信する。また、衛星信号処理装置２６は、衛星信号に含まれる測位情報を用いて位置情報を得る。
【００２１】
電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系３０は、電子密度推定装置３２を有し、受信した信号に基づき衛星位置、信号通過経路の電離層総電子数（ＴＥＣ）、周波数間バイアス等を算出する。受信したデータ・処理結果は、データサーバ系へＬＡＮ経由で伝送され保存される。
【００２２】
インターネットデータ処理系４０は、ルータ４２、ＧＥＯＮＥＴ収集データ処理装置４４、及び外部インターネット網４６で構成されている。ルータ４２は、スイッチングハブでもよい。インターネットデータ処理系４０は、公開されている電離層関連の情報や国土地理院が公開しているＧＰＳ観測データ（ＧＥＯＮＥＴデータ）や国際的に観測結果を公開しているＩＧＳ（Ｉｎｔｅｒｎａｔｉｏｎａｌ ＧＰＳ Ｓｅｒｖｉｃｅ ｆｏｒ Ｇｅｏｄｙｎａｍｉｃｓ）データ等をインターネット経由で収集する装置である。インターネットデータ処理系４０により収集された結果は、電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系３０で処理される。ルータ４２は、セキュリティを考慮して設けられ、ファイアウォールとする。
【００２３】
図２は、本発明の実施例１の電子密度推定装置３２の詳細なブロック図である。図２に示すように、電子密度推定装置３２は、衛星ＴＥＣ算出部３３、ＴＥＣバイアス補正部３４、周波数間バイアス推定部３５、３次元電子密度推定空間設定部３６、電子密度モデル算出部３７、暫定電子密度推定部３８、及び３次元電子密度推定部３９で構成されている。
【００２４】
衛星ＴＥＣ算出部３３、ＴＥＣバイアス補正部３４、及び周波数間バイアス推定部３５は、本発明の総電子数算出部に対応し、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。具体的には、衛星ＴＥＣ算出部３３は、衛星信号受信機２４から測位衛星１０ａ、１０ｂ、１０ｃの各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。
【００２５】
周波数間バイアス推定部３５は、周波数間バイアスを算出し、ＴＥＣバイアス補正部３４に出力する。周波数間バイアス推定部３５による周波数間バイアスの推定方法は、どのようなものでもよく、従来手法を用いてよい。
【００２６】
ＴＥＣバイアス補正部３４は、周波数間バイアス推定部３５により推定された周波数間バイアスを用いて補正を行い、最終的な総電子数を算出する。
【００２７】
３次元電子密度推定空間設定部３６は、本発明の空間設定部に対応し、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する。
【００２８】
電子密度モデル算出部３７は、上述したようなＩＲＩモデルやＧａｌｌａｇｈｅｒモデル等の電離層電子密度モデル関数に基づき電子密度モデルを算出する。
【００２９】
暫定電子密度推定部３８は、本発明の第１電子密度推定部に対応し、総電子数算出部（ＴＥＣバイアス補正部３４）により算出された総電子数に基づいて３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する。本実施例において、暫定電子密度推定部３８は、非特許文献５や非特許文献に記載されているようなＭＡＲＴの手法を用いて電子密度を推定するものとする。
【００３０】
３次元電子密度推定部３９は、本発明の第２電子密度推定部に対応し、暫定電子密度推定部３８により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する。
【００３１】
次に、上述のように構成された本実施の形態の作用を説明する。図３は、本実施例の電子密度推定装置３２の動作を示すフローチャート図である。まず、衛星ＴＥＣ算出部３３は、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき衛星信号の通過経路の総電子数を算出する（ステップＳ１）。具体的には、衛星ＴＥＣ算出部３３は、広範囲に設置した衛星受信機収集データからＧＰＳ等の衛星のＴＥＣを計算する。ここで、本実施例の衛星ＴＥＣ算出部３３は、衛星から放送される２周波の信号（ＧＰＳ衛星のＬ１周波数（１５７５．４２ＭＨｚ）とＬ２周波数（１２２７．６０ＭＨｚ）の信号）に基づいて信号伝搬経路のＴＥＣを算出するものとする。ＧＰＳ以外でも異なる２周波を放送している場合には同様な手法が使用できる。
【００３２】
ここでは、測位衛星１０ａと衛星信号受信機２４との間の擬似距離について考える。衛星信号通過経路のＴＥＣは、測位衛星１０ａから送信される衛星信号（ここではＬ１とＬ２で示す）に基づき求められる。まず、擬似距離（コード距離、シュードレンジ）とキャリア位相擬似距離（フェーズ距離）は、以下のように表すことができる。
【数１】
【００３３】
ここで、ρは、遅延時間による距離（擬似距離）を示す。また、Φは、キャリア位相距離を示す。ｒは、真の距離を示す。ｃは光速である。また、δｔｕは、受信機時刻誤差を示し、δｔｓは、衛星時刻誤差を示す。なお、本実施例において、右下の添え字は基本的に地上（受信機等）に関連する項を示し、右上の添え字は基本的に上空（衛星等）に関連する項を示す。δｔｕ、Ｌ１ｏｒＬ２ ｂｉａｓは、受信機周波数依存ハードウェア依存バイアスを示し、δｔｓＬ１ｏｒＬ２ ｂｉａｓは、衛星周波数依存ハードウェア依存バイアスを示す。さらに、Ｉは、電離層伝搬遅延量を示し、Ｔは、対流圏伝搬遅延量を示す。ＮＬ１ｏｒＬ２，ａｍｂは、整数不確定値を示し、εは、観測誤差を示す。
【００３４】
衛星ＴＥＣ算出部３３は、２つの周波数の観測値の差（遅延時間の差）をとることにより、ＴＥＣを求めることができる。
【数２】
【００３５】
ここで、ｆは周波数である。また、λは波長を示す。Δｂｉａｓは、周波数間バイアスを示す。また、ＴＥＣｔｒｕｅは、真の総電子数を示す。
【００３６】
次に、ＴＥＣバイアス補正部３４は、周波数間バイアス推定部３５により事前に推定された周波数間バイアスΔｂｉａｓによる補正をＴＥＣに実施する（ステップＳ２）。まず、バイアス補正後の２つの周波数の遅延時間差は、次式により表される。
【数３】
【００３７】
ここで、（６）式には、観測値から分からない不確定値（λＬ１ΔＮＬ１，ａｍｂ−λＬ２ΔＮＬ２，ａｍｂ）が含まれている。そこで、ＴＥＣバイアス補正部３４は、（５）式と（６）式とを組み合わせて、以下に示す式により整数不確定値を削除する。この手法は、非特許文献３にも記載されている。
【数４】
【００３８】
ここで、添え字ｋは、データの番号を示す。また、Ｍは、連続的に収集できたサンプル数の合計を示す。さらに、ＴＥＣｋ，ｕは、衛星ｕの伝搬経路上の総電子数（個／ｍ２）を示す。
【００３９】
電子密度モデル算出部３７は、電離層電子密度モデル関数に基づき電子密度モデルを算出し、算出した電子密度モデルを３次元電子密度推定空間設定部３６に出力する。なお、電子密度モデル算出部３７による電離層モデルは、例えば非特許文献１，２，４等に記載されている。また、ステップＳ１，Ｓ２（特にＳ１）は、本発明の総電子数算出ステップに対応する。
【００４０】
次に、３次元電子密度推定空間設定部３６は、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する。この３次元電子密度推定空間設定部３６の動作は、本発明の空間設定ステップに対応する。また、３次元電子密度推定空間設定部３６は、推定する３次元空間の電離層電子密度モデル値あるいは推定値を設定する（ステップＳ３）。
【００４１】
具体的には、３次元電子密度推定空間設定部３６は、衛星から衛星受信機が配置されている地点まで３次元空間における電離層を、緯度・経度・高さ方向に分割し、電離層モデル値を使用して分割した領域内あるいは格子点の位置の電子密度を設定する。非特許文献１，２等に記載されたＩＲＩモデルは、高さ１０００ｋｍまでが有効範囲であるため、衛星が存在する空間までは設定できない。そこで、３次元電子密度推定空間設定部３６は、例えば非特許文献４に記載されたＧＣＰＭ（Ｇｌｏｂａｌ Ｃｏｒｅ Ｐｌａｓｍａ Ｍｏｄｅｌ）を用いる手法や、ＩＲＩモデルの約４００ｋｍ〜１０００ｋｍを外挿する手法を採用することができる。
【００４２】
図４は、本実施例の電子密度推定装置３２内の３次元電子密度推定空間設定部３６による３次元空間設定を説明する図である。本実施例の３次元電子密度推定空間設定部３６は、図４に示すように、３次元空間をＸ１〜Ｘ１２の１２個の領域（ブロック）に分割しているが、実際には膨大な数のブロック数となる。図４において、測位衛星Ｓ１〜Ｓ４と衛星信号受信機Ｒｘ１〜Ｒｘ３との間における信号通過経路は、矢印により表されている。図４の衛星・ＴＥＣ・電子密度の関係は、次式により表される。
【数５】
【００４３】
（９）式において、Ｘ１〜Ｘ１２は、空間を分割した領域（メッシュ）の電子密度を示し、現時点では不明である。また、ＴＥＣｉｊは、ｉ衛星・ｊ受信機間で決まる総電子数を示す。例えば、図４中の太い矢印ＴＥＣ２１は、測位衛星Ｓ２と衛星信号受信機Ｒｘ１との間の信号通過経路における総電子数を示す。ＴＥＣバイアス補正部３４により出力されたＴＥＣの値は、（９）式の右辺に代入される。
【００４４】
ここで、（９）式の行列式をＡＸ＝Ｅと表すことにする。行列Ａの中のゼロでない要素（１で表されている要素）は、ｉ衛星・ｊ受信機間の信号が通過したメッシュを示す。
【００４５】
実際にモデルにより設定した電子密度の値は、メッシュ内の１点の値である。また、電波信号が各メッシュを通過する際に、電波はメッシュを縦に通過する場合のみならず斜めに通過する場合もある。この効果を考慮し、３次元電子密度推定空間設定部３６は、行列Ａの１の代わりに次式で示される値を代入する。
【数６】
【００４６】
ここで、ｄＨｉ，ｊは、ｉ，ｊメッシュの分割した高さ方向の大きさを示す。また、変換係数ＳＦは、斜めから入射した信号を高さ方向に入射した信号に補正する係数である。θは、衛星受信機が受信した衛星電波の入射仰角を示す。図５は、本実施例の電子密度推定装置３２内の３次元電子密度推定空間設定部３６による補正係数ＳＦ設定を説明する図である。
【数７】
【００４７】
（１２）式の左辺は、ｉ，ｊメッシュ内の垂直方向のＴＥＣを示す。また、右辺のＴＥＣｓは、ｉ，ｊメッシュ内を通過する電波方向のＴＥＣを示す。（１３）式内のＲは、図５に示すように、地球半径を示す。また、（１３）式内のｈは、地表面からｉ，ｊメッシュまでの高さを示す。
【００４８】
次に、暫定電子密度推定部３８は、ＭＡＲＴあるいはＧＳＶＤ手法により、衛星信号通過経路の電子密度を推定する（ステップＳ４）。具体的には、暫定電子密度推定部３８は、総電子数算出部（ＴＥＣバイアス補正部３４）により算出された総電子数に基づいて３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する。この暫定電子密度推定部３８の動作は、本発明の第１電子密度推定ステップに対応する。
【００４９】
なお、本実施例において、暫定電子密度推定部３８は、ＭＡＲＴあるいはＧＳＶＤ手法により電子密度を推定するが、本発明を実現するためには、必ずしもいずれかの手法を使う必要はなく、衛星信号受信機が得た観測値に基づいて暫定的に各領域の電子密度を推定できればよい。
【００５０】
ＭＡＲＴの手法は、上述したように非特許文献５，６等に記載されているが、ここで簡潔に説明する。
【数８】
【００５１】
（１４）式は、ＭＡＲＴのアルゴリズムを示す式である。ａｉは、（１０）式のＡ行列ｉ行ベクトルを表す。また、ａｉｊは、Ａ行列のｉｊ成分を表す。Ｘｊは、Ｘ列行列のｊ成分を表す。さらに、Ｅｉは、（９）式のＥ行列のｉ成分を表している。また、λは、収束の度合いを決めるパラメータである。
【００５２】
暫定電子密度推定部３８は、（１４）式に示すＭＡＲＴのアルゴリズムを使用して、（１４）式の値が収束するまで繰り返し処理を実施し、３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する。その際に、暫定電子密度推定部３８は、ＴＥＣバイアス補正部３４により出力されたＴＥＣの値をＥ行列の値に使用するとともに、電子密度モデル算出部３７により算出されたモデル値をＸの初期値として使用する。
【００５３】
なお、暫定電子密度推定部３８は、信号が全く通過していない空間の電子密度を推定することはできない。
【００５４】
ＧＳＶＤ手法を採用する場合には、暫定電子密度推定部３８は、以下に示す方式で電子密度を求める。なお、ＧＳＶＤ手法については、非特許文献７に記載がある。
【数９】
【００５５】
（１５）式のＬは、Ｔｉｋｈｏｎｏｖ行列と呼ばれる。０次のＴｉｋｈｏｎｏｖ正規化では、Ｌを単位行列とする。αは、正規化パラメータで、α＞０とする。また、Ｌは、正値２次形式行列（固有値が全て正）である。
【００５６】
（１５）式を満足する電子密度Ｘは、次式で表すことができる。
【数１０】
【００５７】
次に、３次元電子密度推定部３９は、高さ毎に球面調和関数により電子密度を近似する（ステップＳ５）。すなわち、３次元電子密度推定部３９は、暫定電子密度推定部３８により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する。この３次元電子密度推定部３９の動作は、本発明の第２電子密度推定ステップに対応する。
【００５８】
なお、３次元電子密度推定部３９は、その作用効果において本発明の特徴を示すものであり、従来手法と異なる点である。
【００５９】
（９）式のＸは、空間のメッシュを一次元で表現したものとなっている。３次元電子密度推定部３９は、暫定電子密度推定部３８により推定された電子密度の中から、同じ高さの電子密度だけを抽出する。
【数１１】
【００６０】
（１７）式において、ｘ＝ｓｉｎ（θ）である。Ｖは、ある高さにおける電子密度を示す。ただし、θはメッシュの緯度を示し、φはメッシュの経度を示す。また、Ｓｎ，Ｓｎｍは、Ｓｃｈｕｍｉｄｔ Ｓｅｍｉｎｏｒｍａｌａｉｚｅｄルジャンドル陪関数である。（２０）式にルジャンドル関数を示す。
【数１２】
【００６１】
Ｎｐは、ルジャンドル関数の次数を示す。３次元電子密度推定部３９は、同じ高さの電子密度を（１７）式を使用して近似する。これにより、ＭＡＲＴ手法で求められなかったメッシュの電子密度の値を求めることができる。またＧＳＶＤで求めた値を修正することができる。
【００６２】
すなわち、３次元電子密度推定部３９は、ある緯度・経度（θ，φ）について、高さ方向のメッシュ分（１７）式で電子密度を近似し、高さＨに対する係数（α（Ｈ）、β（Ｈ））を求める。この係数（α（Ｈ）、β（Ｈ））が求まることにより、３次元電子密度推定部３９は、緯度・経度・高さを指定し、（１７）式を用いて、その高さＨにおける電子密度Ｖ（ｘ，φ）を算出することができる。
【００６３】
なお、非特許文献８や非特許文献９に記載された従来手法は、ある高さのみに着目しており、高さ毎の電子密度分布を求めていない。また、非特許文献１０に記載された手法は、高さ方向に対する電子密度分布を求めているが、緯度・経度方向の広がりがある空間の電子密度を求めていない。
【００６４】
電子密度モデル算出部３７による電子密度モデル値と暫定電子密度推定部３８による電子密度の値とが大きくかけ離れている場合に、３次元電子密度推定部３９による近似は、適切でない場合がある。そこで、本実施例の電子密度推定装置は、３次元電子密度推定部３９により求められた任意の場所の電子密度を初期値としてステップＳ３からの処理を繰り返すことにより３次元電子密度推定部３９の近似精度を向上させることができる。
【００６５】
具体的には、３次元電子密度推定部３９は、（２２）式の値Ｅがしきい値ε以下になったか否かを判断する（ステップＳ６）。ここで、（２２）式で求めている差分は、暫定電子密度推定部３８により推定され、衛星信号が通過しているメッシュに対応した電子密度Ｘｋと、３次元電子密度推定部３９が（１７）式で求めた電子密度Ｖのうち、このメッシュに対応した電子密度Ｖとの差の２乗を全て加算したものである。
【数１３】
【００６６】
すなわち、３次元電子密度推定部３９は、推定した電子密度分布に基づく電子密度と暫定電子密度推定部３８により推定された電子密度との差分を領域毎に算出し、領域毎に算出した差分の２乗の合計値が所定のしきい値ε以下である場合に、推定した電子密度分布を真の値であると判断する。
【００６７】
（２２）式の値Ｅがしきい値εよりも大きい場合には、ステップＳ３に戻り、３次元電子密度推定空間設定部３６が３次元空間に３次元電子密度推定部３９により推定された値Ｖを設定する。したがって、暫定電子密度推定部３８は、３次元電子密度推定部３９により領域毎に算出された差分の２乗の合計値が所定のしきい値εよりも大きい場合に、３次元電子密度推定部３９により推定された電子密度分布に基づいて再度領域毎の電子密度を推定する。
【００６８】
以上のように、３次元電子密度推定空間設定部３６、暫定電子密度推定部３８、及び３次元電子密度推定部３９においてステップＳ３からステップＳ６までの処理が繰り返され、（２２）式の値Ｅがしきい値ε以下になった場合に、３次元電子密度推定部３９は、真値に近い電子密度を推定することができたとして、繰り返し処理を終了する。
【００６９】
上述したように、３次元電子密度推定部３９は、（１７）式を用いて同一高さに対する滑らかに近似された電子密度を求めることができる。次に、３次元電子密度推定部３９は、緯度・経度を指定した場合の高さ方向について、任意の高さの電子密度を求めるための処理を行う（ステップＳ７）。本実施例において、３次元電子密度推定部３９は、任意の高さの電子密度を求めるために、（１７）式により求められる関数値を使用して内挿して求める。すなわち、３次元電子密度推定部３９は、推定した高さ毎の電子密度分布を用いて隣接する領域間を内挿して任意の高さにおける電子密度分布を推定する。
【００７０】
なお、内挿の方法は、高さ方向についてのみの１次元的内挿あるいは求めたい場所を取り囲むメッシュを使用した３次元的内挿により求める。
【００７１】
まず、緯度・経度で指定した高さ方向のみの内挿について説明する。３次元電子密度推定部３９は、（１７）式に求めたい場所の緯度（θ）、経度（φ）を代入し、高さ方向のメッシュにおける電子密度（Ｖ（ｘ，φ，ｈ１）、Ｖ（ｘ，φ，ｈ２）、…、Ｖ（ｘ，φ，ｈｍａｘ））を求める。ここで、ｈ１はメッシュの最下層の高さを示し、ｈｍａｘはメッシュの最上層の高さを示す。このＶを使用したスプライン補間等、隣接した区間で１階微分及び２階微分が連続した方法で補間し、３次元電子密度推定部３９は、求めたい高さの電子密度を求める。
【００７２】
次に、求めたい場所を取り囲むメッシュによる内挿について説明する。３次元電子密度推定部３９は、全てのメッシュ点を（１７）式で求め（Ｖ（ｘ，φ，ｈ１）、Ｖ（ｘ，φ，ｈ２）、…、Ｖ（ｘ，φ，ｈｍａｘ））、３次元の内挿を行う。このＶを使用したスプライン補間等、隣接した区間で１階微分及び２階微分が連続した方法で補間し、３次元電子密度推定部３９は、求めたい高さの電子密度を求める。
【００７３】
上述のとおり、本発明の実施例１の形態に係る電子密度推定装置及び電子密度推定方法によれば、３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も行うことができる。
【００７４】
すなわち、本実施例の電子密度推定装置及び電子密度推定方法は、３次元電子密度推定部３９が球面調和関数を用いて高さ毎の電子密度分布を近似するので、水平方向に対して滑らかに連結された電子密度分布を得ることができ、ＧＰＳ等の衛星で観測できない空間を含む任意の緯度・経度における電子密度を推定することができる。
【００７５】
さらに、本実施例の電子密度推定装置及び電子密度推定方法は、３次元電子密度推定部３９において推定した電子密度分布に基づく電子密度と暫定電子密度推定部３８により推定された電子密度との差分を領域毎に算出し、領域毎に算出した差分の２乗の合計値が所定のしきい値ε以下になったか否かを判断基準とし、３次元電子密度推定部３９により求められた任意の場所の電子密度を初期値としてステップＳ３からＳ６までの処理を繰り返すことにより３次元電子密度推定部３９の近似精度を向上させることができる。
【００７６】
また、本実施例の電子密度推定装置及び電子密度推定方法は、３次元電子密度推定部３９で推定した高さ毎の電子密度分布を用いて隣接する領域間を内挿して任意の高さにおける電子密度分布を推定するので、水平方向のみならず高さ方向に対しても滑らかに連結された電子密度分布を得ることができ、従来手法において再現が困難だった電離層が有する特異な高さ方向の電子密度ピーク（Ｅ層やＦ層）を再現することができる。
【産業上の利用可能性】
【００７７】
本発明に係る電子密度推定装置及び電子密度推定方法は、衛星から送信される複数の異なる周波数信号による通過経路の電離層電子密度分布を推定する電子密度推定装置に利用可能である。
【符号の説明】
【００７８】
１０ａ、１０ｂ、１０ｃ 測位衛星
２０ 衛星信号処理系
２２ アンテナ
２４ 衛星信号受信機
２６ 衛星信号処理装置
３０ 電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系
３２ 電子密度推定装置
３３ 衛星ＴＥＣ算出部
３４ ＴＥＣバイアス補正部
３５ 周波数間バイアス推定部
３６ ３次元電子密度推定空間設定部
３７ 電子密度モデル算出部
３８ 暫定電子密度推定部
３９ ３次元電子密度推定部
４０ インターネットデータ処理系
４２ ルータ
４４ ＧＥＯＮＥＴ収集データ処理装置
４６ 外部インターネット網
５０ 通信回線
Ｓ１，Ｓ２，Ｓ３，Ｓ４ 測位衛星
Ｒｘ１，Ｒｘ２，Ｒｘ３ 衛星信号受信機
【技術分野】
【０００１】
本発明は、衛星から送信される複数の異なる周波数信号による通過経路の電離層電子密度分布を推定する電子密度推定装置及び電子密度推定方法に関する。
【背景技術】
【０００２】
衛星ナビゲーションシステム「Ｇａｌｉｌｅｏ」やＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）、ＱＺＳＳ（Ｑｕａｓｉ Ｚｅｎｉｔｈ Ｓａｔｅｌｌｉｔｅｓ Ｓｙｓｔｅｍ；準天頂衛星システム）等の衛星から送信される信号は、地球上の観測点の測位を行うのに利用される。これらの衛星により放送される複数周波数信号の擬似距離差は、電波の伝搬経路上の電子数に関係していることが知られている。
【０００３】
衛星から地上に向けて送信される信号は、電離層を通過する際に伝搬遅延量を生じる。この電離層伝搬遅延量は、伝搬する信号の周波数に依存している。そのため、衛星からの複数の周波数信号を受信することにより、信号の通過経路における総電子数（ＴＥＣ：Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ）を得ることができる。電波の伝搬経路上の総電子数を推定する場合には、複数周波数信号の送信・受信時刻による擬似距離（コード距離）差を使用する方法や複数周波数のキャリア位相擬似距離（フェーズ距離）の差を利用する方法が知られている。
【０００４】
非特許文献１や非特許文献２には、電離層における電子密度モデル関数であるＩＲＩ（Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ）モデルが記載されている。さらに、非特許文献４には、やはり電離層電子密度モデル関数の一つであるＧａｌｌａｇｈｅｒのモデルが記載されている。これらのモデル関数を用いて計算によりＴＥＣを求めることもできる。
【０００５】
非特許文献５や非特許文献６には、トモグラフィ手法による衛星ＴＥＣを用いた３次元電子密度推定が記載されており、ＭＡＲＴ（Ｍｕｌｔｉｐｌｉｃａｔｉｖｅ Ａｌｇｅｂｒａｉｃ Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ Ｔｅｃｈｎｉｑｕｅ）と呼ばれる代数的再構成法が適用されている。一方、非特許文献７には、ＧＳＶＤ（Ｇｅｎｅｒａｌｉｚｅｄ Ｓｉｎｇｕｌａｒ Ｖａｌｕｅ Ｄｅｃｏｍｐｏｓｉｔｉｏｎ）を用いたトモグラフィック再構成による電子密度推定が記載されている。
【０００６】
また、非特許文献８及び非特許文献９には、電離層を球核の薄い層とみなした電子密度計算について記載されている。さらに、非特許文献１０には、３次元方向の電子密度分布を高さ方向に計算した場合について記載されている。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ ２０００”，Ｒａｄｉｏ．Ｓｃｉｅｎｃｅ，Ｖｏｌ．３６，Ｎｕｍｂｅｒ２，ＰＰ２６１−２７５，Ｍａｒｃｈ／Ａｐｒｉｌ，２００１
【非特許文献２】Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ，ｅｔ，ａｌ．，：“Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ １９９０”，Ｎｏｖｅｍｂｅｒ，１９９０．
【非特許文献３】Ａ．Ｋｏｍｊａｔｈｙ：‘Ｇｌｏｂａｌ Ｉｏｎｏｓｐｈｅｒｉｃ Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ Ｍａｐｐｉｎｇ Ｕｓｉｎｇ ｔｈｅ Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ’ ，ＵＮＢ，Ｔｅｃｈｎｉｃａｌ Ｒｅｐｏｒｔ Ｎｏ．１８８，Ｓｅｐ．１９９７．
【非特許文献４】Ｇａｌｌａｇｈｅｒ，Ｄ．Ｌ．，Ｐ．Ｄ．Ｃｒａｖｅｎ，ａｎｄ Ｒ．Ｈ．Ｃｏｍｆｏｒｔ，Ｇｌｏｂａｌ ｃｏｒｅ ｐｌａｓｍａ ｍｏｄｅｌ，Ｊ．Ｇｅｏｐｈｙｓ．Ｒｅｓ．１０５，Ａ８，１８，８１９−１８，８３３，２０００．
【非特許文献５】Ｒａｙｍｕｎｄ，Ｔ．Ｄ．，Ａｕｓｔｉｎ，Ｊ．Ｒ．，Ｆｒａｎｋｅ，Ｓ．Ｊ．，Ｌｉｕ，Ｃ．Ｈ．，Ｋｌｏｂｕｃｈａｒ，Ｊ．Ａ．，ａｎｄ Ｓｔａｌｋｅｒ，Ｊ．，：“Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｃｏｍｕｔｅｒｉｚｅｄ Ｔｏｍｏｇｒａｐｈｙ ｔｏ ｔｈｅ Ｉｎｖｅｓｔｉｇａｔｉｏｎ ｏｆ Ｉｏｎｏｓｐｈｅｒｉｃ Ｓｔｒｕｃｔｕｒｅｓ”，Ｒａｄｉｏ Ｓｃｉ．，２５（５），７７１−７８９，１９９０．
【非特許文献６】Ｓｔｅｆａｎ Ｓｃｈｌｕｅｒ，ｅｔ ａｌ．：“Ｍｏｎｉｔｏｒｉｎｇ ｔｈｅ ３ Ｄｉｍｅｎｓｉｏｎａｌ Ｉｏｎｏｓｐｈｅｒｉｃ Ｅｌｅｃｔｒｏｎ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｂａｓｅｄ ｏｎ ＧＰＳ Ｍｅａｓｕｒｅｍｅｎｔｓ”
【非特許文献７】Ｋ．Ｂｈｕｙａｎ，ｅｔ ａｌ：“Ｔｏｍｏｇｒａｐｈｉｃｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｏｆ ｔｈｅ Ｉｏｎｏｓｐｈｅｒｅ ｕｓｉｎｇ ｇｅｎｅｒａｌｉｚｅｄ ｓｉｎｇｕｌａｒ ｖａｌｕｅ ｄｅｃｏｍｐｏｓｉｔｉｏｎ” ，Ｃｕｒｒｅｎｔ Ｓｃｉ．，Ｖｏｌ１．８３，Ｎｏ．９，１０ Ｎｏｖ．，２００２．
【非特許文献８】Ｏｔｕｋａ，Ｙ．，ｅｔ．ａｌ．：“Ａ ｎｅｗ ｔｅｃｈｎｉｑｕｅ ｆｏｒ ｍａｐｐｉｎｇ ｏｆ ｔｏｔａｌ ｅｌｅｃｔｒｏｎ ｃｏｎｔｅｎｔ ｕｓｉｎｇ ＧＰＳ ｎｅｔｗｏｒｋ ｉｎ Ｊａｐａｎ，”，Ｅａｒｔｈ Ｐｌａｎｅｔｓ Ｓｐａｃｅ，１１１−１２０，２００１．
【非特許文献９】Ｍａ，Ｇ．，ａｎｄ Ｔ．Ｍａｒｕｙａｍａ（２００３），Ｄｅｒｉｖａｔｉｏｎ ｏｆ ＴＥＣ ａｎｄ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｉｎｓｔｒｕｍｅｎｔａｌ ｂｉａｓｅｓ ｆｒｏｍ ＧＥＯＮＥＴ ｉｎ Ｊａｐａｎ，Ａｎｎ．Ｇｅｏｐｈｙｓ．，２１，２０８３−２０９３．
【非特許文献１０】五十嵐来良，斉藤昭則，大塚雄一：“標準電離層モデル（ＩＲＩ）で導出したＴＥＣとＧＰＳ受信機網で観測したＴＥＣの比較解析”，第１回電離圏の利用と影響に関するシンポジウム講和集，ｐｐ２４−１〜８，２００３．
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、非特許文献５や非特許文献６に記載された手法は、観測した衛星の数が限定されるため、限られた空間の領域にしか適用できないという短所がある。特に、非特許文献５に記載された手法は、空間平均を使用しており、電離層が有する特異な高さ方向の電子密度の谷（Ｅ層とＦ層の間にある密度が小さくなる部分）を消してしまうため、電離層の高さ方向変化を捉えることができない。また、非特許文献６には、空間平均についての記載がされていない。
【０００９】
非特許文献８及び非特許文献９に記載された手法は、上述したように電離層を球核の薄い層とみなして電子密度を計算しているため、３次元高さ方向の分布を知ることができない。また、非特許文献１０は、３次元方向の電子密度分布を高さ方向に計算しているが、水平方向の広がりについての言及が無く、緯度・経度方向の広がりがある空間の電子密度を求めていない。これらにより、従来の文献においては、ＧＰＳ等の衛星で観測できない空間部分の取り扱い方法が不十分である。
【００１０】
すなわち、従来の技術は、３次元の電離層電子密度分布を十分推定しきれていない。従来技術における手法は、空間平均を実施して電子密度を推定するため、本来なら通過する電波であっても反射したと判断してしまう場合がある。また、推定された空間とそれ以外の空間との差が大きい場合に、空間的な電子密度分布変化に不連続が生じてしまう。このように、従来手法は、電子密度変化に不連続が生じた結果、電子密度変化率を使用して電波伝搬経路を推定する際に大きな誤差を生じてしまうという問題がある。
【００１１】
本発明は上述した従来技術の問題点を解決するもので、３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も可能とする電子密度推定装置及び電子密度推定方法を提供することを課題とする。
【課題を解決するための手段】
【００１２】
本発明に係る電子密度推定装置は、上記課題を解決するために、複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定装置であって、前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定部と、前記総電子数算出部により算出された総電子数に基づいて前記空間設定部により分割された領域毎の電子密度を推定する第１電子密度推定部と、前記第１電子密度推定部により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定部とを備えることを特徴とする。
【００１３】
本発明に係る電子密度推定方法は、上記課題を解決するために、複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定方法であって、前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出ステップと、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定ステップと、前記総電子数算出ステップにより算出された総電子数に基づいて前記空間設定ステップにより分割された領域毎の電子密度を推定する第１電子密度推定ステップと、前記第１電子密度推定ステップにより推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定ステップとを備えることを特徴とする。
【発明の効果】
【００１４】
本発明によれば、３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も可能とする電子密度推定装置及び電子密度推定方法を提供することができる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施例１の形態の電子密度推定装置及び測位システムの構成を示すブロック図である。
【図２】本発明の実施例１の形態の電子密度推定装置の詳細なブロック図である。
【図３】本発明の実施例１の形態の電子密度推定装置の動作を示すフローチャート図である。
【図４】本発明の実施例１の形態の電子密度推定装置内の３次元電子密度推定空間設定部による３次元空間設定を説明する図である。
【図５】本発明の実施例１の形態の電子密度推定装置内の３次元電子密度推定空間設定部による補正係数設定を説明する図である。
【発明を実施するための形態】
【００１６】
以下、本発明の電子密度推定装置及び電子密度推定方法の実施の形態について図面を参照して詳細に説明する。
【実施例１】
【００１７】
以下、本発明の実施例について図面を参照しながら説明する。図１は、本発明の実施例１の電子密度推定装置３２及び測位システムの構成を示すブロック図である。この電子密度推定装置３２は、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される衛星信号を受信する１以上の衛星信号受信機２４を有し且つ衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する。
【００１８】
まず、本実施の形態の構成を説明する。本実施例の電子密度推定装置３２及び測位システムは、図１に示すように、測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系３０、及びインターネットデータ処理系４０で構成されている。測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃ、及び衛星信号処理系２０は、本発明の測位システムに対応する。なお、インターネットデータ処理系４０も含めて測位システムとすることも可能であるが、インターネットデータ処理系４０は、必ずしも必須のものではなく、付加的なものである。また、衛星信号処理系２０、電離層ＴＥＣ・バイアス推定処理系３０、及びインターネットデータ処理系４０は、お互いに通信回線５０で接続されている。
【００１９】
測位衛星１０ａ、測位衛星１０ｂ、測位衛星１０ｃは、ＧＰＳ、Ｇａｌｉｌｅｏ、準天頂衛星（ＱＺＳＳ）等の航法衛星であり、複数の周波数の衛星信号を送信する。
【００２０】
衛星信号処理系２０は、アンテナ２２、衛星信号受信機２４、及び衛星信号処理装置２６から構成される。衛星信号受信機２４は、本発明の受信機に対応し、複数の測位衛星（測位衛星１０ａ、１０ｂ、１０ｃ）から送信される複数の周波数の衛星信号をアンテナ２２を介して受信する。また、衛星信号処理装置２６は、衛星信号に含まれる測位情報を用いて位置情報を得る。
【００２１】
電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系３０は、電子密度推定装置３２を有し、受信した信号に基づき衛星位置、信号通過経路の電離層総電子数（ＴＥＣ）、周波数間バイアス等を算出する。受信したデータ・処理結果は、データサーバ系へＬＡＮ経由で伝送され保存される。
【００２２】
インターネットデータ処理系４０は、ルータ４２、ＧＥＯＮＥＴ収集データ処理装置４４、及び外部インターネット網４６で構成されている。ルータ４２は、スイッチングハブでもよい。インターネットデータ処理系４０は、公開されている電離層関連の情報や国土地理院が公開しているＧＰＳ観測データ（ＧＥＯＮＥＴデータ）や国際的に観測結果を公開しているＩＧＳ（Ｉｎｔｅｒｎａｔｉｏｎａｌ ＧＰＳ Ｓｅｒｖｉｃｅ ｆｏｒ Ｇｅｏｄｙｎａｍｉｃｓ）データ等をインターネット経由で収集する装置である。インターネットデータ処理系４０により収集された結果は、電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系３０で処理される。ルータ４２は、セキュリティを考慮して設けられ、ファイアウォールとする。
【００２３】
図２は、本発明の実施例１の電子密度推定装置３２の詳細なブロック図である。図２に示すように、電子密度推定装置３２は、衛星ＴＥＣ算出部３３、ＴＥＣバイアス補正部３４、周波数間バイアス推定部３５、３次元電子密度推定空間設定部３６、電子密度モデル算出部３７、暫定電子密度推定部３８、及び３次元電子密度推定部３９で構成されている。
【００２４】
衛星ＴＥＣ算出部３３、ＴＥＣバイアス補正部３４、及び周波数間バイアス推定部３５は、本発明の総電子数算出部に対応し、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。具体的には、衛星ＴＥＣ算出部３３は、衛星信号受信機２４から測位衛星１０ａ、１０ｂ、１０ｃの各々に対する擬似距離に基づき衛星信号の通過経路の総電子数を算出する。
【００２５】
周波数間バイアス推定部３５は、周波数間バイアスを算出し、ＴＥＣバイアス補正部３４に出力する。周波数間バイアス推定部３５による周波数間バイアスの推定方法は、どのようなものでもよく、従来手法を用いてよい。
【００２６】
ＴＥＣバイアス補正部３４は、周波数間バイアス推定部３５により推定された周波数間バイアスを用いて補正を行い、最終的な総電子数を算出する。
【００２７】
３次元電子密度推定空間設定部３６は、本発明の空間設定部に対応し、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する。
【００２８】
電子密度モデル算出部３７は、上述したようなＩＲＩモデルやＧａｌｌａｇｈｅｒモデル等の電離層電子密度モデル関数に基づき電子密度モデルを算出する。
【００２９】
暫定電子密度推定部３８は、本発明の第１電子密度推定部に対応し、総電子数算出部（ＴＥＣバイアス補正部３４）により算出された総電子数に基づいて３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する。本実施例において、暫定電子密度推定部３８は、非特許文献５や非特許文献に記載されているようなＭＡＲＴの手法を用いて電子密度を推定するものとする。
【００３０】
３次元電子密度推定部３９は、本発明の第２電子密度推定部に対応し、暫定電子密度推定部３８により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する。
【００３１】
次に、上述のように構成された本実施の形態の作用を説明する。図３は、本実施例の電子密度推定装置３２の動作を示すフローチャート図である。まず、衛星ＴＥＣ算出部３３は、１以上の受信機の各々から複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき衛星信号の通過経路の総電子数を算出する（ステップＳ１）。具体的には、衛星ＴＥＣ算出部３３は、広範囲に設置した衛星受信機収集データからＧＰＳ等の衛星のＴＥＣを計算する。ここで、本実施例の衛星ＴＥＣ算出部３３は、衛星から放送される２周波の信号（ＧＰＳ衛星のＬ１周波数（１５７５．４２ＭＨｚ）とＬ２周波数（１２２７．６０ＭＨｚ）の信号）に基づいて信号伝搬経路のＴＥＣを算出するものとする。ＧＰＳ以外でも異なる２周波を放送している場合には同様な手法が使用できる。
【００３２】
ここでは、測位衛星１０ａと衛星信号受信機２４との間の擬似距離について考える。衛星信号通過経路のＴＥＣは、測位衛星１０ａから送信される衛星信号（ここではＬ１とＬ２で示す）に基づき求められる。まず、擬似距離（コード距離、シュードレンジ）とキャリア位相擬似距離（フェーズ距離）は、以下のように表すことができる。
【数１】
【００３３】
ここで、ρは、遅延時間による距離（擬似距離）を示す。また、Φは、キャリア位相距離を示す。ｒは、真の距離を示す。ｃは光速である。また、δｔｕは、受信機時刻誤差を示し、δｔｓは、衛星時刻誤差を示す。なお、本実施例において、右下の添え字は基本的に地上（受信機等）に関連する項を示し、右上の添え字は基本的に上空（衛星等）に関連する項を示す。δｔｕ、Ｌ１ｏｒＬ２ ｂｉａｓは、受信機周波数依存ハードウェア依存バイアスを示し、δｔｓＬ１ｏｒＬ２ ｂｉａｓは、衛星周波数依存ハードウェア依存バイアスを示す。さらに、Ｉは、電離層伝搬遅延量を示し、Ｔは、対流圏伝搬遅延量を示す。ＮＬ１ｏｒＬ２，ａｍｂは、整数不確定値を示し、εは、観測誤差を示す。
【００３４】
衛星ＴＥＣ算出部３３は、２つの周波数の観測値の差（遅延時間の差）をとることにより、ＴＥＣを求めることができる。
【数２】
【００３５】
ここで、ｆは周波数である。また、λは波長を示す。Δｂｉａｓは、周波数間バイアスを示す。また、ＴＥＣｔｒｕｅは、真の総電子数を示す。
【００３６】
次に、ＴＥＣバイアス補正部３４は、周波数間バイアス推定部３５により事前に推定された周波数間バイアスΔｂｉａｓによる補正をＴＥＣに実施する（ステップＳ２）。まず、バイアス補正後の２つの周波数の遅延時間差は、次式により表される。
【数３】
【００３７】
ここで、（６）式には、観測値から分からない不確定値（λＬ１ΔＮＬ１，ａｍｂ−λＬ２ΔＮＬ２，ａｍｂ）が含まれている。そこで、ＴＥＣバイアス補正部３４は、（５）式と（６）式とを組み合わせて、以下に示す式により整数不確定値を削除する。この手法は、非特許文献３にも記載されている。
【数４】
【００３８】
ここで、添え字ｋは、データの番号を示す。また、Ｍは、連続的に収集できたサンプル数の合計を示す。さらに、ＴＥＣｋ，ｕは、衛星ｕの伝搬経路上の総電子数（個／ｍ２）を示す。
【００３９】
電子密度モデル算出部３７は、電離層電子密度モデル関数に基づき電子密度モデルを算出し、算出した電子密度モデルを３次元電子密度推定空間設定部３６に出力する。なお、電子密度モデル算出部３７による電離層モデルは、例えば非特許文献１，２，４等に記載されている。また、ステップＳ１，Ｓ２（特にＳ１）は、本発明の総電子数算出ステップに対応する。
【００４０】
次に、３次元電子密度推定空間設定部３６は、電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する。この３次元電子密度推定空間設定部３６の動作は、本発明の空間設定ステップに対応する。また、３次元電子密度推定空間設定部３６は、推定する３次元空間の電離層電子密度モデル値あるいは推定値を設定する（ステップＳ３）。
【００４１】
具体的には、３次元電子密度推定空間設定部３６は、衛星から衛星受信機が配置されている地点まで３次元空間における電離層を、緯度・経度・高さ方向に分割し、電離層モデル値を使用して分割した領域内あるいは格子点の位置の電子密度を設定する。非特許文献１，２等に記載されたＩＲＩモデルは、高さ１０００ｋｍまでが有効範囲であるため、衛星が存在する空間までは設定できない。そこで、３次元電子密度推定空間設定部３６は、例えば非特許文献４に記載されたＧＣＰＭ（Ｇｌｏｂａｌ Ｃｏｒｅ Ｐｌａｓｍａ Ｍｏｄｅｌ）を用いる手法や、ＩＲＩモデルの約４００ｋｍ〜１０００ｋｍを外挿する手法を採用することができる。
【００４２】
図４は、本実施例の電子密度推定装置３２内の３次元電子密度推定空間設定部３６による３次元空間設定を説明する図である。本実施例の３次元電子密度推定空間設定部３６は、図４に示すように、３次元空間をＸ１〜Ｘ１２の１２個の領域（ブロック）に分割しているが、実際には膨大な数のブロック数となる。図４において、測位衛星Ｓ１〜Ｓ４と衛星信号受信機Ｒｘ１〜Ｒｘ３との間における信号通過経路は、矢印により表されている。図４の衛星・ＴＥＣ・電子密度の関係は、次式により表される。
【数５】
【００４３】
（９）式において、Ｘ１〜Ｘ１２は、空間を分割した領域（メッシュ）の電子密度を示し、現時点では不明である。また、ＴＥＣｉｊは、ｉ衛星・ｊ受信機間で決まる総電子数を示す。例えば、図４中の太い矢印ＴＥＣ２１は、測位衛星Ｓ２と衛星信号受信機Ｒｘ１との間の信号通過経路における総電子数を示す。ＴＥＣバイアス補正部３４により出力されたＴＥＣの値は、（９）式の右辺に代入される。
【００４４】
ここで、（９）式の行列式をＡＸ＝Ｅと表すことにする。行列Ａの中のゼロでない要素（１で表されている要素）は、ｉ衛星・ｊ受信機間の信号が通過したメッシュを示す。
【００４５】
実際にモデルにより設定した電子密度の値は、メッシュ内の１点の値である。また、電波信号が各メッシュを通過する際に、電波はメッシュを縦に通過する場合のみならず斜めに通過する場合もある。この効果を考慮し、３次元電子密度推定空間設定部３６は、行列Ａの１の代わりに次式で示される値を代入する。
【数６】
【００４６】
ここで、ｄＨｉ，ｊは、ｉ，ｊメッシュの分割した高さ方向の大きさを示す。また、変換係数ＳＦは、斜めから入射した信号を高さ方向に入射した信号に補正する係数である。θは、衛星受信機が受信した衛星電波の入射仰角を示す。図５は、本実施例の電子密度推定装置３２内の３次元電子密度推定空間設定部３６による補正係数ＳＦ設定を説明する図である。
【数７】
【００４７】
（１２）式の左辺は、ｉ，ｊメッシュ内の垂直方向のＴＥＣを示す。また、右辺のＴＥＣｓは、ｉ，ｊメッシュ内を通過する電波方向のＴＥＣを示す。（１３）式内のＲは、図５に示すように、地球半径を示す。また、（１３）式内のｈは、地表面からｉ，ｊメッシュまでの高さを示す。
【００４８】
次に、暫定電子密度推定部３８は、ＭＡＲＴあるいはＧＳＶＤ手法により、衛星信号通過経路の電子密度を推定する（ステップＳ４）。具体的には、暫定電子密度推定部３８は、総電子数算出部（ＴＥＣバイアス補正部３４）により算出された総電子数に基づいて３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する。この暫定電子密度推定部３８の動作は、本発明の第１電子密度推定ステップに対応する。
【００４９】
なお、本実施例において、暫定電子密度推定部３８は、ＭＡＲＴあるいはＧＳＶＤ手法により電子密度を推定するが、本発明を実現するためには、必ずしもいずれかの手法を使う必要はなく、衛星信号受信機が得た観測値に基づいて暫定的に各領域の電子密度を推定できればよい。
【００５０】
ＭＡＲＴの手法は、上述したように非特許文献５，６等に記載されているが、ここで簡潔に説明する。
【数８】
【００５１】
（１４）式は、ＭＡＲＴのアルゴリズムを示す式である。ａｉは、（１０）式のＡ行列ｉ行ベクトルを表す。また、ａｉｊは、Ａ行列のｉｊ成分を表す。Ｘｊは、Ｘ列行列のｊ成分を表す。さらに、Ｅｉは、（９）式のＥ行列のｉ成分を表している。また、λは、収束の度合いを決めるパラメータである。
【００５２】
暫定電子密度推定部３８は、（１４）式に示すＭＡＲＴのアルゴリズムを使用して、（１４）式の値が収束するまで繰り返し処理を実施し、３次元電子密度推定空間設定部３６により分割された領域毎の電子密度を推定する。その際に、暫定電子密度推定部３８は、ＴＥＣバイアス補正部３４により出力されたＴＥＣの値をＥ行列の値に使用するとともに、電子密度モデル算出部３７により算出されたモデル値をＸの初期値として使用する。
【００５３】
なお、暫定電子密度推定部３８は、信号が全く通過していない空間の電子密度を推定することはできない。
【００５４】
ＧＳＶＤ手法を採用する場合には、暫定電子密度推定部３８は、以下に示す方式で電子密度を求める。なお、ＧＳＶＤ手法については、非特許文献７に記載がある。
【数９】
【００５５】
（１５）式のＬは、Ｔｉｋｈｏｎｏｖ行列と呼ばれる。０次のＴｉｋｈｏｎｏｖ正規化では、Ｌを単位行列とする。αは、正規化パラメータで、α＞０とする。また、Ｌは、正値２次形式行列（固有値が全て正）である。
【００５６】
（１５）式を満足する電子密度Ｘは、次式で表すことができる。
【数１０】
【００５７】
次に、３次元電子密度推定部３９は、高さ毎に球面調和関数により電子密度を近似する（ステップＳ５）。すなわち、３次元電子密度推定部３９は、暫定電子密度推定部３８により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する。この３次元電子密度推定部３９の動作は、本発明の第２電子密度推定ステップに対応する。
【００５８】
なお、３次元電子密度推定部３９は、その作用効果において本発明の特徴を示すものであり、従来手法と異なる点である。
【００５９】
（９）式のＸは、空間のメッシュを一次元で表現したものとなっている。３次元電子密度推定部３９は、暫定電子密度推定部３８により推定された電子密度の中から、同じ高さの電子密度だけを抽出する。
【数１１】
【００６０】
（１７）式において、ｘ＝ｓｉｎ（θ）である。Ｖは、ある高さにおける電子密度を示す。ただし、θはメッシュの緯度を示し、φはメッシュの経度を示す。また、Ｓｎ，Ｓｎｍは、Ｓｃｈｕｍｉｄｔ Ｓｅｍｉｎｏｒｍａｌａｉｚｅｄルジャンドル陪関数である。（２０）式にルジャンドル関数を示す。
【数１２】
【００６１】
Ｎｐは、ルジャンドル関数の次数を示す。３次元電子密度推定部３９は、同じ高さの電子密度を（１７）式を使用して近似する。これにより、ＭＡＲＴ手法で求められなかったメッシュの電子密度の値を求めることができる。またＧＳＶＤで求めた値を修正することができる。
【００６２】
すなわち、３次元電子密度推定部３９は、ある緯度・経度（θ，φ）について、高さ方向のメッシュ分（１７）式で電子密度を近似し、高さＨに対する係数（α（Ｈ）、β（Ｈ））を求める。この係数（α（Ｈ）、β（Ｈ））が求まることにより、３次元電子密度推定部３９は、緯度・経度・高さを指定し、（１７）式を用いて、その高さＨにおける電子密度Ｖ（ｘ，φ）を算出することができる。
【００６３】
なお、非特許文献８や非特許文献９に記載された従来手法は、ある高さのみに着目しており、高さ毎の電子密度分布を求めていない。また、非特許文献１０に記載された手法は、高さ方向に対する電子密度分布を求めているが、緯度・経度方向の広がりがある空間の電子密度を求めていない。
【００６４】
電子密度モデル算出部３７による電子密度モデル値と暫定電子密度推定部３８による電子密度の値とが大きくかけ離れている場合に、３次元電子密度推定部３９による近似は、適切でない場合がある。そこで、本実施例の電子密度推定装置は、３次元電子密度推定部３９により求められた任意の場所の電子密度を初期値としてステップＳ３からの処理を繰り返すことにより３次元電子密度推定部３９の近似精度を向上させることができる。
【００６５】
具体的には、３次元電子密度推定部３９は、（２２）式の値Ｅがしきい値ε以下になったか否かを判断する（ステップＳ６）。ここで、（２２）式で求めている差分は、暫定電子密度推定部３８により推定され、衛星信号が通過しているメッシュに対応した電子密度Ｘｋと、３次元電子密度推定部３９が（１７）式で求めた電子密度Ｖのうち、このメッシュに対応した電子密度Ｖとの差の２乗を全て加算したものである。
【数１３】
【００６６】
すなわち、３次元電子密度推定部３９は、推定した電子密度分布に基づく電子密度と暫定電子密度推定部３８により推定された電子密度との差分を領域毎に算出し、領域毎に算出した差分の２乗の合計値が所定のしきい値ε以下である場合に、推定した電子密度分布を真の値であると判断する。
【００６７】
（２２）式の値Ｅがしきい値εよりも大きい場合には、ステップＳ３に戻り、３次元電子密度推定空間設定部３６が３次元空間に３次元電子密度推定部３９により推定された値Ｖを設定する。したがって、暫定電子密度推定部３８は、３次元電子密度推定部３９により領域毎に算出された差分の２乗の合計値が所定のしきい値εよりも大きい場合に、３次元電子密度推定部３９により推定された電子密度分布に基づいて再度領域毎の電子密度を推定する。
【００６８】
以上のように、３次元電子密度推定空間設定部３６、暫定電子密度推定部３８、及び３次元電子密度推定部３９においてステップＳ３からステップＳ６までの処理が繰り返され、（２２）式の値Ｅがしきい値ε以下になった場合に、３次元電子密度推定部３９は、真値に近い電子密度を推定することができたとして、繰り返し処理を終了する。
【００６９】
上述したように、３次元電子密度推定部３９は、（１７）式を用いて同一高さに対する滑らかに近似された電子密度を求めることができる。次に、３次元電子密度推定部３９は、緯度・経度を指定した場合の高さ方向について、任意の高さの電子密度を求めるための処理を行う（ステップＳ７）。本実施例において、３次元電子密度推定部３９は、任意の高さの電子密度を求めるために、（１７）式により求められる関数値を使用して内挿して求める。すなわち、３次元電子密度推定部３９は、推定した高さ毎の電子密度分布を用いて隣接する領域間を内挿して任意の高さにおける電子密度分布を推定する。
【００７０】
なお、内挿の方法は、高さ方向についてのみの１次元的内挿あるいは求めたい場所を取り囲むメッシュを使用した３次元的内挿により求める。
【００７１】
まず、緯度・経度で指定した高さ方向のみの内挿について説明する。３次元電子密度推定部３９は、（１７）式に求めたい場所の緯度（θ）、経度（φ）を代入し、高さ方向のメッシュにおける電子密度（Ｖ（ｘ，φ，ｈ１）、Ｖ（ｘ，φ，ｈ２）、…、Ｖ（ｘ，φ，ｈｍａｘ））を求める。ここで、ｈ１はメッシュの最下層の高さを示し、ｈｍａｘはメッシュの最上層の高さを示す。このＶを使用したスプライン補間等、隣接した区間で１階微分及び２階微分が連続した方法で補間し、３次元電子密度推定部３９は、求めたい高さの電子密度を求める。
【００７２】
次に、求めたい場所を取り囲むメッシュによる内挿について説明する。３次元電子密度推定部３９は、全てのメッシュ点を（１７）式で求め（Ｖ（ｘ，φ，ｈ１）、Ｖ（ｘ，φ，ｈ２）、…、Ｖ（ｘ，φ，ｈｍａｘ））、３次元の内挿を行う。このＶを使用したスプライン補間等、隣接した区間で１階微分及び２階微分が連続した方法で補間し、３次元電子密度推定部３９は、求めたい高さの電子密度を求める。
【００７３】
上述のとおり、本発明の実施例１の形態に係る電子密度推定装置及び電子密度推定方法によれば、３次元の高さ方向及び水平方向の電離層電子密度分布をより正確に推定し、衛星で観測できない空間部分の推定も行うことができる。
【００７４】
すなわち、本実施例の電子密度推定装置及び電子密度推定方法は、３次元電子密度推定部３９が球面調和関数を用いて高さ毎の電子密度分布を近似するので、水平方向に対して滑らかに連結された電子密度分布を得ることができ、ＧＰＳ等の衛星で観測できない空間を含む任意の緯度・経度における電子密度を推定することができる。
【００７５】
さらに、本実施例の電子密度推定装置及び電子密度推定方法は、３次元電子密度推定部３９において推定した電子密度分布に基づく電子密度と暫定電子密度推定部３８により推定された電子密度との差分を領域毎に算出し、領域毎に算出した差分の２乗の合計値が所定のしきい値ε以下になったか否かを判断基準とし、３次元電子密度推定部３９により求められた任意の場所の電子密度を初期値としてステップＳ３からＳ６までの処理を繰り返すことにより３次元電子密度推定部３９の近似精度を向上させることができる。
【００７６】
また、本実施例の電子密度推定装置及び電子密度推定方法は、３次元電子密度推定部３９で推定した高さ毎の電子密度分布を用いて隣接する領域間を内挿して任意の高さにおける電子密度分布を推定するので、水平方向のみならず高さ方向に対しても滑らかに連結された電子密度分布を得ることができ、従来手法において再現が困難だった電離層が有する特異な高さ方向の電子密度ピーク（Ｅ層やＦ層）を再現することができる。
【産業上の利用可能性】
【００７７】
本発明に係る電子密度推定装置及び電子密度推定方法は、衛星から送信される複数の異なる周波数信号による通過経路の電離層電子密度分布を推定する電子密度推定装置に利用可能である。
【符号の説明】
【００７８】
１０ａ、１０ｂ、１０ｃ 測位衛星
２０ 衛星信号処理系
２２ アンテナ
２４ 衛星信号受信機
２６ 衛星信号処理装置
３０ 電離層ＴＥＣ・バイアス推定・３次元電子密度推定処理系
３２ 電子密度推定装置
３３ 衛星ＴＥＣ算出部
３４ ＴＥＣバイアス補正部
３５ 周波数間バイアス推定部
３６ ３次元電子密度推定空間設定部
３７ 電子密度モデル算出部
３８ 暫定電子密度推定部
３９ ３次元電子密度推定部
４０ インターネットデータ処理系
４２ ルータ
４４ ＧＥＯＮＥＴ収集データ処理装置
４６ 外部インターネット網
５０ 通信回線
Ｓ１，Ｓ２，Ｓ３，Ｓ４ 測位衛星
Ｒｘ１，Ｒｘ２，Ｒｘ３ 衛星信号受信機
【特許請求の範囲】
【請求項１】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定装置であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、
電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定部と、
前記総電子数算出部により算出された総電子数に基づいて前記空間設定部により分割された領域毎の電子密度を推定する第１電子密度推定部と、
前記第１電子密度推定部により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定部と、
を備えることを特徴とする電子密度推定装置。
【請求項２】
前記第２電子密度推定部は、推定した電子密度分布に基づく電子密度と前記第１電子密度推定部により推定された電子密度との差分を前記領域毎に算出し、領域毎に算出した差分の２乗の合計値が所定のしきい値以下である場合に推定した電子密度分布を真の値であると判断し、
前記第１電子密度推定部は、前記第２電子密度推定部により領域毎に算出された差分の２乗の合計値が所定のしきい値よりも大きい場合に、前記第２電子密度推定部により推定された電子密度分布に基づいて再度前記領域毎の電子密度を推定することを特徴とする請求項１記載の電子密度推定装置。
【請求項３】
前記第２電子密度推定部は、推定した高さ毎の電子密度分布を用いて隣接する領域間を内挿して任意の高さにおける電子密度分布を推定することを特徴とする請求項１又は請求項２記載の電子密度推定装置。
【請求項４】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定方法であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出ステップと、
電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定ステップと、
前記総電子数算出ステップにより算出された総電子数に基づいて前記空間設定ステップにより分割された領域毎の電子密度を推定する第１電子密度推定ステップと、
前記第１電子密度推定ステップにより推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定ステップと、
を備えることを特徴とする電子密度推定方法。
【請求項１】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定装置であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出部と、
電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定部と、
前記総電子数算出部により算出された総電子数に基づいて前記空間設定部により分割された領域毎の電子密度を推定する第１電子密度推定部と、
前記第１電子密度推定部により推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定部と、
を備えることを特徴とする電子密度推定装置。
【請求項２】
前記第２電子密度推定部は、推定した電子密度分布に基づく電子密度と前記第１電子密度推定部により推定された電子密度との差分を前記領域毎に算出し、領域毎に算出した差分の２乗の合計値が所定のしきい値以下である場合に推定した電子密度分布を真の値であると判断し、
前記第１電子密度推定部は、前記第２電子密度推定部により領域毎に算出された差分の２乗の合計値が所定のしきい値よりも大きい場合に、前記第２電子密度推定部により推定された電子密度分布に基づいて再度前記領域毎の電子密度を推定することを特徴とする請求項１記載の電子密度推定装置。
【請求項３】
前記第２電子密度推定部は、推定した高さ毎の電子密度分布を用いて隣接する領域間を内挿して任意の高さにおける電子密度分布を推定することを特徴とする請求項１又は請求項２記載の電子密度推定装置。
【請求項４】
複数の測位衛星から送信される衛星信号を受信する１以上の受信機を有し、前記衛星信号に含まれる測位情報を用いて位置情報を得る測位システムを利用する電子密度推定方法であって、
前記１以上の受信機の各々から前記複数の測位衛星の各々に対する擬似距離およびキャリア位相擬似距離に基づき前記衛星信号の通過経路の総電子数を算出する総電子数算出ステップと、
電子密度を推定する３次元空間を設定し、当該３次元空間を複数の領域に分割する空間設定ステップと、
前記総電子数算出ステップにより算出された総電子数に基づいて前記空間設定ステップにより分割された領域毎の電子密度を推定する第１電子密度推定ステップと、
前記第１電子密度推定ステップにより推定された電子密度に基づいて、同じ高さの電子密度を抽出し、球面調和関数を用いて高さ毎の電子密度分布を近似して推定する第２電子密度推定ステップと、
を備えることを特徴とする電子密度推定方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図２】
【図３】
【図４】
【図５】
【公開番号】特開２０１１−１３７６９８（Ｐ２０１１−１３７６９８Ａ）
【公開日】平成２３年７月１４日（２０１１．７．１４）
【国際特許分類】
【出願番号】特願２００９−２９７３０９（Ｐ２００９−２９７３０９）
【出願日】平成２１年１２月２８日（２００９．１２．２８）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【Ｆターム（参考）】
【公開日】平成２３年７月１４日（２０１１．７．１４）
【国際特許分類】
【出願日】平成２１年１２月２８日（２００９．１２．２８）
【出願人】（０００００３０７８）株式会社東芝 (54,554)
【Ｆターム（参考）】
[ Back to top ]