JP4861131B2 - 電離層電子密度算出装置 - Google Patents

Description

本発明は、電離層中の任意の位置における電子密度を算出する電離層電子密度算出装置に関する。

電波による通信を行う場合、通信の効率を良好に保つために、通信波の伝搬経路をできるだけ正確に決める必要がある。ところで、電離層の影響を受けるＨＦ帯などの周波数の信号により通信を行う場合、通信波が電離層を通過する際に、空間に分布している電子により通信波の伝搬経路が屈折する。このため、伝搬経路の算出に際しては、電離層の影響を考慮しなければならない。

電離層を通過する通信波の伝搬経路の屈折の仕方は、電離層の電子密度の分布により決まるため、伝搬経路の算出の過程において、電離層中の所望の位置における電子密度を算出する必要がある。電離層中の任意の位置における電子密度は、既存の電離層電子密度分布モデルを用いて求めることができる。この電離層電子密度分布モデルに、電子密度を算出したい日時、位置などのデータを入力することにより、その位置の電子密度を算出することができる。

電離層電子密度分布モデルとしては、非特許文献１に記載されているＩＲＩ（Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ）モデルやＢｅｎｔモデルなどがある。
Ｄｉｅｔｅｒ Ｂｉｌｉｔｚａ，ｅｔ．ａｌ．，"Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｒｅｆｅｒｅｎｃｅ Ｉｏｎｏｓｐｈｅｒｅ １９９０"，Ｎｏｖｅｍｂｅｒ，１９９０

しかし、上記ＩＲＩモデルなどは、月平均レベルの電離層電子密度分布モデルであるため、ある時間における電離層の状態を正確に示すことはできない。このため、ＩＲＩモデルなどでは、現実の電離層の状態に即した電子密度を算出できないことがあった。

本発明は上記に鑑みてなされたもので、現実の電離層の状態に即して、電離層中の所望の位置における電子密度を算出することができる電離層電子密度算出装置を提供することを目的とする。

上記目的を達成するため、本発明の電離層電子密度算出装置は、送信される位置および送信周波数が既知であり、電離層を通過して到来する電波を受信するアレイアンテナと、このアレイアンテナからの出力信号から前記電波の到来方位角、到来仰角を第１の到来方位角、第１の到来仰角として算出する第１の方向算出手段と、電離層の電子密度分布を示す電離層電子密度分布モデルに入力する複数のパラメータをその初期値から変化させるための変化量を乱数により生成する変化量生成手段と、この変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記アレイアンテナで受信した前記電波の伝搬経路を算出し、この伝搬経路に基づいて前記電波の到来方位角、到来仰角を第２の到来方位角、第２の到来仰角として算出する第２の方向算出手段と、航法衛星から送信される複数周波数の衛星信号を受信する衛星信号受信手段と、この衛星信号受信手段で受信した前記複数周波数の衛星信号の電離層での周波数による伝搬遅延量の差を用いて、前記衛星信号の通過経路における総電子数を第１の総電子数として算出する第１の総電子数算出手段と、前記変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記衛星信号の通過経路における総電子数を第２の総電子数として算出する第２の総電子数算出手段と、前記第１の到来方位角と前記第２の到来方位角との差である第１の誤差を算出し、前記第１の到来仰角と前記第２の到来仰角との差である第２の誤差を算出し、前記第１の総電子数と前記第２の総電子数との差である第３の誤差を算出するとともに、前記第１〜第３の誤差を合わせた値の大きさの指標である誤差指標を算出する誤差指標算出手段と、この誤差指標算出手段で算出した前記誤差指標が所定の閾値以下であるか否かを判断し、前記誤差指標が前記所定の閾値以下のとき、当該誤差指標の算出に用いた前記第２の到来方位角、前記第２の到来仰角、および前記第２の総電子数の算出に用いた変化量の分だけ所望の位置における初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記所望の位置の電子密度を算出し、前記誤差指標が前記所定の閾値より大きいときは、前記変化量生成手段に新たな変化量を生成させる電子密度算出手段とを備えることを特徴とする。

また、本発明の電離層電子密度算出装置は、送信される位置および送信周波数が既知であり、電離層を通過して到来する電波を受信するアレイアンテナと、このアレイアンテナからの出力信号から前記電波の到来方位角、到来仰角を第１の到来方位角、第１の到来仰角として算出する第１の方向算出手段と、電離層の電子密度分布を示す電離層電子密度分布モデルに入力する複数のパラメータをその初期値から変化させるための変化量を、前記複数のパラメータの数より多い所定の個数だけ乱数により生成する第１の変化量生成手段と、この第１の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記アレイアンテナで受信した前記電波の伝搬経路を算出し、この伝搬経路に基づいて前記電波の到来方位角、到来仰角を第２の到来方位角、第２の到来仰角として算出する処理を、前記所定の個数の変化量のそれぞれについて行う第２の方向算出手段と、航法衛星から送信される複数周波数の衛星信号を受信する衛星信号受信手段と、この衛星信号受信手段で受信した前記複数周波数の衛星信号の電離層での周波数による伝搬遅延量の差を用いて、前記衛星信号の通過経路における総電子数を第１の総電子数として算出する第１の総電子数算出手段と、前記第１の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記衛星信号の通過経路における総電子数を第２の総電子数として算出する処理を、前記所定の個数の変化量のそれぞれについて行う第２の総電子数算出手段と、前記第１の到来方位角と前記第２の到来方位角との差である第１の誤差を算出し、前記第１の到来仰角と前記第２の到来仰角との差である第２の誤差を算出し、前記第１の総電子数と前記第２の総電子数との差である第３の誤差を算出するとともに、前記第１〜第３の誤差を合わせた値の大きさの指標である誤差指標を算出する処理を、前記所定の個数の変化量のそれぞれについて行う第１の誤差指標算出手段と、この第１の誤差指標算出手段において前記所定の個数の変化量のそれぞれに対応して算出された第１の誤差指標に基づいて、前記複数のパラメータを前記変化量に応じて変化させたときに第１の誤差指標がどのように変化するかを決める変数を、最小二乗法を用いて算出する変数算出手段と、この変数算出手段で算出した前記変数を用いて、前記各パラメータをその初期値から変化させるための変化量を乱数により生成する第２の変化量生成手段と、この第２の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記アレイアンテナで受信した前記電波の伝搬経路を算出し、この伝搬経路に基づいて前記電波の到来方位角、到来仰角を第３の到来方位角、第３の到来仰角として算出する第３の方向算出手段と、前記第２の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記衛星信号の通過経路における総電子数を第３の総電子数として算出する第３の総電子数算出手段と、前記第１の到来方位角と前記第３の到来方位角との差である第４の誤差を算出し、前記第１の到来仰角と前記第３の到来仰角との差である第５の誤差を算出し、前記第１の総電子数と前記第３の総電子数との差である第６の誤差を算出するとともに、前記第４〜第６の誤差を合わせた値の大きさの指標である第２の誤差指標を算出する第２の誤差指標算出手段と、この第２の誤差指標算出手段で算出した前記第２の誤差指標が所定の閾値以下であるか否かを判断し、前記第２の誤差指標が前記所定の閾値以下のとき、当該第２の誤差指標の算出に用いた前記第３の到来方位角、前記第３の到来仰角、および前記第３の総電子数の算出に用いた変化量の分だけ所望の位置における初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記所望の位置の電子密度を算出し、前記第２の誤差指標が前記所定の閾値より大きいときは、第２の変化量生成手段に新たな変化量を生成させる電子密度算出手段とを備えることを特徴とする。

本発明の電離層電子密度算出装置によれば、現実の電離層の状態に即して、電離層中の所望の位置における電子密度を算出することができる。

以下、本発明を実施するための最良の形態について、図面を参照して説明する。

図１は本発明の実施の形態に係る電離層電子密度算出装置の構成を示すブロック図である。図１に示すように本実施の形態に係る電離層電子密度算出装置は、衛星観測部１と、電離層観測部２と、インターネットデータ処理部３と、信号処理部４と、データサーバ部５とを備える。

衛星観測部１は、衛星信号受信アンテナ１１と、衛星信号受信部１２と、衛星信号処理部１３とを備える。衛星信号受信部１２は、衛星信号受信アンテナ１１を介して受信したＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）衛星、Ｇａｌｉｌｅｏ衛星、Ｇｌｏｎａｓｓ衛星、準天頂衛星などの航法衛星から送信される複数周波数の衛星信号を衛星信号処理部１３へ供給する。

衛星信号処理部１３は、衛星信号受信部１２からの複数周波数の衛星信号から、航法衛星の位置を算出するとともに、各周波数の伝搬遅延量の違いを用いて、衛星信号の通過経路に存在する総電子数（ＴＥＣ：Ｔｏｔａｌ Ｅｌｅｃｔｒｏｎ Ｃｏｎｔｅｎｔ）を算出し、得られた算出結果をネットワーク７を介して信号処理部４およびデータサーバ部５へ送信する。また、衛星信号処理部１３は、衛星信号から得られる各種データをネットワーク７を介してデータサーバ部５へ送信する。

なお、電離層を通過する際に生じる周波数による伝搬遅延量の違いにより衛星信号の通過経路の総電子数（ＴＥＣ）を算出する手法については、例えば、 “ＧＬＯＢＡＬ ＰＯＳＩＴＩＯＮＩＧ ＳＹＳＴＥＭ Ｓｉｇｎａｌｓ，Ｍｅａｓｕｒｅｍｅｎｔｓ，ａｎｄ Ｐｅｒｆｏｒｍａｎｃｅ．”（Ｐｒａｔａｐ Ｍｉｓｒａ，ＰｅｒＥｎｇｅ，Ｇａｎｇａ−Ｊａｍｕｎａ Ｐｒｅｓｓ，２００１）に記載されている。

電離層観測部２は、イオノゾンデ用アンテナ２１と、イオノゾンデ２２と、イオノゾンデ収集データ処理部２３とを備える。イオノゾンデ２２は、イオノゾンデ用アンテナ２１を介して複数の周波数の観測信号を電離層に送信し、送信した観測信号が電離層で反射された反射信号を受信して、観測信号の往復時間等のデータを収集する。

イオノゾンデ収集データ処理部２３は、イオノゾンデ２２で収集したデータから、電離層の高さ方向の電子密度分布情報（Ｅ層電子密度、Ｅ層ピーク電子密度高度、Ｆ１層電子密度、Ｆ１層ピーク電子密度高度、Ｆ２層電子密度、Ｆ２層ピーク電子密度高度、臨海（プラズマ）周波数など）を算出し、得られた算出結果をネットワーク７を介して信号処理部４およびデータサーバ部５へ送信する。また、イオノゾンデ収集データ処理部２３は、イオノゾンデ２２で収集したデータをネットワーク７を介してデータサーバ部５へ送信する。

インターネットデータ処理部３は、ＧＥＯＮＥＴ収集データ処理部３１と、中継部３２とを備える。ＧＥＯＮＥＴ収集データ処理部３１は、電離層に関連する情報（太陽黒点数、太陽フラックス強度Ｆ１０．７等）をインターネット６を介して取得する。

また、ＧＥＯＮＥＴ収集データ処理部３１は、一般に公開されているＧＰＳ観測データを、国土地理院が管理するＧＰＳ受信観測網（ＧＥＯＮＥＴ）や、国際的なＧＰＳ観測データを公開しているＩＧＳ（Ｉｎｔｅｒｎａｔｉｏｎａｌ ＧＰＳ Ｓｅｒｖｉｃｅ ｆｏｒ Ｇｅｏｄｙｎａｍｉｃｓ）からインターネット６経由で取得する。

中継部３２は、スイッチングハブあるいはルータにより構成され、ＧＥＯＮＥＴ収集データ処理部３１で取得した各種データをネットワーク７を介してデータサーバ部５へ送信する。インターネットデータ処理部３は外部とのつながりあるため、ＧＥＯＮＥＴ収集データ処理部３１および中継部３２は、ファイアオール機能を有するものとする。

信号処理部４は、複数のアンテナ４１ａ，４１ｂ，…，４１ｎからなるアレイアンテナ４１と、アンプ４２ａ，４２ｂ，…，４２ｎと、周波数変換部４３ａ，４３ｂ，…，４３ｎと、Ａ／Ｄ変換部４４ａ，４４ｂ，…，４４ｎと、演算部４５と、表示操作部４６とを備える。

アンテナ４１ａ，４１ｂ，…，４１ｎは、既存の放送局から送信され、電離層を通過して到来する放送信号（電波）を受信して受信信号を出力する。アンプ４２ａ，４２ｂ，…，４２ｎは、それぞれアンテナ４１ａ，４１ｂ，…，４１ｎからの信号を増幅し、周波数変換部４３ａ，４３ｂ，…，４３ｎは、アンプ４２ａ，４２ｂ，…，４２ｎで増幅された信号をベースバンドの信号に変換する。Ａ／Ｄ変換部４４ａ，４４ｂ，…，４４ｎは、周波数変換部４３ａ，４３ｂ，…，４３ｎからの信号をＡ／Ｄ変換してデジタル受信信号を出力する。

演算部４５は、各Ａ／Ｄ変換部４４ａ，４４ｂ，…，４４ｎからのデジタル受信信号から、ＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ Ｓｉｇｎａｌ Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法や独立成分分析(ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ)の手法などを用いて、既存の放送局からの電波の到来方位角および到来仰角を算出する。

また、演算部４５は、電離層電子密度分布モデルを用いて、既存の放送局からの電波の伝搬経路をレイトレーシング手法により算出し、算出した伝搬経路に基づいて、電波の到来方位角および到来仰角を算出する。なお、レイトレーシング手法の詳細については、例えば、“Ｔｈｅ Ｐｒｏｐａｇａｔｉｏｎ ｏｆ ｒａｄｉｏ ｗａｖｅｓ”（Ｋ．Ｇ．Ｂｕｄｅｎ，Ｃａｍｂｒｉｄｇｅ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，１９８８）に記載されている。

また、演算部４５は、衛星信号処理部１３で算出した各航法衛星の位置、および電離層電子密度分布モデルを用いて、各航法衛星から衛星信号受信アンテナ１１の設置位置までの衛星信号の通過経路上の電子密度を積算することにより、各通過経路の総電子数（ＴＥＣ）を算出する。

また、演算部４５は、Ａ／Ｄ変換部４４ａ〜４４ｎからのデジタル受信信号から算出した電波の到来方位角および到来仰角、レイトレーシング手法を用いて求めた電波の到来方位角および到来仰角、衛星信号処理部１３で算出した衛星信号の通過経路の総電子数、電離層電子密度分布モデルを用いて算出した衛星信号の通過経路の総電子数を用いて、電離層電子密度分布モデルに入力する各パラメータを初期値から変化させるための修正パラメータ変化量を算出し、算出した修正パラメータ変化量を用いて所望の位置の電子密度を算出する。

表示操作部４６は、ユーザによる入力信号の受付、また、演算部４５で得られた演算結果の表示を行う。

データサーバ部５は、記録処理部５１と、記録部５２とを備える。記録処理部５１は、衛星観測部１、電離層観測部２、インターネットデータ処理部３、信号処理部４で得られた各種演算結果や各種データをネットワーク７を介して受け取り、記録部５２に記録する処理を行う。

次に、本実施の形態に係る電離層電子密度算出装置の動作を説明する。

アンテナ４１ａ〜４１ｎは、図２に示すように既存の放送局から送信され、電離層を通過して到来する放送信号（電波）を受信して受信信号を出力する。アンプ４２ａ〜４２ｎは、それぞれアンテナ４１ａ〜４１ｎからの信号を増幅し、周波数変換部４３ａ〜４３ｎは、アンプ４２ａ〜４２ｎで増幅された信号をベースバンドの信号に変換する。そして、Ａ／Ｄ変換部４４ａ〜４４ｎは、周波数変換部４３ａ〜４３ｎからの信号をＡ／Ｄ変換してデジタル受信信号を演算部４５に出力する。

そして、演算部４５は、Ａ／Ｄ変換部４４ａ〜４４ｎからのデジタル受信信号から、ＭＵＳＩＣ法や独立成分分析の手法などを用いて、既存の放送局からの電波の到来方位角および到来仰角を算出する。以下、このようにアレイアンテナ４１で受信して得られた信号を用いて算出した到来方位角、到来仰角を適宜それぞれ観測方位角、観測仰角と呼ぶ。

また、演算部４５は、電離層電子密度分布モデルを使用し、放送局の位置から適当な送信方位角と送信仰角とを設定して、放送局からの電波の伝搬経路をレイトレーシング手法により算出する。なお、放送局の位置と、放送局からの送信周波数は既知であるとする。このとき算出した伝搬経路から得られる受信位置が、実際に受信を行ったアンテナ４１ａ〜４１ｎの設置位置になるよう、送信側の送信方位角と送信仰角とを調整する。調整後、レイトレーシング手法により求めた伝搬経路から、受信位置で受信される信号の到来方位角と到来仰角が求められる。

一方、衛星観測部１の衛星信号受信アンテナ１１は、図２に示すように、航法衛星９ａ，９ｂ，…から送信される複数周波数の衛星信号を観測できる場所に設置されており、衛星信号受信部１２は、各航法衛星９ａ，９ｂ，…から送信される複数周波数の衛星信号を、衛星信号受信アンテナ１１を介して受信する。

そして、衛星信号処理部１３は、衛星信号受信部１２で受信した複数周波数の衛星信号から、航法衛星の位置を算出するとともに、衛星信号の周波数による伝搬遅延量の違いを用いて、衛星信号の通過経路１０ａ，１０ｂ，…に存在する総電子数（ＴＥＣ）を求める。この際、インターネットで公開されている衛星信号の観測データがあれば、インターネットデータ処理部３でその観測データを取得し、取得した観測データも同様に使用して衛星信号の通過経路における総電子数を求める。以下、このように衛星信号受信アンテナ１１で受信した衛星信号またはインターネットデータ処理部３で取得した観測データを用いて算出した衛星信号の通過経路の総電子数を適宜観測総電子数と呼ぶ。

また、演算部４５は、衛星信号処理部１３で算出した各航法衛星の位置、および電離層電子密度分布モデルを用いて、各航法衛星９ａ，９ｂ，…から衛星信号受信アンテナ１１の設置位置（衛星観測位置）までの衛星信号の通過経路１０ａ，１０ｂ，…上の電子密度を積算することにより、各通過経路１０ａ，１０ｂ，…の総電子数（ＴＥＣ）を算出する。

電離層電子密度分布モデルでは、電離層のいくつかの特徴を示すパラメータを用いて電子密度の分布を算出している。その中で代表的なものを表１に示す。上述したような、電離層電子密度分布モデルを用いたレイトレーシング手法による到来方位角、到来仰角の算出や、電離層電子密度分布モデルを用いた衛星信号の通過経路の総電子数算出の際に、各計算箇所で、表１に示す各パラメータが電離層電子密度分布モデルから求められる。

表１において、Ｐは一般的な位置を示し、ＳＬＴ（Ｓｕｎ Ｌｏｃａｌ Ｔｉｍｅ）は、位置Ｐにおけるローカル時刻を示している。このＳＬＴは、太陽が南中した時を昼間の１２時とする時刻であり、経度１５°毎に１時間の時差が生じる。

ここで、表１に示した各パラメータを使用して、電離層電子密度分布モデルを用いたレイトレーシング手法による到来方位角、到来仰角、および衛星信号の通過経路における総電子数を求める関数をＦ（ｘ_１，ｘ_２，…，ｘ_７）とする。ここで、ｘ_１，ｘ_２，…，ｘ_７は、放送局位置、送信周波数、受信時間、送信仰角、送信方位角、電離層パラメータ（Ｅ層高度、Ｆ層高度、太陽黒点数、太陽フラックスＦ１０．７）など、既存の放送局から放送信号受信位置までの伝搬経路を計算するのに必要なパラメータと、電離層電子密度分布モデルを用いて衛星信号の通過経路における総電子数を計算するのに必要となる航法衛星位置と衛星観測位置などである。

この関数Ｆ（ｘ_１，ｘ_２，…，ｘ_７）の出力は、観測した放送局ごとに電離層電子密度分布モデルを使用してレイトレーシング手法により算出した到来方位角、到来仰角、および航法衛星ごとに電離層電子密度分布モデルを用いて算出した衛星信号の通過経路における総電子数である。

次いで、演算部４５は、放送局ごとに観測方位角とレイトレーシング手法により算出した到来方位角との差である第１の誤差、および観測仰角とレイトレーシング手法により算出した到来仰角との差である第２の誤差を算出する。また、演算部４５は、航法衛星ごとに、観測総電子数と電離層電子密度分布モデルを用いて算出した衛星信号の通過経路における総電子数との差である第３の誤差を算出する。

ここで、本実施の形態の電離層電子密度算出装置では、複数の放送局からの放送信号や、複数の航法衛星からの衛星信号を観測できるため、全体の誤差を表す指標としての誤差指標が必要になる。演算部４５は、算出した第１〜第３の誤差から、以下の（１）〜（４）のいずれかを誤差指標として算出する。（１）第１〜第３の誤差の絶対値の和（２）第１〜第３の誤差の絶対値の平方和（３）第１〜第３の誤差の絶対値の平方和の平方根（４）第１〜第３の誤差の分散の和
そして、演算部４５は、電離層電子密度分布モデルにより所望の位置の電子密度を算出するために用いる修正パラメータ変化量を算出する。以下、修正パラメータ変化量を算出する処理について説明する。

上述したように、電離層電子密度分布モデルを用いたレイトレーシング手法による到来方位角、到来仰角の算出や、電離層電子密度分布モデルを用いた衛星信号の通過経路における総電子数算出の際に、位置Ｐ、ローカル時刻ＳＬＴで表１の各パラメータが求められる。ここで算出された各パラメータの値を初期値とする。

まず、演算部４５は、各パラメータの初期値Ｘからの変化量ΔＸを、それぞれ乱数で決める。

上記（数式１）において、ΔＸ＝（Δｘ_１，Δｘ_２，…，Δｘ_７）^ｔが乱数により決められた値（変化量）である。このとき、この変化量ΔＸは、位置Ｐとローカル時刻ＳＬＴに依存していない。

そして、演算部４５は、この変化量ΔＸの分だけ初期値から変化させたパラメータＸ´を用いて、電離層電子密度分布モデルを用いたレイトレーシング手法により電波の到来方位角、到来仰角を算出するとともに、電離層電子密度分布モデルを用いて衛星信号の通過経路における総電子数を算出する。

そして、演算部４５は、ここで算出した到来方位角、到来仰角、および衛星信号の通過経路における総電子数と、観測方位角、観測仰角、および観測総電子数を用いて、以下の（数式２）〜（数式５）により、誤差指標ΔＦを算出する。ここでは、誤差指標として第１〜第３の誤差の分散の和を用いている。

ここで、σ_ＡＺ ^２は第１の誤差の分散、σ_ＥＬ ^２は第２の誤差の分散、σ_ＴＥＣ ^２は第３の誤差の分散、ＡＺ_{ｋ，ｍｏｄｅｌ}は電離層電子密度分布モデルを用いたレイトレーシング手法により算出した到来方位角、ＥＬ_{ｋ，ｍｏｄｅｌ}は電離層電子密度分布モデルを用いたレイトレーシング手法により算出した到来仰角、ＴＥＣ_{ｋ，ｍｏｄｅｌ}は電離層電子密度分布モデルを用いて算出した衛星信号の通過経路における総電子数、ＡＺ_{ｋ，ｏｂｓ}は観測方位角、ＥＬ_{ｋ，ｏｂｓ}は観測仰角、ＴＥＣ_{ｋ，ｏｂｓ}は観測総電子数、Ｍ_ｂｒａｄは１回の誤差指標作成の際に電離層電子密度分布モデルを用いたレイトレーシング手法による到来方位角、到来仰角の算出に使える放送局の数、Ｍ_ｓａｔは１回の誤差指標作成の際に電離層電子密度分布モデルを用いた総電子数の算出に使える航法衛星の数、Ｆ０は観測した放送信号および衛星信号から得られる算出結果（観測方位角、観測仰角、観測総電子数）を示す。

演算部４５は、算出した誤差指標ΔＦと予め設定した所定の閾値ε_{ｔｈｒｅｓｈ}とを比較し、誤差指標ΔＦが０以上で閾値ε_{ｔｈｒｅｓｈ}以下のとき、つまり、以下の（数式６）の条件を満たすとき、そのときのパラメータＸ´を初期値からの修正パラメータとする。

そして、このときの変化量ΔＸを修正パラメータ変化量として採用し、衛星信号や放送信号の通過していない任意の位置の電離層電子密度を算出する際にも、この修正パラメータ変化量を各パラメータの初期値に加算して用いる。なお、この際、電離層観測部２のイオノゾンデ収集データ処理部２３で算出した電離層の高さ方向の電子密度分布情報や、インターネットデータ処理部３で取得した電離層に関連する情報（太陽黒点数、太陽フラックス強度Ｆ１０．７等）を各パラメータの初期値として用いてもよい。

誤差指標ΔＦが閾値ε_{ｔｈｒｅｓｈ}より大きいとき、演算部４５は、新たな変化量ΔＸを乱数で発生して、変化量ΔＸの分だけ初期値から変化させた各パラメータを用いて、電離層電子密度分布モデルを用いたレイトレーシング手法により到来方位角、到来仰角を算出し、電離層電子密度分布モデルを用いて衛星信号の通過経路における総電子数を算出する。そして、演算部４５は、算出した到来方位角、到来仰角、および衛星信号の通過経路における総電子数を用いて誤差指標ΔＦを算出し、この誤差指標ΔＦを所定の閾値ε_{ｔｈｒｅｓｈ}と比較する。演算部４５は、この処理を誤差指標ΔＦが所定の閾値ε_{ｔｈｒｅｓｈ}以下になるまで繰り返し行う。

このように本実施の形態によれば、電離層電子密度分布モデルに入力するパラメータを変化させるための変化量を乱数で発生し、逐次誤差指標を算出するというモンテカルロ法的な手法により、誤差指標が所定の閾値以下になるときの変化量を求め、このときの変化量を修正パラメータ変化量として採用し、この修正パラメータ変化量を用いて任意の位置の電子密度を算出するので、現実の電離層の状態に即した電子密度を算出することができる。

（変形例１）
変形例１では、（数式１）の変化量ΔＸを位置Ｐとローカル時刻ＳＬＴに依存して変化させる。

変形例１における演算部４５の処理について説明する。まず、上記実施の形態と同様に、演算部４５は、表１の各パラメータの初期値からの変化量を、それぞれ乱数で決める。このとき、この変化量の適用できる位置を決める。例えば、放送信号を受信した位置Ｐｏとその時のローカル時刻ＳＬＴｏとする。レイトレーシング手法による到来方位角、到来仰角の算出や、電離層電子密度分布モデルを用いた衛星信号の通過経路における総電子数算出の際に、各算出位置に適用できるよう、変化量を算出した位置の違い、ローカル時刻を考慮し、算出する位置における変化量を補正する。補正の方法としては、以下の（数式７）に示す日変化を計算する関数ＨＰＯＬを使用する。

ここで、Ｙｕは一日の中で変化する量の最大値、Ｙｂは一日の中で変化する量の最小値、ＳＲはローカル時刻の日の出時刻[時]、ＳＳはローカル時刻の日没時刻[時]、ＳＬＴは計算を行う位置のローカル時刻である。

なお、上記（数式７）の日変化を計算する関数ＨＰＯＬの詳細については、例えば、非特許文献１に記載されている。

次いで、演算部４５は、乱数で発生した変化量ΔＸが適用される位置Ｐｏにおける日の出時刻ＳＲ_０と日没時刻ＳＳ_０とから、以下の（数式８）により、真昼のローカル時刻Ｔ_ｄａｙと深夜時刻Ｔ_{ｎｉｇｈｔ}とを求める。

次に、（数式７）のＹｕ，Ｙｂの近似値として、Ｙｕを日中Ｔ_ｄａｙのパラメータの値とし、Ｙｂを深夜Ｔ_{ｎｉｇｈｔ}のパラメータの値とする。そして、演算部４５は、表１に示す７つのパラメータそれぞれを、電離層電子密度分布モデルを使用し、位置Ｐｏにおいて求める。このときに生成される各パラメータの値をそれぞれの、ＹｕとＹｂに相当する値とする。

同様に、演算部４５は、所望の位置Ｐｋにおける日の出時刻、日没時刻から、位置Ｐｋにおける（数式７）のＹｕ，Ｙｂの近似値を、電離層電子密度分布モデルにより求める。そして、演算部４５は、位置Ｐｏと位置ＰｋのＹｕ，Ｙｂの値と、位置Ｐｏにおける変化量とを用いて、所望の位置Ｐｋにおける変化量を次の（数式９）により求める。演算部４５は、この処理を表１に示す各パラメータに対して行う。

ここで、Δｘ_ｊ ^Ｐｋは位置Ｐｋにおける各パラメータの変化量、Δｘ_ｊ ^Ｐｏは位置Ｐｏにおける各パラメータの変化量、（Ｙｕ−Ｙｂ）_ｊ ^Ｐｋは位置ＰｋにおけるＹｕとＹｂとの差、（Ｙｕ−Ｙｂ）_ｊ ^Ｐｏは位置ＰｏにおけるＹｕとＹｂとの差、ｊは表１に示す各パラメータの識別子である。

このようにして算出した位置Ｐｋにおける変化量を用いて、各パラメータを補正し、この補正した各パラメータを電離層電子密度分布モデルへ代入して、位置Ｐｋにおける電子密度を求める。この操作を、レイトレーシング手法による到来方位角、到来仰角の算出や、電離層電子密度分布モデルを用いた衛星信号の通過経路における総電子数算出の際に、電子密度を求めたい位置ごとに行う。

そして、上記実施の形態と同様に、演算部４５は、算出した到来方位角、到来仰角、および総電子数と、観測方位角、観測仰角、および観測総電子数を用いて、上述の（数式２）〜（数式５）により、誤差指標ΔＦを算出し、この誤差指標ΔＦを所定の閾値ε_{ｔｈｒｅｓｈ}と比較する。

演算部４５は、この誤差指標ΔＦが閾値ε_{ｔｈｒｅｓｈ}以下になるまで、変化量ΔＸを乱数で発生させ、以下同様の処理を行う。そして、演算部４５は、誤差指標ΔＦが閾値ε_{ｔｈｒｅｓｈ}以下になった際の変化量ΔＸを（数式９）に従う日変化を考慮して修正し、この修正した変化量を、任意の位置における表１に示す各パラメータに適用し、電子密度を算出する。

このように変形例１では、各パラメータの変化量ΔＸを、各パラメータの日変化を考慮して修正するので、より正確に現実の電離層の状態に即した電子密度を算出することができる。

（変形例２）
上記実施の形態では、単純に乱数を発生させてパラメータの変化量を決め、何度か繰り返し試行することにより適切な変化量の値を見つけ出す。しかしながら、この方法は非常に強力ではあるが、繰り返し試行の回数が多くなり、処理に長時間を要するおそれがある。

変形例２では、以下に述べる最小二乗法を用いて誤差指標が小さくなる方向にパラメータの変化量を調整することにより、試行回数を少なくする。

変形例２における演算部４５の処理について説明する。まず、演算部４５は、各パラメータの初期値Ｘからの変化量ΔＸを、数回から数十回程度（パラメータ数以上の回数Ｍ）、乱数により決める。このとき、まったく初期値を指定せずに乱数により変化量ΔＸを発生させてもよいし、パラメータごとに決めた範囲を、ステップごとに変化させて変化量ΔＸを発生させてもよい。なお、このとき、変化量ΔＸは、位置Ｐとローカル時刻ＳＬＴに依存していない。

そして、演算部４５は、変化量ΔＸを変化させるごとに、変化量ΔＸの分だけ初期値から変化させたパラメータを用いて、電離層電子密度分布モデルを用いたレイトレーシング手法により電波の到来方位角、到来仰角を算出するとともに、電離層電子密度分布モデルを用いて衛星信号の通過経路における総電子数を算出する。

そして、演算部４５は、算出した到来方位角、到来仰角、および衛星信号の通過経路における総電子数と、観測方位角、観測仰角、および観測総電子数とを用いて、上述の（数式２）〜（数式５）により、誤差指標ΔＦを算出する。これを発生した変化量ΔＸの個数分（Ｍ回）だけ行う。

誤差指標ΔＦはパラメータ（ｘ_１，ｘ_２，…，ｘ_７）の関数であり、変化量ΔＸを加算する前のパラメータの初期値（ｘ_１０，ｘ_２０，…，ｘ_７０）近傍でテーラー展開すると、以下の（数式１０）のようになる。

ここで、ΔＦ０は、変化量ΔＸを加算する前の各パラメータの初期値に基づいて算出した到来方位角、到来仰角、衛星信号の通過経路における総電子数と、観測方位角、観測仰角、および観測総電子数とを用いて算出した誤差指標である。変化量ΔＸを加算する前の各パラメータが適切であれば、ΔＦ０はゼロに近い値になっているはずである。

レイトレーシング手法により電波の到来方位角、到来仰角を算出し、電離層電子密度分布モデルを用いて衛星信号の通過経路における総電子数を算出し、算出した到来方位角、到来仰角、および衛星信号の通過経路における総電子数を用いて誤差指標ΔＦを算出する処理をＭ回試行することにより、（数式１０）は複数できる。

そして、演算部４５は、変化量ΔＸを加算する前の各パラメータの初期値近傍で、各パラメータを変化させたとき、誤差指標ΔＦがどのように変化するかを決める変数Ｇを、以下の（数式１２）のように、最小二乗法により求める。

これにより、各パラメータの初期値を使った際に、各パラメータを変化させた場合、どれだけ誤差指標ΔＦが変化するかという変化率を得たことになる。また、この変数Ｇは、（ｘ_１，ｘ_２，…，ｘ_７）が張る空間の誤差局面のグラジエントになっている。

ここで、求められた変数Ｇを用いて、誤差指標ΔＦを表すと、次の（数式１３）を得る。

ここで、閾値ε_threshはあらかじめ設定した値である。上記（数式１３）の条件を満たすには、誤差指標ΔＦが最小になるような変化量ΔＸ＝（Δｘ_１，Δｘ_２，…，Δｘ_７）^ｔを決定すればよい。

変数ｇ_ｋは、誤差指標ΔＦのグラジエントであり、変数ｇ_ｋの符号により、Δｘ_ｋが増加または減少した際にΔＦが減少するようにすればよいので、次の（数式１４）のような関係が成り立つ。

そして、演算部４５は、上記（数式１４）の条件を満足する変化量ΔＸ＝（Δｘ_１，Δｘ_２，…，Δｘ_７）^ｔを乱数により決める。このとき、まったく初期値を指定せずに乱数により変化量ΔＸを発生させてもよいし、パラメータごとに決めた範囲を、ステップごとに変化させて変化量ΔＸを発生させてもよい。

そして、上記実施の形態と同様に、演算部４５は、変化量ΔＸの分だけ初期値から変化させたパラメータを用いて、電離層電子密度分布モデルを用いたレイトレーシング手法により電波の到来方位角、到来仰角を算出するとともに、電離層電子密度分布モデルを用いて衛星信号の通過経路における総電子数を算出する。

そして、演算部４５は、算出した到来方位角、到来仰角、および衛星信号の通過経路における総電子数と、観測方位角、観測仰角、および観測総電子数とを用いて、上述の（数式２）〜（数式５）により、誤差指標ΔＦを算出し、この誤差指標ΔＦを所定の閾値ε_{ｔｈｒｅｓｈ}と比較する。

演算部４５は、この誤差指標ΔＦが閾値ε_{ｔｈｒｅｓｈ}以下になるまで、ΔＸを乱数で発生させ、以下同様の処理を行う。そして、演算部４５は、誤差指標ΔＦが閾値ε_{ｔｈｒｅｓｈ}以下になった際の変化量ΔＸを修正パラメータ変化量として採用し、衛星信号や放送信号の通過していない任意の位置の電離層電子密度を算出する際にも、この修正パラメータ変化量を各パラメータの初期値に加算して用いる。

このように変形例２では、最小二乗法を用いて誤差指標が小さくなる方向にパラメータの変化量を調整することで、修正パラメータ変化量を決定するまでの試行回数を低減し、処理時間の短縮を図ることができる。

また、変形例２においても、変形例１と同様に、レイトレーシング手法による到来方位角、到来仰角の算出や、電離層電子密度分布モデルを用いた衛星信号の通過経路における総電子数算出の際に、変化量ΔＸを（数式９）に従う日変化を考慮して修正するようにしてもよい。このようにすることで、より正確に現実の電離層の状態に即した電子密度を算出することができる。

なお、上記実施の形態は、あくまでも本発明の説明のためのものであり、本発明の範囲を制限するものではない。したがって、当業者であれば、これらの各要素または全要素を含んだ各種の実施形態を採用することが可能であるが、これらの実施形態も本発明の範囲に含まれる。

本発明の実施の形態に係る電離層電子密度算出装置の構成を示すブロック図である。 衛星信号の通過経路における総電子数の算出について説明するための図である。

符号の説明

１ 衛星観測部
２ 電離層観測部
３ インターネットデータ処理部
４ 信号処理部
５ データサーバ部
１１ 衛星信号受信アンテナ
１２ 衛星信号受信部
１３ 衛星信号処理部
２１ イオノゾンデ用アンテナ
２２ イオノゾンデ
２３ イオノゾンデ収集データ処理部
３１ ＧＥＯＮＥＴ収集データ処理部
３２ 中継部
４１ アレイアンテナ
４１ａ，４１ｂ，…，４１ｎ アンテナ
４２ａ，４２ｂ，…，４２ｎ アンプ
４３ａ，４３ｂ，…，４３ｎ 周波数変換部
４４ａ，４４ｂ，…，４４ｎ Ａ／Ｄ変換部
４５ 演算部
４６ 表示操作部
５１ 記録処理部
５２ 記録部

Claims (8)

1. 送信される位置および送信周波数が既知であり、電離層を通過して到来する電波を受信するアレイアンテナと、
   このアレイアンテナからの出力信号から前記電波の到来方位角、到来仰角を第１の到来方位角、第１の到来仰角として算出する第１の方向算出手段と、
   電離層の電子密度分布を示す電離層電子密度分布モデルに入力する複数のパラメータをその初期値から変化させるための変化量を乱数により生成する変化量生成手段と、
   この変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記アレイアンテナで受信した前記電波の伝搬経路を算出し、この伝搬経路に基づいて前記電波の到来方位角、到来仰角を第２の到来方位角、第２の到来仰角として算出する第２の方向算出手段と、
   航法衛星から送信される複数周波数の衛星信号を受信する衛星信号受信手段と、
   この衛星信号受信手段で受信した前記複数周波数の衛星信号の電離層での周波数による伝搬遅延量の差を用いて、前記衛星信号の通過経路における総電子数を第１の総電子数として算出する第１の総電子数算出手段と、
   前記変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記衛星信号の通過経路における総電子数を第２の総電子数として算出する第２の総電子数算出手段と、
   前記第１の到来方位角と前記第２の到来方位角との差である第１の誤差を算出し、前記第１の到来仰角と前記第２の到来仰角との差である第２の誤差を算出し、前記第１の総電子数と前記第２の総電子数との差である第３の誤差を算出するとともに、前記第１〜第３の誤差を合わせた値の大きさの指標である誤差指標を算出する誤差指標算出手段と、
   この誤差指標算出手段で算出した前記誤差指標が所定の閾値以下であるか否かを判断し、前記誤差指標が前記所定の閾値以下のとき、当該誤差指標の算出に用いた前記第２の到来方位角、前記第２の到来仰角、および前記第２の総電子数の算出に用いた変化量の分だけ所望の位置における初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記所望の位置の電子密度を算出し、前記誤差指標が前記所定の閾値より大きいときは、前記変化量生成手段に新たな変化量を生成させる電子密度算出手段と
   を備えることを特徴とする電離層電子密度算出装置。
2. 前記変化量生成手段は、前記各パラメータの時刻による変化に応じて前記変化量を変化させることを特徴とする請求項１に記載の電離層電子密度算出装置。
3. 前記誤差指標は、前記第１〜第３の誤差の絶対値の和、前記第１〜第３の誤差の絶対値の平方和、前記第１〜第３の誤差の絶対値の平方和の平方根、および前記第１〜第３の誤差の分散の和のうちのいずれかであることを特徴とする請求項１または２に記載の電離層電子密度算出装置。
4. 送信される位置および送信周波数が既知であり、電離層を通過して到来する電波を受信するアレイアンテナと、
   このアレイアンテナからの出力信号から前記電波の到来方位角、到来仰角を第１の到来方位角、第１の到来仰角として算出する第１の方向算出手段と、
   電離層の電子密度分布を示す電離層電子密度分布モデルに入力する複数のパラメータをその初期値から変化させるための変化量を、前記複数のパラメータの数より多い所定の個数だけ乱数により生成する第１の変化量生成手段と、
   この第１の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記アレイアンテナで受信した前記電波の伝搬経路を算出し、この伝搬経路に基づいて前記電波の到来方位角、到来仰角を第２の到来方位角、第２の到来仰角として算出する処理を、前記所定の個数の変化量のそれぞれについて行う第２の方向算出手段と、
   航法衛星から送信される複数周波数の衛星信号を受信する衛星信号受信手段と、
   この衛星信号受信手段で受信した前記複数周波数の衛星信号の電離層での周波数による伝搬遅延量の差を用いて、前記衛星信号の通過経路における総電子数を第１の総電子数として算出する第１の総電子数算出手段と、
   前記第１の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記衛星信号の通過経路における総電子数を第２の総電子数として算出する処理を、前記所定の個数の変化量のそれぞれについて行う第２の総電子数算出手段と、
   前記第１の到来方位角と前記第２の到来方位角との差である第１の誤差を算出し、前記第１の到来仰角と前記第２の到来仰角との差である第２の誤差を算出し、前記第１の総電子数と前記第２の総電子数との差である第３の誤差を算出するとともに、前記第１〜第３の誤差を合わせた値の大きさの指標である誤差指標を算出する処理を、前記所定の個数の変化量のそれぞれについて行う第１の誤差指標算出手段と、
   この第１の誤差指標算出手段において前記所定の個数の変化量のそれぞれに対応して算出された第１の誤差指標に基づいて、前記複数のパラメータを前記変化量に応じて変化させたときに第１の誤差指標がどのように変化するかを決める変数を、最小二乗法を用いて算出する変数算出手段と、
   この変数算出手段で算出した前記変数を用いて、前記各パラメータをその初期値から変化させるための変化量を乱数により生成する第２の変化量生成手段と、
   この第２の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記アレイアンテナで受信した前記電波の伝搬経路を算出し、この伝搬経路に基づいて前記電波の到来方位角、到来仰角を第３の到来方位角、第３の到来仰角として算出する第３の方向算出手段と、
   前記第２の変化量生成手段で生成した変化量の分だけ初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記衛星信号の通過経路における総電子数を第３の総電子数として算出する第３の総電子数算出手段と、
   前記第１の到来方位角と前記第３の到来方位角との差である第４の誤差を算出し、前記第１の到来仰角と前記第３の到来仰角との差である第５の誤差を算出し、前記第１の総電子数と前記第３の総電子数との差である第６の誤差を算出するとともに、前記第４〜第６の誤差を合わせた値の大きさの指標である第２の誤差指標を算出する第２の誤差指標算出手段と、
   この第２の誤差指標算出手段で算出した前記第２の誤差指標が所定の閾値以下であるか否かを判断し、前記第２の誤差指標が前記所定の閾値以下のとき、当該第２の誤差指標の算出に用いた前記第３の到来方位角、前記第３の到来仰角、および前記第３の総電子数の算出に用いた変化量の分だけ所望の位置における初期値から変化した前記各パラメータ、および前記電離層電子密度分布モデルを用いて、前記所望の位置の電子密度を算出し、前記第２の誤差指標が前記所定の閾値より大きいときは、第２の変化量生成手段に新たな変化量を生成させる電子密度算出手段と
   を備えることを特徴とする電離層電子密度算出装置。
5. 前記第１および第２の変化量生成手段は、前記各パラメータの時刻による変化に応じて前記変化量を変化させることを特徴とする請求項４に記載の電離層電子密度算出装置。
6. 前記第１の誤差指標は、前記第１〜第３の誤差の絶対値の和、前記第１〜第３の誤差の絶対値の平方和、前記第１〜第３の誤差の絶対値の平方和の平方根、および前記第１〜第３の誤差の分散の和のうちのいずれかであり、前記第２の誤差指標は、前記第４〜第６の誤差の絶対値の和、前記第４〜第６の誤差の絶対値の平方和、前記第４〜第６の誤差の絶対値の平方和の平方根、および前記第４〜第６の誤差の分散の和のうちのいずれかであることを特徴とする請求項４または５に記載の電離層電子密度算出装置。
7. 複数周波数の観測信号を電離層に送信し、前記観測信号が電離層で反射された反射信号を受信して、電離層の高さ方向の電子密度分布情報を取得する電離層観測手段を備え、
   前記電子密度算出手段は、前記電離層観測手段で取得した前記電子密度分布情報を前記各パラメータの所望の位置における初期値の一部として用いることを特徴とする請求項１〜６のいずれか１項に記載の電離層電子密度算出装置。
8. 複数の観測地点で複数の航法衛星を観測して得られる観測データを外部から取得するデータ取得手段を備え、
   前記第１の総電子数算出手段は、前記データ取得手段で取得した前記観測データを用いて前記各航法衛星から送信される衛星信号の通過経路における総電子数を前記第１の総電子数として算出することを特徴とする請求項１〜７のいずれか１項に記載の電離層電子密度算出装置。

Images