"複数のアンテナと、前記複数のアンテナから得られる複数のアンテナ受信信号に基づき推定した到来波数と到来信号を用いてレイトレーシングによる電波伝搬経路推定を行って到来方向を求め、該求めた到来方向を初期値として到来方位および仰角を測定する信号処理装置と、を備えることを特徴とする方位仰角測定システム。前記信号処理装置は、アンテナ受信信号に基づき、IRIモデルまたは該IRIモデルをGPS信号で修正したモデルにより電離層電子密度を推定して電波伝搬経路推定を行うことを特徴とする請求項1記載の方位仰角測定システム。上述したAkaikesandRissanenscriteriaの方法を用いる前提条件は、信号が独立性を有し、ノイズが正規分布していることであるが、実際には、独立性がないコヒーレントマルチパス信号が存在しているため、コヒーレントマルチパス信号を含む受信信号の信号数は不明確である。また、コヒーレント信号の到来方向、特に仰角については、初期値として適切な値を設定することは難しい。その結果、到来信号の方位および仰角を精密に測定できないという問題がある。本発明が解決しようとする課題は、到来信号の方位および仰角を精密に測定できる方位仰角測定システムを提供することにある。実施形態に係る方位仰角測定システムによれば、複数のアンテナと、複数のアンテナから得られる複数のアンテナ受信信号に基づき推定した到来波数と到来信号を用いてレイトレーシングによる電波伝搬経路推定を行って到来方向を求め、該求めた到来方向を初期値として到来方位および仰角を測定する信号処理装置を備えることを特徴とする。以下、本発明の実施例について、図面を参照しながら詳細に説明する。図1は、第1の実施形態に係る方位仰角測定システムの構成を示すブロック図である。方位仰角測定システムは、信号処理系1、衛星信号処理系2、電離層観測系3、インターネットデータ処理系4およびデータサーバ系5を備えている。信号処理系1は、アレイ信号処理を行う部分であり、p個pは2以上の整数の通信電波受信用のアンテナ111~11p、p個のアンプ121~12p、p個の周波数変換部131~13p、p個のデジタイザ141~14p、信号処理装置15および表示部16を備えている。p個のアンテナ111~11pは、外部から到来する電波を受信してアナログの電気信号に変換し、p個のアンプ121~12pにそれぞれ送る。p個のアンプ121~12pは、p個のアンテナ111~11pからそれぞれ送られてくる電気信号を増幅し、p個の周波数変換部131~13pにそれぞれ送る。p個の周波数変換部131~13pは、p個のアンプ121~12pからそれぞれ送られてくる増幅された電気信号の周波数を変換し、p個のデジタイザ141~14pにそれぞれ送る。p個のデジタイザ141~14pは、p個の周波数変換部131~13pからそれぞれ送られてくる周波数変換されたアナログ信号をデジタル信号に変換し、アンテナ受信信号として信号処理装置15に送る。信号処理装置15は、この方位仰角測定システムの主要部であり、到来信号の到来方向測定を実施する。この信号処理装置15は、到来方向推定のために、p個のデジタイザ141~14pからそれぞれ送られてくるp個のアンテナ受信信号に基づき電離層電子密度を推定し、電波伝搬経路推定を行う。電離層電子密度の推定には、IRIInternationalReferenceIonosphereモデルやこのIRIモデルをGPSGlobalPositioningSystem衛星等からのGPS信号で修正したTECTotalElectronContentsを用いることができる。また、信号処理装置15は、電離層観測系3やインターネットデータ観測系4から送られてくる電離層実測データを使用した電離層電子密度の推定も実行する。この信号処理装置15で行う到来方向推定、電離層電子密度推定および電波伝搬経路推定は、1つのCPUCentralProcessingUnitで行うこともできるし、複数のCPUで並行して行うこともできる。電波伝搬経路推定の場合、決められた計算開始場所と到来方向を、複数のCPUで同時に計算させ、到来方向測定の結果を総合するように構成することができる。表示部16は、ディスプレイを含むパーソナルコンピュータから構成されており、信号処理装置15における処理結果を表示する。この表示部16は、パーソナルコンピュータを複数用意して並行に処理を行うように構成することができる。衛星信号処理系2は、衛星信号受信アンテナ21、衛星信号受信器22および衛星信号処理装置23を備えている。衛星信号受信アンテナ21は、GPS衛星またはGalileo衛星等から送られてくる複数の電波を受信し、電気信号に変換して衛星信号受信器22に送る。衛星信号受信器22は、衛星信号受信アンテナ21から送られてくる複数の電気信号の増幅等を行い、衛星信号処理装置23に送る。衛星信号処理装置23は、衛星信号受信器22から送られてくる複数の信号から擬似距離、位相および軌道暦を取得して信号処理系1およびデータサーバ系5へ送る。電離層観測系3は、イオノゾンデ用アンテナ31、イオノゾンデ32およびイオノゾンデ収集データ処理装置33を備えている。イオノゾンデ用アンテナ31は、パルス電波を上空に放射し、電離層で反射された電波を受信して電気信号に変換し、イオノゾンデ32に送る。イオノゾンデ32は、イオノゾンデ用アンテナ31から送られてきた電気信号から、パルス電波の放射から反射波を受信するまでの遅延時間を求めて電離気体の分布を表すイオノグラムを作成し、イオノゾンデ収集データ処理装置33に送る。イオノゾンデ収集データ処理装置33は、イオノゾンデ32から送られてくるイオノグラムから電離層の臨海周波数等を算出し、信号処理系1およびデータサーバ系5に送る。インターネットデータ処理系4は、ネットワーク接続機器41およびGEONET収集データ処理装置42を備えている。ネットワーク接続機器41は、スイッチングハブまたはルータから構成されており、このインターネットデータ処理系4を他の系統に接続する。GEONET収集データ処理装置42は、例えば、ディスプレイを含むパーソナルコンピュータから構成されており、電離層観測情報または衛星観測情報等を入手し、信号処理系1およびデータサーバ系5に送る。また、GEONET収集データ処理装置42は、外部のインターネット網に接続されている。データサーバ系5は、ディスプレイを含むパーソナルコンピュータからなるデータ収集サーバ51と記録装置52を備えている。データ収集サーバ51は、上述した信号処理系1、衛星信号処理系2、電離層観測系3およびインターネットデータ処理系4で処理されたデータを受信して記録装置52に記録するとともに、記録装置52に記録されているデータを読み出して信号処理系1、衛星信号処理系2、電離層観測系3およびインターネットデータ処理系4に送信し、また、記録装置52に記録されているデータを削除する。次に、上記のように構成される第1の実施形態に係る方位仰角測定システムの動作を説明する。図2は、第1の実施形態に係る方位仰角測定システムにおいて実行される到来方向初期値推定処理を示すフローチャートである。なお、以下に示す処理は、主に信号処理系1において実行される。到来方向初期値推定処理では、まず、受信信号共分散行列作成固有値固有ベクトル計算が実行されるステップS11。以下、詳細に説明する。信号処理系1のアンテナ111~11pで受信された信号sは、下記1式により表される。ここで、ここでE[]は、期待値を計算することを意味し、Rは観測された信号の共分散行列、Esは信号固有ベクトル、Λは固有値を対角要素とする行列、σはノイズの標準偏差、Iは単位行列を表しており、Hはエルミート転置を示している。m1~mpは固有値を示している。このうち、m1~mMは信号に対応し、mM1~mpはノイズに対応した固有値である。次いで、到来信号数の推定が行われるステップS12。この到来信号数の推定は、統計的な性質を利用したAkaikesAICまたはRissanensMDLによる方法非特許文献3参照、あるいは、7式の固有値が予め設定した閾値より大きい数として設定する。AICまたはMDLによる方法は、下記8式または9式のkを変化させ、AICまたはMDLの値が最小となった際のkを推定信号数とする方法である。ここで、m1>m2>...>mMk=0,1,,p-1である。次いで、到来方向の推定が行われるステップS13。到来方向の推定方式には、ビームフォーマ法、CAPON法、MUSIC法またはESPRIT法非特許文献4参照などが存在し、何れの方法を用いてもよい。以下では、この中のMUSIC法について説明する。ここで、En:ノイズ固有ベクトルθφ:ピークサーチ仰角方位角である。MUSIC法は、共分散行列の雑音固有値と信号に対応する固有値が直行する性質を利用する方法である。到来方向をサーチする方位仰角を変化させてステアリングベクトルを算出し、10式がピークとなる方向を検出する。次いで、ステップS12で推定された到来信号の信号数が複数であるかどうかが調べられるステップS14。このステップS14において、信号数が複数であることが判断されると、信号抽出が行われるステップS15。信号抽出の方法には、MSN法、DCMP法、ICA独立成分分析法、ゼロフォーシングZF法、MMSE法非特許文献4参照などが存在し、何れの方法を用いてもよい。以下では、MSN法とZF法について説明する。これらの方法では、抽出したい信号の到来方向に対するステアリングベクトルaθk,φkk=1~sと共分散行列が使用される。11式のmkおよびσは、対応する固有値およびノイズ標準偏差をそれぞれ示している。MSN法による抽出方法を11式および12式に示す。11式および12式のEは、固有値mkの固有ベクトルを示しており、12式のs0kは計算された初期値の受信信号を示し、Xはアンテナ受信信号を示している。ZF法による抽出方法を13式に示す。以上の信号抽出処理が終了すると、ステップS16の処理に進む。なお、上記ステップS14において、信号数が複数でないことが判断された場合は、ステップS15の処理はスキップされ、ステップS16の処理に進む。ステップS16においては、電波伝搬経路推定および減衰量推定が行われる。電波伝搬経路推定では、抽出されたk信号の伝搬経路をレイトレーシングにより推定する処理が行われる。なお、レイトレーシングによる電波伝搬経路推定については、前田憲一,後藤三男:電波伝播,岩波全書1953年2月、KGBudenThePropagationofradiowavesCambridgeUniversityPress1988またはIwaneKimuraEffectsofIonsonWhistlerModeRayTracingRadioScienceVolINo3269283March1966などに説明されている。レイトレーシングによる電波伝搬経路推定では、受信した信号のホップhop数は不明であるため、ホップ数分のホップ数伝搬経路推定が行われる。このとき、伝搬モードとして、左偏波Oモードおよび右偏波Xモードの両方、または、何れか一方について推定が行われる図3、図4参照。推定された伝搬経路の中で、地表に到達した位置計算位置から受信位置まで逆に電波伝搬経路推定を実行する。図4は、計算位置2から、1ホップ、2ホップ、3ホップ、の逆計算を示している。このときも、OモードおよびXモードの両方、または、何れか一方について逆計算が実施される。上述した電波伝搬経路推定を行うことにより、伝搬過程で生じる減衰量を推定できる前田憲一,後藤三男:電波伝播,岩波全書1953年2月参照。下記14式~16式は、近似的な減衰量計算方法を示している。ここで、E:受信電界強度E0:送信源の電界強度R1:地表-電離層間の伝搬距離の総和R2:電離層通過距離の総和Γ:電離層による減衰量q:電子の電荷c:光速度me:電子質量ε0:真空の誘電率f:電波の周波数STEP:計算ステップ総数ν0:電離層最下部における平均衝突回数α:高さに対する変化率zk:ステップkにおける高さNk:高さzkにおける電子密度nk:高さzkにおける屈折率Δl:ステップkとステップk1の距離図5は、上述した電波伝搬経路推定および減衰量推定の処理の詳細を示すフローチャートである。電波伝搬経路推定および減衰量推定の処理では、まず、逆計算仰角初期設定が行われるステップS21。次いで、計算位置を開始位置としたレイトレーシングが実施されるステップS22。次いで、電離層貫通があったかどうかが調べられるステップS23。このステップS23において、電離層貫通があったことが判断されると、ステップS28の処理に進む。一方、ステップS23において、電離層貫通がなかったことが判断されると、次いで、計算位置と受信位置との差の絶対値がεより小さいかどうかが調べられるステップS24。このステップS24において、計算位置と受信位置との差の絶対値がεより小さくないことが判断されると、逆計算仰角が変更されステップS25、その後、ステップS22の処理に戻って上述した処理が繰り返される。一方、上記ステップS24において、計算位置と受信位置との差の絶対値がεより小さいことが判断されると、計算された到来仰角方位角が保存されるステップS26。次いで、減衰量が計算されて保存されるステップS27。その後、ステップS28の処理に進む。ステップS28においては、OモードおよびXモードが完了したかどうかが調べられる。このステップS28において、OモードおよびXモードが完了していないことが判断されると、OモードからXモードまたはXモードからOモードへのモード変更が行われステップS29、その後、ステップS22の処理に戻って上述した処理が繰り返される。上記ステップS28において、OモードおよびXモードが完了したことが判断されると、次いで、ホップが最大ホップ数より大きくなったがどうかが調べられるステップS30。このステップS30において、ホップが最大ホップ数より大きくなっていないことが判断されると、ホップ数がインクリメント1されるステップS31。次いで、逆計算仰角初期設定が行われステップS32、その後、ステップS22の処理に戻って上述した処理が繰り返される。上記ステップS30において、ホップが最大ホップ数より大きくなったことが判断されると、次いで、全部の逆計算開始位置に対する処理が完了したかどうかが調べられるステップS33。このステップS33において、全部の逆計算開始位置に対する処理が完了していないことが判断されると、ステップS22の処理に戻って上述した処理が繰り返される。一方、ステップS33において、全部の逆計算開始位置に対する処理が完了したことが判断されると、電波伝搬経路推定および減衰量推定の処理は終了する。図2に示すフローチャートのステップS16における電波伝搬経路推定および減衰量推定の処理が終了すると、次いで、到来方向初期値設定が行われるステップS17。このステップS17では、上記ステップS16で求められた到来信号の到来方向θu,φu=θhk,φhkh:逆計算ホップ数、k:計算位置番号u=1~mと減衰量が使用される。ここで、mは、求められたマルチパスの数である。12式で求められた信号s0は、m個のコヒーレントマルチパスの和として表される。各アンテナ素子で受信した信号x0uは、信号uの到来方向に対応するステアリングベクトルを掛けて求められる。以上の処理の後、到来方位仰角精密測定が行われる。すなわち、ステップS17で求められた到来信号s0の到来方向推定の初期値を基に、さらに詳細な到来方位仰角推定が行われる。推定方法の一例としてSAGE法を使用した例を以下に示す。複数のマルチパスから構成される信号s0は、24式で表される。ここで、25式のkは繰り返し回数を示している。SAGE法の処理は、EExpectation-StepとMMaximization-Stepといった大きく2つのステップから構成される。まずE-Stepが実施され、次にM-Stepが実施され、これらのステップが、予め設定された評価値が閾値以下になるまで繰り返される。[E-Step]まず、複数信号から構成されている信号のうち、信号uの各アンテナ受信信号が26式により更新される。[M-Step]次いで、E-Stepで更新された信号uの共分散行列が計算される。ここで、27式のE[]は、括弧内の平均をとる操作である。この共分散行列を利用して、ビームフォーミングアルゴリズムにより到来方位が求められる非特許文献2参照。次いで、M-Stepで計算されたu信号の到来方向を用いてu信号のステアリングベクトルのみが更新される。次いで、更新された、アレーマニフォールド混合行列Aが生成される。次いで、これを用いてZF法により信号の分離が行われる。このとき、29式の括弧内のεIは、行列が非正則の場合に正則化するためのパラメータである。εがゼロの場合、正則でない場合に、十分小さい値0.001~0.00001が設定される。繰り返し処理の終了条件は、E-Stepの変化が少なくなった場合である。収束条件の一例を以下に示す。30式のεは、予め設定された閾値である。以上の処理により28式で求められた方位角および仰角が精密方位仰角とされる。このように、第1の実施形態に係る方位仰角測定システムによれば、レイトレーシングによる電波伝搬経路推定を用いることにより、到来波数と到来信号を推定し、SAGEアルゴリズムに代表されるアルゴリズムに、初期設定を行うことができる。その結果、ランダムに初期信号を設定した場合よりも効率的に処理を行うことができ、より実際の信号に近い結果を得ることができる。その結果、到来信号の方位および仰角を精密に測定できる。以上のように、いくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。"