(57)【要約】【課題】  ２次元で電波到来方向を推定でき、しかもコヒーレントやインコヒーレントの区別なく従来技術に比較して高精度で電波到来方向の推定を行う。【解決手段】  方形の電子制御導波器アレーアンテナ装置の各可変リアクタンス素子による放射パターンの各方向ベクトルに対する重み付けを表す等価的なベクトルからなる等価ウエイト行列を校正し、複数の受信信号間の相関を表す相関行列Ｒｙの固有値及び固有ベクトルを計算し、計算された相関行列Ｒｙとその固有値及び固有ベクトルとに基づいて、素子間結合を打ち消した素子空間に変換された相関行列Ｒｓを計算し、計算された相関行列Ｒｓの固有値及び固有ベクトルを固有値分解法を用いて計算し、計算された相関行列Ｒｓの固有値及び固有ベクトルに基づいてＭＯＤＥ法を用いて１次元又は２次元の電波到来方向に対応するパラメータを根とする多項式を解くことにより電波到来方向の角度を推定して計算する。【選択図】図１
2005062040.tif 000002
【特許請求の範囲】
【請求項１】
無線信号を受信するための励振素子と、少なくとも３本の非励振素子とを格子形状で配置し、上記各非励振素子にそれぞれ接続された可変リアクタンス素子を備え、上記各可変リアクタンス素子のリアクタンス値を変化させることにより、上記各非励振素子をそれぞれ導波器又は反射器として動作させて指向特性を変化させるアレーアンテナを用いた電波到来方向探知方法であって、
上記アレーアンテナの互いに異なる複数の放射パターンの各可変リアクタンス素子のリアクタンス値の複数のセットとによりそれぞれ設定された互いに異なる複数の放射パターンの状態において上記アレーアンテナで受信し、上記受信された各受信信号を検出し、上記検出された複数の受信信号間の相関を示す相関行列Ｒｙを計算し、上記計算された相関行列Ｒｙに基づいて固有値分解法を用いて各方位角毎の受信信号に対する固有値及び固有ベクトルを計算し、上記各方位角の方向毎の最大値の固有値に対応する固有ベクトルを用いて、上記アレーアンテナの各可変リアクタンス素子による放射パターンの各方向ベクトルに対する重み付けを表す等価的なベクトルからなる等価ウエイト行列を校正するステップと、
上記互いに異なる可変リアクタンス素子のリアクタンス値の複数のセットをそれぞれ設定したときに上記アレーアンテナによって受信される各受信信号を検出し、上記複数の受信信号間の相関を表す相関行列Ｒｙを計算し、上記計算された各相関行列Ｒｙを固有値分解して受信信号に対する上記各相関行列Ｒｙの固有値及び固有ベクトルを計算し、上記計算された相関行列Ｒｙとその固有値及び固有ベクトルとに基づいて、素子間結合を打ち消した所定の素子空間に変換された相関行列Ｒｓを計算し、計算された相関行列Ｒｓの固有値及び固有ベクトルを固有値分解法を用いて計算するステップと、
上記計算された相関行列Ｒｓの固有値及び固有ベクトルに基づいてＭＯＤＥ法を用いて、１次元又は２次元の電波到来方向に対応するパラメータを根とする多項式を解くことにより電波到来方向の角度を推定して計算するステップとを含むことを特徴とする電波到来方向探知方法。
【請求項２】
上記電波到来方向の角度を推定して計算するステップは、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する１元多項式を解くことにより電波到来方向の角度を推定して計算することを特徴とする請求項１記載の電波到来方向探知方法。
【請求項３】
上記電波到来方向の角度を推定して計算するステップは、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する２元多項式における各係数ベクトルの関数であるパラメータ行列Ｂの偏りをなくすための正規化行列Ｗｎをコレスキ分解してなる行列Ｗｎ１／２を定義し、Ｗｎ＝（ＢＨＢ）−１としたときに、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨを含む尤度関数が実質的に最小となるように、上記２つの多項式の各係数ベクトルを反復演算により計算して、２元多項式を解くことにより電波到来方向の２つの角度を推定して計算することを特徴とする請求項１記載の電波到来方向探知方法。
【請求項４】
無線信号を受信するための励振素子と、少なくとも３本の非励振素子とを格子形状で配置し、上記各非励振素子にそれぞれ接続された可変リアクタンス素子を備え、上記各可変リアクタンス素子のリアクタンス値を変化させることにより、上記各非励振素子をそれぞれ導波器又は反射器として動作させて指向特性を変化させるアレーアンテナを用いた電波到来方向探知装置であって、
上記アレーアンテナの互いに異なる複数の放射パターンの各可変リアクタンス素子のリアクタンス値の複数のセットとによりそれぞれ設定された互いに異なる複数の放射パターンの状態において上記アレーアンテナで受信し、上記受信された各受信信号を検出し、上記検出された複数の受信信号間の相関を示す相関行列Ｒｙを計算し、上記計算された相関行列Ｒｙに基づいて固有値分解法を用いて各方位角毎の受信信号に対する固有値及び固有ベクトルを計算し、上記各方位角の方向毎の最大値の固有値に対応する固有ベクトルを用いて、上記アレーアンテナの各可変リアクタンス素子による放射パターンの各方向ベクトルに対する重み付けを表す等価的なベクトルからなる等価ウエイト行列を校正する手段と、
上記互いに異なる可変リアクタンス素子のリアクタンス値の複数のセットをそれぞれ設定したときに上記アレーアンテナによって受信される各受信信号を検出し、上記複数の受信信号間の相関を表す相関行列Ｒｙを計算し、上記計算された各相関行列Ｒｙを固有値分解して受信信号に対する上記各相関行列Ｒｙの固有値及び固有ベクトルを計算し、上記計算された相関行列Ｒｙとその固有値及び固有ベクトルとに基づいて、素子間結合を打ち消した所定の素子空間に変換された相関行列Ｒｓを計算し、計算された相関行列Ｒｓの固有値及び固有ベクトルを固有値分解法を用いて計算する手段と、
上記計算された相関行列Ｒｓの固有値及び固有ベクトルに基づいてＭＯＤＥ法を用いて、１次元又は２次元の電波到来方向に対応するパラメータを根とする多項式を解くことにより電波到来方向の角度を推定して計算する手段とを備えたことを特徴とする電波到来方向探知装置。
【請求項５】
上記電波到来方向の角度を推定して計算する手段は、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する１元多項式を解くことにより電波到来方向の角度を推定して計算することを特徴とする請求項４記載の電波到来方向探知装置。
【請求項６】
上記電波到来方向の角度を推定して計算する手段は、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する２元多項式における各係数ベクトルの関数であるパラメータ行列Ｂの偏りをなくすための正規化行列Ｗｎをコレスキ分解してなる行列Ｗｎ１／２を定義し、Ｗｎ＝（ＢＨＢ）−１としたときに、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨを含む尤度関数が実質的に最小となるように、上記２つの多項式の各係数ベクトルを反復演算により計算して、２元多項式を解くことにより電波到来方向の２つの角度を推定して計算することを特徴とする請求項４記載の電波到来方向探知装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数のアンテナ素子を備えて指向特性を変化させることができるアレーアンテナを用いた電波到来方向探知方法に関し、特に、指向特性を適応的に変化させることができる電子制御導波器アレーアンテナ装置（Electronically Steerable Passive Array Radiator Antenna）を用いた電波到来方向探知方法及び装置に関する。
【背景技術】
【０００２】
アレーアンテナに入射する信号の到来角（Direction-of-Arrival：ＤＯＡ）の推定は、ディジタルアレー処理における興味ある、かつ重要な問題である。この問題を解決するために、いくつかの方法が提案されている。特に、非特許文献６に記載のＭＵＳＩＣ（MUltiple SIgnal Classification）アルゴリズムは、漸近的に不偏推定量に近づく結果を提供するものとして広く知られている。
【０００３】
ところで、特許文献１などにおいて提案されている電子制御導波器アレーアンテナ装置は、無線信号が給電される励振素子と、この励振素子から所定の間隔だけ離れて設けられ、無線信号が給電されない少なくとも１個の非励振素子と、この非励振素子に接続された可変リアクタンス素子とから成るアレーアンテナを備え、上記可変リアクタンス素子のリアクタンス値を変化させることにより、上記アレーアンテナの指向特性を変化させることができる。この電子制御導波器アレーアンテナ装置は、従来のフェーズドアレーアンテナに比べてハードウェア構成が簡単なため、魅力的なアンテナであり、近年、無線通信や高分解能方位推定などへの応用が精力的に研究されている。電子制御導波器アレーアンテナ装置は給電が１系統であり、各時刻に受信できるデータは１つのみである。この点が各時刻において全ての素子の受信データが観測可能な従来のフェーズドアレーと原理的に異なる点である。
【０００４】
電子制御導波器アレーアンテナ装置の原理から自明なように、異なる方向に指向性ビームを形成するリアクタンス値セットを求め、それらによりビームを切り替えることにより、簡易なビーム形成器が実現できる。ただし、電磁結合を利用してビームを形成するという電子制御導波器アレーアンテナ装置の原理上、素子間隔を広げることができずアンテナ開口が比較的小さなものとなる。従って、鋭いペンシルビームを形成し、分解能を向上することは困難であった。
【０００５】
小さなアンテナ開口において高分解能性を実現する手法としては、ＭＵＳＩＣ法（非特許文献２参照。）、ＥＳＰＲＩＴ（Estimation of Signal Parameters via Rotational Invariance Techniques）法（非特許文献３参照。）やＭＯＤＥ（method of Direction Estimation）法（非特許文献４参照。なお、ＷＳＦ（Weighted Subspace Fitting）法とも呼ばれる（非特許文献５参照。））などのスーパーレゾリューション法の適用が考えられる。プラポウ（Ｐｌａｐｏｕｓ）らは送信信号が時間依存性が無いものと仮定しＭＵＳＩＣ法が電子制御導波器アレーアンテナ装置に適用可能であることを示した（非特許文献６参照。）。時間依存性の無い送信信号という仮定は非常に限定されたものであるが、この概念は周期的に繰り返すプリアンブル信号を有する送信波形信号に拡張されている（非特許文献７参照。）。現実には、リアクタンス値を高精度に制御することは困難であるが、等価ウエイト行列の校正により、その問題点の解決が可能であり、ＭＵＳＩＣ法による高分解能方位角推定が可能であることが実験により実証されている（非特許文献８参照。）。ＭＵＳＩＣ法自体はコヒーレント波の分離検出が不可能であるが、平田等は得られた等価ウエイト行列を用いた逆変換（等価素子空間変換）を施すことにより、ＳＳＰ（Spatial Smoothing Processing）（非特許文献９参照。）の適用が可能となり、従ってコヒーレント波が分離検出できることを示した（非特許文献１０参照。）。
【０００６】
【特許文献１】特開２００１−２４４３１号公報。
【非特許文献１】T. Ohira et al., "Electrically steerable passive array radiator antennas for low-cost analog adaptive beamforming", IEEE International Conference Phased Array System Technology, pp.101-104, Dana Point, California, U.S.A., May 2000。
【非特許文献２】R. O. Schmidt, "Multiple emitter location and signal parameter estimation", IEEE Transactions on Antennas and Propagation, Vol.AP-34, No.3, pp.276-280, Mar.1986
【非特許文献３】R. Roy et al., "ESPRIT - Estimation of signal parameters via rotational invariance techniques", IEEE Transactions on Acoustic Speech and Signal Processing, Vol.ASSP-37, No.7, pp.984-995, July 1989。
【非特許文献４】P. Stoica et al., "Novel eigenanalysis method for direction estimation", IEE Proceeding, Vol.137, Pt.F, No.1, pp.19-26, February 1990。
【非特許文献５】H.Krim et al., "Two decades of array signal processing research", IEEE Signal Processing Magazine, Vol.13, No.4, pp.67-94, July 1996。
【非特許文献６】C. Plapous et al., “電子制御導波器アレーアンテナ装置によるリアクタンスドメインＭＵＳＩＣ法”, 電子情報通信学会研究技術報告，電子情報通信学会発行，Ｖｏｌ．ＲＣＳ２００２−１４７，ｐｐ．１−８，２００２年８月。
【非特許文献７】E. Taillefer et al., "Reactance-domain MUSIC for ESPAR antennas (Experiment)", Proceedings of Wireless Communication Network Conference, WCNC2003, pp.98-102, New Orleans, U.S.A., March 2003。
【非特許文献８】平田明史ほか，“等価ウエイト行列をキャリブレートした電子制御導波器アレーアンテナ装置の高分解能到来方向推定”,電子情報通信学会研究技術報告，電子情報通信学会発行，Ｖｏｌ．Ａ・Ｐ２００２−２２６，ｐｐ．１−８，２００３年３月。
【非特許文献９】T. J. Shan et al., "On spatial smoothing for direction-of-arrival estimation of coherent waves", IEEE Transactions on Acoustic Speech and Signal Processing, Vol.ASSP-33, No.4, pp.806-811, August1985。
【非特許文献１０】平田明史ほか，“`電子制御導波器アレーアンテナ装置を用いたリアクタンスドメインＭＵＳＩＣ法によるコヒーレント２波の到来方向推定実験”，電子情報通信学会研究技術報告，電子情報通信学会発行，Ｖｏｌ．Ａ・Ｐ２００３−２４，ｐｐ．５９−６４，２００３年３月。
【非特許文献１１】中澤達也ほか，“３次元ＭＯＤＥ法を用いた多重波パラメータ推定”，電子情報通信学会研究技術報告，電子情報通信学会発行，Ｖｏｌ．Ａ・Ｐ２００１−８８，ｐｐ．１０１−１０７，２００１年１０月。
【非特許文献１２】M. P. Clark et al., "Two-dimensional model analysis based on maximum likelihood", IEEE Transactions on Signal Processing, Vol.42, No.6, pp.1443-1451, June 1994。
【非特許文献１３】山田寛喜ほか，“アレーアンテナの素子間相互結合校正法とその精度に関する比較検討”，電子情報通信学会研究技術報告，電子情報通信学会発行，Ｖｏｌ．Ａ・Ｐ２００２−２１８，ｐｐ．１７９−１８６，２００３年３月。
【非特許文献１４】C. M. S. See, "Sensor array calibration in the presence of mutual coupling and unknown sensor gains and phases", Electronics Letters, Vol.30, No.3, pp.373-375, March 1994。
【発明の開示】
【発明が解決しようとする課題】
【０００７】
以上の検討は、全て中央に励振素子を有する円形７素子のモノポール電子制御導波器アレーアンテナ装置においてなされたものであり、ＤＯＡ推定も方位方向のみの１次元であった。また、従来技術では、インコヒーレントな２波を識別して到来方向を推定する場合において、精度はいまだ低いという問題点があった。
【０００８】
本発明の目的は以上の問題点を解決し、２次元で電波到来方向を推定でき、しかもコヒーレントやインコヒーレントの区別なく従来技術に比較して高精度で電波到来方向の推定を行うことができる電波到来方向探知方法及び装置を提供することにある。
【課題を解決するための手段】
【０００９】
第１の発明に係る電波到来方向探知方法は、無線信号を受信するための励振素子と、少なくとも３本の非励振素子とを格子形状で配置し、上記各非励振素子にそれぞれ接続された可変リアクタンス素子を備え、上記各可変リアクタンス素子のリアクタンス値を変化させることにより、上記各非励振素子をそれぞれ導波器又は反射器として動作させて指向特性を変化させるアレーアンテナを用いた電波到来方向探知方法であって、
上記アレーアンテナの互いに異なる複数の放射パターンの各可変リアクタンス素子のリアクタンス値の複数のセットとによりそれぞれ設定された互いに異なる複数の放射パターンの状態において上記アレーアンテナで受信し、上記受信された各受信信号を検出し、上記検出された複数の受信信号間の相関を示す相関行列Ｒｙを計算し、上記計算された相関行列Ｒｙに基づいて固有値分解法を用いて各方位角毎の受信信号に対する固有値及び固有ベクトルを計算し、上記各方位角の方向毎の最大値の固有値に対応する固有ベクトルを用いて、上記アレーアンテナの各可変リアクタンス素子による放射パターンの各方向ベクトルに対する重み付けを表す等価的なベクトルからなる等価ウエイト行列を校正するステップと、
上記互いに異なる可変リアクタンス素子のリアクタンス値の複数のセットをそれぞれ設定したときに上記アレーアンテナによって受信される各受信信号を検出し、上記複数の受信信号間の相関を表す相関行列Ｒｙを計算し、上記計算された各相関行列Ｒｙを固有値分解して受信信号に対する上記各相関行列Ｒｙの固有値及び固有ベクトルを計算し、上記計算された相関行列Ｒｙとその固有値及び固有ベクトルとに基づいて、素子間結合を打ち消した所定の素子空間に変換された相関行列Ｒｓを計算し、計算された相関行列Ｒｓの固有値及び固有ベクトルを固有値分解法を用いて計算するステップと、
上記計算された相関行列Ｒｓの固有値及び固有ベクトルに基づいてＭＯＤＥ法を用いて、１次元又は２次元の電波到来方向に対応するパラメータを根とする多項式を解くことにより電波到来方向の角度を推定して計算するステップとを含むことを特徴とする。
【００１０】
また、上記電波到来方向探知方法において、上記電波到来方向の角度を推定して計算するステップは、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する１元多項式を解くことにより電波到来方向の角度を推定して計算することを特徴とする。
【００１１】
さらに、上記電波到来方向探知方法において、上記電波到来方向の角度を推定して計算するステップは、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する２元多項式における各係数ベクトルの関数であるパラメータ行列Ｂの偏りをなくすための正規化行列Ｗｎをコレスキ分解してなる行列Ｗｎ１／２を定義し、Ｗｎ＝（ＢＨＢ）−１としたときに、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨを含む尤度関数が実質的に最小となるように、上記２つの多項式の各係数ベクトルを反復演算により計算して、２元多項式を解くことにより電波到来方向の２つの角度を推定して計算することを特徴とする。
【００１２】
第２の発明に係る電波到来方向探知装置は、無線信号を受信するための励振素子と、少なくとも３本の非励振素子とを格子形状で配置し、上記各非励振素子にそれぞれ接続された可変リアクタンス素子を備え、上記各可変リアクタンス素子のリアクタンス値を変化させることにより、上記各非励振素子をそれぞれ導波器又は反射器として動作させて指向特性を変化させるアレーアンテナを用いた電波到来方向探知装置であって、
上記アレーアンテナの互いに異なる複数の放射パターンの各可変リアクタンス素子のリアクタンス値の複数のセットとによりそれぞれ設定された互いに異なる複数の放射パターンの状態において上記アレーアンテナで受信し、上記受信された各受信信号を検出し、上記検出された複数の受信信号間の相関を示す相関行列Ｒｙを計算し、上記計算された相関行列Ｒｙに基づいて固有値分解法を用いて各方位角毎の受信信号に対する固有値及び固有ベクトルを計算し、上記各方位角の方向毎の最大値の固有値に対応する固有ベクトルを用いて、上記アレーアンテナの各可変リアクタンス素子による放射パターンの各方向ベクトルに対する重み付けを表す等価的なベクトルからなる等価ウエイト行列を校正する手段と、
上記互いに異なる可変リアクタンス素子のリアクタンス値の複数のセットをそれぞれ設定したときに上記アレーアンテナによって受信される各受信信号を検出し、上記複数の受信信号間の相関を表す相関行列Ｒｙを計算し、上記計算された各相関行列Ｒｙを固有値分解して受信信号に対する上記各相関行列Ｒｙの固有値及び固有ベクトルを計算し、上記計算された相関行列Ｒｙとその固有値及び固有ベクトルとに基づいて、素子間結合を打ち消した所定の素子空間に変換された相関行列Ｒｓを計算し、計算された相関行列Ｒｓの固有値及び固有ベクトルを固有値分解法を用いて計算する手段と、
上記計算された相関行列Ｒｓの固有値及び固有ベクトルに基づいてＭＯＤＥ法を用いて、１次元又は２次元の電波到来方向に対応するパラメータを根とする多項式を解くことにより電波到来方向の角度を推定して計算する手段とを備えたことを特徴とする。
【００１３】
また、上記電波到来方向探知装置において、上記電波到来方向の角度を推定して計算する手段は、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する１元多項式を解くことにより電波到来方向の角度を推定して計算することを特徴とする。
【００１４】
さらに、上記電波到来方向探知装置において、上記電波到来方向の角度を推定して計算する手段は、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨにおける行列Ｂと、到来方向との関係を有する２元多項式における各係数ベクトルの関数であるパラメータ行列Ｂの偏りをなくすための正規化行列Ｗｎをコレスキ分解してなる行列Ｗｎ１／２を定義し、Ｗｎ＝（ＢＨＢ）−１としたときに、信号の固有値からなる信号部分空間Ｅｓと直交する空間である雑音部分空間を張る行列Ｂ（ＢＨＢ）−１ＢＨを含む尤度関数が実質的に最小となるように、上記２つの多項式の各係数ベクトルを反復演算により計算して、２元多項式を解くことにより電波到来方向の２つの角度を推定して計算することを特徴とする。
【発明の効果】
【００１５】
従って、本発明に係る電波到来方向探知方法及び装置によれば、２次元で電波到来方向を推定でき、しかもコヒーレントやインコヒーレントの区別なく従来技術に比較して高精度で電波到来方向の推定を行うことができる。
【発明を実施するための最良の形態】
【００１６】
以下、本発明に係る実施形態について図面を参照して説明する。なお、同様の構成要素については同一の符号を付している。
【００１７】
本発明者らは、特許文献１などに開示された円形の電子制御導波器アレーアンテナ装置は平面アレーアンテナであるので、原理的には２次元推定も可能であるということに着目して、コヒーレント波の２次元到来方向推定を目的とした場合、方形配列の電子制御導波器アレーアンテナ装置を採用すると、各種アルゴリズムの適用が可能となることを発見した。本実施形態では、コヒーレント波を直接分離できる最尤推定手法であるＭＯＤＥ法を採用し、一例として中央のアンテナ素子を励振素子とした９素子（３×３）モノポール方形の電子制御導波器アレーアンテナ装置を用いた２次元のＤＯＡ推定法について説明する。アンテナ素子をモノポールとしているため、本質的に天頂（仰角＝９０度）方向の推定は不向きであるが、ＳＳＰを用いることなくマルチパス波の分離推定が可能となる。ここでは、モーメント法による解析を行い、電子制御導波器アレーアンテナ装置を用いた２次元ＭＯＤＥ法の基本特性を明らかにする。
【００１８】
図１は、９素子の方形電子制御導波器アレーアンテナ装置であるアレーアンテナ装置１００を用いた、本発明の実施形態に係る電波到来方向探知装置の構成を示すブロック図であり、図２は図１のアレーアンテナ装置１００の詳細な構成を示す断面図である。
【００１９】
図１において、アレーアンテナ装置１００は、接地導体１１上に３行３列で直交するｘｙ座標系に対して平行な位置に格子状に設けられた９本のアンテナ素子から構成され、ここで、９本のアンテナ素子は好ましくは、ｘ軸方向の間隔Δｘでかつｙ軸方向の間隔Δｙで互いに離間してｘ軸及びｙ軸に対して設けられている。これら９本のアンテナ素子のうち、ｘｙ座標系の原点Ｏに配置された中央のアンテナ素子を励振素子Ａ５とし、その周囲のアンテナ素子を非励振素子Ａ１乃至Ａ９（Ａ５を除く。）として用いる。励振素子Ａ５及び各非励振素子Ａ１乃至Ａ９（Ａ５を除く。）の長さは、例えば約λ／４（但し、λは所望波の波長である。）になるように構成される。また、接地導体１１は、４Δｘ以上のｘ軸方向の長さと４Δｙ以上のｙ軸方向の長さを有し、実質的に電気的に無限大の大きさを有する矩形形状を有する。励振素子Ａ５の給電点は同軸ケーブル９を介して無線受信機４の低雑音増幅器（ＬＮＡ）１に接続され、また、非励振素子Ａ１乃至Ａ９（Ａ５を除く。）はそれぞれ可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）に接続され、これら可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）のリアクタンス値はリアクタンス値コントローラ１０からのリアクタンス値信号によって設定される。
【００２０】
図２のアレーアンテナ装置１００の縦断面図において、励振素子Ａ５は接地導体１１と電気的に絶縁され、各非励振素子Ａ１乃至Ａ９（Ａ５を除く。）はそれぞれ、可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）を介して、接地導体１１に対して高周波的に接地される。可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）は、例えば、制御電圧（又はバイアス電圧）が印加されることによってそのリアクタンス値が変化する可変容量ダイオードであって、制御電圧はリアクタンス値コントローラ１０からのリアクタンス値信号を介して印加される。リアクタンス値コントローラ１０は、デジタルシグナルプロセッサをベースとするコントローラであって、リアクタンス値テーブルメモリ１３内に予め設定されたディジタル電圧値を参照し、内蔵した８個のＤ／Ａ変換器（図示せず。）を使って上記ディジタル電圧値をアナログの制御電圧値に変換し、この制御電圧値をリアクタンス値信号としてアレーアンテナ装置１００の各可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）に印加することにより、アレーアンテナ装置１００上で、対応する各指向性ビームパターンが形成される。
【００２１】
可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）の動作を説明すると、例えば励振素子Ａ５と非励振素子Ａ１乃至Ａ９（Ａ５を除く。）の長手方向の長さが実質的に同一であるとき、例えば、可変リアクタンス素子１２−１がインダクタンス性（Ｌ性）を有するときは、可変リアクタンス素子１２−１は延長コイルとなり、非励振素子Ａ１の電気長が励振素子Ａ５に比較して長くなり、反射器として働く。一方、例えば、可変リアクタンス素子１２−１がキャパシタンス性（Ｃ性）を有するときは、可変リアクタンス素子１２−１は短縮コンデンサとなり、非励振素子Ａ１の電気長が励振素子Ａ５に比較して短くなり、導波器として働く。また、他の可変リアクタンス素子１２−２乃至１２−９（１２−５を除く。）に接続された非励振素子Ａ２乃至Ａ９（Ａ５を除く。）についても同様に動作する。従って、図１のアレーアンテナ装置１００において、各非励振素子Ａ１乃至Ａ９（Ａ５を除く。）に接続された可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）のリアクタンス値を変化させることにより、アレーアンテナ装置１００の平面指向性特性を変化させることができる。
【００２２】
図１のアレーアンテナの制御装置において、アレーアンテナ装置１００は、図１の左上に図示するように、ｘｙ平面上でｘ軸からの右回りの角度で定義される方位角φと、ｚ軸方向である天頂方向からの仰角（以下、天頂仰角という。）θの到来方向で無線信号を受信し、上記受信された無線信号である受信信号は、励振素子Ａ５に接続された同軸ケーブル９から出力される。出力された受信信号は、無線受信機４の低雑音増幅器１を介してダウンコンバータ２に入力され、ダウンコンバータ２は入力される受信信号を所定の中間周波数の中間周波信号に周波数変換した後、Ａ／Ｄ変換器３に出力する。Ａ／Ｄ変換器３は、入力されるアナログの中間周波信号をディジタルの中間周波信号に変換した後、電波到来方向探知コンピュータ２０に出力する。さらに、電波到来方向探知コンピュータ２０は、入力される中間周波信号に基づいて、その中間周波信号の相関行列Ｒｙ（後述する数４参照。）を計算し、本発明に係る２次元ＭＯＤＥ法を用いて受信された無線信号の到来角を計算し、その結果をＣＲＴディスプレイ２１に出力して表示する。ここで、電波到来方向探知コンピュータ２０は、本実施形態では、図３に示すように、等価ウエイト行列を校正して計算する処理（ステップＳ１）を実行した後、固有値分解処理（ステップＳ２）を実行し、さらに、到来角推定処理（ステップＳ３）を実行することにより天頂仰角θと方位角φの到来方向を推定して計算してＣＲＴディスプレイ２１に表示する。
【００２３】
なお、上記ステップＳ１の処理を実行するときに、電波到来方向探知コンピュータ２０からの制御信号に基づいて、アレーアンテナ装置１００を励振素子Ａ５を中心軸として機械的に回転させる回転機構３０が設けられている。
【００２４】
次いで、本実施形態における問題の定式化について以下に説明する。図１のアレーアンテナ装置１００において、ｄ個の入射波がそれぞれ（θｋ，φｋ）方向（ｋ＝１，２，…，ｄ）から到来している場合を考える。各入射波の複素振幅はｓｋ（ｔ）とする。もし各アンテナ素子Ａ１乃至Ａ９での受信信号が観測可能であれば、ｘｙ座標系の座標（ｘｉ，ｙｊ）に位置する素子の受信信号ｒ（ｘｉ，ｙｊ，ｔ）は次式のように表される。なお、当該明細書において、数式がイメージ入力された墨付き括弧の数番号と、数式が文字入力された大括弧の数式番号とを混在して用いており、また、当該明細書での一連の数式番号として「式（１）」の形式を用いて数式番号を式の最後部に付与して（なお、簡単化のため付与していない数式も存在する。）用いることとする。
【００２５】
【数１】
2005062040.tif 000003
【００２６】
ここで、ｆｃは無線信号の搬送波周波数であり、ｃは光速である。本実施形態に係るアレーアンテナ装置１００では、図１に示すように、中央に位置する１つのアンテナ素子である励振素子Ａ５から出力される受信信号ｒ（ｘ２，ｙ２）を観測し、他の非励振素子Ａ１乃至Ａ９（Ａ５を除く。）には可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）が接続されている。このとき、アレーアンテナ装置１００による受信信号ｙｍ（ｔ）は、次式で表される。
【００２７】
［数式１］
ｙｍ（ｔ）＝ｗｍＴｒ（ｔ）＋ｎｍ（ｔ）      （２）
【００２８】
ただし、ｒ（ｔ）は次式で表されるベクトルである。
【００２９】
［数式２］
ｒ（ｔ）
＝［ｒ（ｘ１，ｙ１，ｔ），ｒ（ｘ１，ｙ２，ｔ），
ｒ（ｘ１，ｙ３，ｔ），…，
ｒ（ｘ３，ｙ１，ｔ），ｒ（ｘ３，ｙ２，ｔ），
ｒ（ｘ３，ｙ３，ｔ）］Ｔ
（３）
【００３０】
ここで、Ｔは行列の転置を表す記号である。ｗｍは非励振素子のリアクタンス値セットにより決定される等価ウエイトベクトルであり、次式で表される。
【００３１】
［数式３］
ｗｍ
＝２Ｒｅ［Ｚ０］
（Ｚ＋ｄｉａｇ［ｊｘｍ１，…，ｊｘｍ４，Ｚ０，ｊｘｍ６，…，ｊｘｍ９］）−１ｕ０
（４）
【００３２】
ここで、Ｒｅ［Ｚ０］は無線受信機４の内部インピーダンスＺ０の実数部であり、Ｚはアレーアンテナ装置１００の素子間相互結合に関するインピーダンス行列である。また、ｕ０は次式で定義される単位ベクトルである。さらに、ｄｉａｇ［・］は対角行列を意味する。また、ｎｍ（ｔ）は受信機で生じる白色ガウス雑音である。
【００３３】
［数式４］
ｕ０＝［０，…，０，１，０，…，０］Ｔ
【００３４】
次いで、本実施形態で用いる等価素子空間変換について以下詳細説明する。本実施形態で適用する２次元ＭＯＤＥ法もスーパーレゾリューション法の一つであり、本来、各素子の受信信号の振幅、位相情報を用いた高分解能到来方向推定手法である。従って、各素子の受信データが得られないアレーアンテナ装置１００では、直接適用することはできない。本発明者らは、各リアクタンス値セット毎に同一波形を繰り返し送信するようなプリアンブル信号を、送信信号に付加し、式（１）に示すような素子間結合を打ち消した素子空間データに変換する。なお、等価素子空間変換については、非特許文献１０において開示している。リアクタンス値セット切り替え間隔（フレーム間隔）をＴとし、サンプリング間隔をΔｔとすると、送信プリアンブル信号は次式のように表される。
【００３５】
［数式５］
ｓｋ（ｋＴ＋ｉΔｔ）＝ｓｋ（ｌＴ＋ｉΔｔ），
ｋ＝１，２，…，Ｍ，
ｌ＝１，２，…，Ｍ，
ｉ＝１，２，…，ｎｍａｘ
（５）
【００３６】
ここで、
［数式６］
ｎｍａｘΔｔ≦Ｔ
である。
【００３７】
Ｍは用いるリアクタンス値セットの数である。このような送信信号を用いてＭ個の独立なリアクタンス値セットにより得られた受信信号をベクトル表示すると、式（１）及び式（２）より、次式が得られる。
【００３８】
［数式７］
ｙ（ｉΔｔ）
＝［ｙ１（ｉΔｔ），ｙ２（ｉΔｔ），…，ｙＭ（ｉΔｔ）］Ｔ
＝ＷＡｓ（ｉΔｔ）＋ｎ      （６）
【００３９】
ここで、等価ウエイト行列Ｗと、ステアリング行列Ａは次式で表される。
【００４０】
【数２】
2005062040.tif 000004
【００４１】
［数式８］
Ａ＝［ａ（θ１，φ１），ａ（θ２，φ２），…，ａ（θｄ，φｄ）］
（７ｂ）
【００４２】
［数式９］
ａ（θｋ，φｋ）
＝［ｅｘｐ（ｊ２πｆｃ（（ｘ１／ｃ）ｓｉｎθｋｃｏｓφｋ
＋（ｙ１／ｃ）ｓｉｎθｋｓｉｎφｋ）），
ｅｘｐ（ｊ２πｆｃ（（ｘ１／ｃ）ｓｉｎθｋｃｏｓφｋ
＋（ｙ２／ｃ）ｓｉｎθｋｓｉｎφｋ）），…，…，
ｅｘｐ（ｊ２πｆｃ（（ｘ３／ｃ）ｓｉｎθｋｃｏｓφｋ
＋（ｙ３／ｃ）ｓｉｎθｋｓｉｎφｋ））］Ｔ
（７ｃ）
【００４３】
［数式１０］
ｓ（ｉΔｔ）
＝［ｓ１（ｉΔｔ），…，ｓｄ（ｉΔｔ）］Ｔ
（７ｄ）
【００４４】
式（６）より、電子制御導波器アレーアンテナ装置においても等価ウエイト行列Ｗが正確に既知であれば、以下の変換により、通常のスーパーレゾリューションアレーアンテナで用いられる等価な素子空間における信号相関行列ＲＳが得られる。
【００４５】
［数式１１］
ＲＳ
＝Ｗ−１（Ｒｙ−σ２Ｉ）（Ｗ−１）Ｈ
＝ＡＳＡＨ      （８ａ）
【００４６】
［数式１２］
Ｒｙ＝Ｅ［ｙｙＨ］      （８ｂ）
【００４７】
［数式１３］
Ｓ＝Ｅ［ｓｓＨ］      （８ｃ）
【００４８】
上記式（８ａ），式（８ｂ），式（８ｃ）が等価素子空間変換を表す式である。ここで、Ｅ［・］はエルゴード性を仮定した所定の時間期間における時間平均値（又はアンサンブル平均値）であり、上付き添字Ｈは行列のエルミート転置を表す。また、σ２は雑音電力である。電波環境の変動が無視できる時間内に、上述のデータが得られれば、この変換により、スーパーレゾリューションアレーアンテナにおいてアダプティブで用いられる様々な信号処理が適用可能となる。
【００４９】
等価ウエイト行列Ｗに対応するリアクタンス値セットのリアクタンス値は、可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）への印加電圧により決定される。また、アンテナ素子間の素子間相互結合は、公知のモーメント法などにより導出可能である。ただし、可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）等の素子のばらつきが存在すると設定値と異なることとなる。スーパーレゾリューション法は、誤差要因による特性劣化が大きいため、そのような場合には、等価ウエイト行列自体を推定（校正）する必要があり（非特許文献１０参照。）、これについては詳細後述する。
【００５０】
電子制御導波器アレーアンテナ装置であるアレーアンテナ装置１００は、上述のごとく、非励振素子のリアクタンスを制御することで、アンテナの指向性を変化させるアレーアンテナである。すなわち、リアクタンス値を変化させて得られたリアクタンスドメインの受信データは、アンテナ指向性の観点からは、成型ビームにより得られた受信データ、すなわちビームスペース領域の受信データに対応する。すなわち、式（６）は通常のスーパーレゾリューションアレーのビームスペース変換式と同形である。
【００５１】
次いで、受信信号ｙ（ｔ）と相関行列Ｒｙとの関係について以下に説明する。上記の等価ウエイト行列Ｗは、アレーアンテナ装置１００において可変リアクタンス素子による放射パターンの各方向ベクトルに対する重み付けを表す等価的な行列であり、これを用いてベクトル表記の受信信号系列は次式で表される。
【００５２】
【数３】
2005062040.tif 000005
【００５３】
このとき、相関行列と固有ベクトルの関係から、複数の受信信号から検出した複数のシーケンス信号間の相関行列Ｒｙは次式で表される。
【００５４】
【数４】
2005062040.tif 000006
【００５５】
ここで、λｉは固有値であり、ｅｌは信号部分空間を張る第ｌの固有ベクトル（ｌ＝１，２，…，Ｌ）である。すなわち、相関行列Ｒｙを固有値分解することにより９個の固有値と９個の固有ベクトルを得ることができる。なお、例えば、アレーアンテナ装置１００のＭＵＳＩＣスペクトラムは等価ウエイト行列Ｗを用いて、次式で表される（非特許文献１０参照。）。
【００５６】
【数５】
2005062040.tif 000007
【００５７】
ここで、
【数６】
2005062040.tif 000008
である。
【００５８】
すなわち、ＭＵＳＩＣ法では、雑音の固有ベクトルからなる行列ＥＮによりＭＵＳＩＣスペクトルを計算でき、これに基づいて到来角を計算できるが、本実施形態では、変形された２次元ＭＯＤＥ法を用いる。
【００５９】
次いで、オムニパターンを含めた９つの指向性パターンを用いて、信号部分空間からアレーアンテナ装置１００の「等価ウエイト行列」の校正処理（図３のステップＳ１）について以下に説明する。
【００６０】
上記インピーダンス行列Ｚの各要素は未知である。鋭いＭＵＳＩＣスペクトラムを求め角度分解能を向上させるためには、このインピーダンス行列Ｚを含んだ等価ウエイト行列Ｗを校正する必要がある。本実施形態では、校正方法として９つの指向性パターンで受信した受信信号系列を用いて、信号の固有値に対応する信号固有ベクトルによって等価ウエイト行列を推定する方法を提案する。校正では到来波は１波のみ（Ｌ＝１）とし、平面波を２４方位から順次に到来させる。
【００６１】
校正の手順の具体的な詳細は後述するが、その概略を説明すると、図４に示すように、まず、受信アンテナであるアレーアンテナ装置１００から見た電波到来方向を０°に固定する（ステップＳ１２でθ＝０゜）。次に、可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）のリアクタンス値を、表１のように９通り（ｍ＝０，１，…，８）に設定し、指向性パターンを回転させて受信信号系列ｙｍ（ｔ）を取得する（ステップＳ１３−Ｓ１６）。ｔはサンプリング時刻を表し、１シンボルにつき１回サンプリングを行う。指向性パターン毎にＰサンプル（＝Ｐシンボル）の信号を受信する。受信信号系列を上記数３のようにベクトル表記し（ステップＳ１８）、上記数４の相関行列Ｒｙを固有値分解した後（ステップＳ１９）、信号の固有値（最大値の固有値）に対応する信号固有ベクトルｅ１（１）を取得する（ステップＳ２０）。
【００６２】
本実施形態で用いるリアクタンス値セットを表１に示す。表１において、−９０Ωは可変リアクタンス素子１２のリアクタンス値最小値であり、−４．７７Ωは可変リアクタンス素子１２のリアクタンス値最大値である。ｍ＝０の放射パターンは正方形配列であるために方位角９０度毎に周期性を有するパターンであり、ｍ＝１乃至８は、それぞれ４５度間隔の方位角での主ビームを有するセクタパターンである。
【００６３】
【表１】
2005062040.tif 000009
【００６４】
次いで、電波到来方向を１５°に固定し（ステップＳ１２）、同様の手順によって信号固有ベクトルｅ１（２）を取得する（ステップＳ１３−Ｓ２０）。以降、電波到来方向を０°から１５°間隔で２４方位に順に設定し、同様の手順を繰り返すと２４個の信号固有ベクトルが得られる（Ｓ１２−Ｓ２２の繰り返しフロー）。電波到来方向は既知であり、その方向の方向ベクトルも既知であるため、次式が成立する。
【００６５】
【数７】
2005062040.tif 000010
【００６６】
ここで、上述の校正手順で計算された各方位角での最大値の固有値に対応する固有ベクトル［ｅ１（１）  ｅ１（２）  …  ｅ１（２４）］を行列Ｅｓと置き換え、各方位角毎のステアリングベクトルからなるステアリング行列［ａ（φ１）  ａ（φ２）  …  ａ（φ２４）］を行列Ａと置き換えると、等価ウエイト行列Ｗは次式で表される（非特許文献１４参照。）。
【００６７】
【数８】
2005062040.tif 000011
【００６８】
上記数８から明らかなように、等価ウエイト行列Ｗを、上述の校正手順で計算された各方位角での最大値の固有値からなる固有ベクトルの行列Ｅｓと、各方位角毎のステアリングベクトルからなるステアリング行列Ａとを用いて校正している。
【００６９】
次いで、アレーセットを用いた素子空間への変換について以下に説明する。上記数３における等価ウエイト行列Ｗは、信号部分空間を用いた校正処理によって精度よく求めることができる。以降のＷは校正処理後の等価ウエイト行列を表す。等価ウエイト行列Ｗはアンテナ素子間の素子間結合を含んでおり、相関行列ＲｙはそのままではＭＯＤＥ法を適用することができない。そこで、等価ウエイト行列Ｗの逆行列によって素子間結合を打ち消した素子空間（エレメントスペース）への変換を次式を用いて行う。
【００７０】
［数式１４］
Ｒｓ＝Ｗ−１（Ｒｙ−σ２）（Ｗ−１）Ｈ
【００７１】
ここで、相関行列Ｒｙを固有値分解することにより固有値及び固有ベクトルを計算することができ、こうして得られるＬ個の信号の固有値は、λｉ（ｉ＝１，２，…，Ｌ）で表すことができ、それ以外の雑音の固有値はλｉ（ｉ＝Ｌ＋１，Ｌ＋２，…，９）で表すことができるので、当該（９−Ｌ）個の雑音の固有値を平均した雑音電力の推定値の２乗値σ２は次式で求められる。
【００７２】
【数９】
2005062040.tif 000012
【００７３】
次いで、本実施形態に係る２次元ＭＯＤＥ法について以下に説明する。等価素子空間変換により、様々なスーパーレゾリューション法の適用が可能となる。特に既に、図１で示した方形アレーではＳＡＧＥ（Space Alternating Generalized Expectation Maximization Algorithm）法やＥＭ（Expectation Maximization Algorithm）法などのアレー形状に依存しないものは言うまでも無く、ＭＵＳＩＣ法、ＥＳＰＲＩＴ法、ＭＯＤＥ法など、非常に多くのアルゴリズムが適用可能となる。
【００７４】
ＭＯＤＥ法は、ＭＵＳＩＣ法やＥＳＰＲＩＴ法とは異なり、ＳＳＰ等の相関抑圧手法を適用することなくコヒーレント波の分離検出が可能な手法である。ＳＳＰを適用するとサブアレー化により実効アレー開口が減少し、それが推定精度の劣化をもたらす。サブアレー化を必要としないＭＯＤＥ法では、そのような劣化が存在しない（非特許文献１１参照。）。ＭＯＤＥ法の２次元への拡張は、非特許文献１２により報告されているＩＱＭＬ（Interactive Quadratic Maximum Likelihood）法の定式化を利用している。以下に、２次元ＭＯＤＥ法の概要を示す。式（８ａ）で得られた等価素子空間相関行列を固有値展開すると、次式を得る。
【００７５】
［数式１５］
ＲＳ＝ＥＳＡＳＥＳＨ      （９ａ）
［数式１６］
ＥＳ＝［ｅ１，…，ｅｄ’］      （９ｂ）
［数式１７］
ΛＳ＝ｄｉａｇ［λ１，…，λｄ’］      （９ｃ）
【００７６】
ここで、λｄ’は行列ＲＳの主要な固有値であり、全ての入射波がインコヒーレントならばｄ’＝ｄであり、コヒーレントな波が存在する場合にはｄ’＜ｄである。なお、行列ＲＳでは既に雑音成分が除去されていることに注意を要する。ＭＯＤＥ法では得られた信号部分空間ＥＳの直交空間、すなわち雑音部分空間を推定することにより、入射波の到来方向を決定する。この尤度関数は次式で与えられ、本実施形態では、この尤度関数が最小となるような行列Ｂを計算する。  
【００７７】
［数式１８］
Ｆｗ
＝Ｔｒ［Ｂ（ＢＨＢ）−１ＢＨＲＳ］
＝Ｔｒ［Ｂ（ＢＨＢ）−１ＢＨＥＳΛＳＥＳＨ］      （１０）
【００７８】
ここで、Ｔｒ［・］は引数の行列のトレース又は対角和である。また、
［数式１９］
Ｂ（ＢＨＢ）−１ＢＨ
が雑音部分空間を張る行列である。図１のｘ軸上に一直線上でアンテナ素子を所定の間隔で配置してなるリニアレーアンテナでの１次元ＤＯＡ推定では、この行列Ｂと方位角φｋの間に次式の関係が成立する。
【００７９】
［性質］行列Ｂと方位角φｋの関係（リニアアレー１次元推定）：
【００８０】
［数式２０］
ｚｋ＝ｅｘｐ（ｊ２π（ｆΔｘ／ｃ）ｃｏｓφｋ），
（ｋ＝１，２，…，ｄ）
【００８１】
上記式を根とする次式の多項式を考える。
【００８２】
【数１０】
2005062040.tif 000013
【００８３】
このときの多項式の係数より、次式のように、係数ベクトルｂ及び係数行列Ｂを定義する。
【００８４】
［数式２１］
ｂ＝［ｂ０，ｂ１，…，ｂｄ］Ｔ      （１２ａ）
【００８５】
【数１１】
2005062040.tif 000014
【００８６】
このとき、行列ＢとモードベクトルＡには次式の関係が成立する。
【００８７】
［数式２２］
ＢＨＡ＝０      （１３）
【００８８】
この関係より、雑音部分空間の形式として式（１２ａ）の形を採用し、係数ベクトルｂを推定することにより到来方向推定が可能となることが分かる。この性質を２次元に拡張するためには入射波の到来天頂仰角θｋ及び到来方位角φｋに関する２つの多項式が必要となる。以下の関係を有する多項式が存在することが知られている。
【００８９】
［性質２］多項式の係数と（θｋ、φｋ）の関係（方形アレー２次元推定）：
次の２つの多項式を定義する。
【００９０】
［数式２３］
Ｕ（ｚ）＝ｕ０ｚｄ＋ｕ１ｚｄ−１＋…＋ｕｄ−１ｚ＋ｕｄ      （１４ａ）
［数式２４］
Ｖ（ｚ）＝ｖ０ｚｄ−１＋ｖ１ｚｄ−２＋…＋ｖｄ−２ｚ＋ｖｄ−１      （１４ｂ）
【００９１】
ここで、多項式Ｕ（ｚ）と多項式Ｖ（ｚ）では項の数が１つだけ違っている理由はＭＯＤＥ法がその性質を利用するものであるからであり、多項式Ｕ（ｚ）と多項式Ｖ（ｚ）はベクトルや行列ではなくスカラー量である。ここで、次式のように定義する。
【００９２】
［数式２５］
ｕ＝［ｕ０，…，ｕｄ］Ｔ
［数式２６］
ｖ＝［ｖ０，…，ｖｄ−１］Ｔ
［数式２７］
μｋ（１）＝ｓｉｎθｋｃｏｓφｋ
［数式２８］
μｋ（２）＝ｓｉｎθｋｓｉｎφｋ
【００９３】
上記のごとく定義すると、次式のように、係数ベクトルセット（ｕ，ｖ）と方位角係数セット（μｋ（２），μｋ（１））を対応付ける写像ψが存在する。
【００９４】
［数式２９］
（ｕ，ｖ）
＝ψ（（ｅｘｐ（ｊ２π（ｆΔｙ／ｃ）μｋ（２）），
ｅｘｐ（ｊ２π（ｆΔｘ／ｃ）μｋ（１））））
（１５）
【００９５】
このとき、式（１４ａ）及び式（１４ｂ）は次式の関係を満たす。
【００９６】
［数式３０］
Ｕ（ｅｘｐ（ｊ２π（ｆΔｙ／ｃ）μｋ（２）））＝０，
ｋ＝１，２，…，ｄ      （１６ａ）
【００９７】
［数式３１］
Ｖ（ｅｘｐ（ｊ２π（ｆΔｙ／ｃ）μｋ（２）））
＝ｅｘｐ（ｊ２π（ｆΔｘ／ｃ）μｋ（１）），
ｋ＝１，２，…，ｄ      （１６ｂ）
【００９８】
ここで、行列Ｕは（Ｌｙ−ｄ）×Ｌｙの行列である。この性質より２次元ＭＯＤＥ法ではθ，φのパラメータの組推定が可能となる。このような多項式の係数ベクトルｕ，ｖを求めるのが次の問題である。まず、式（１０）を次式のように変形することができる。
【００９９】
［数式３２］
Ｆｗ（ｕ，ｖ）＝‖Ｗｎ１／２ＢＨＥＳΛ１／２‖２      （１７ａ）
【０１００】
［数式３３］
Ｗｎ＝（ＢＨＢ）−１      （１７ｂ）
【０１０１】
ここで、行列Ｗｎ１／２は正規化行列Ｗのコレスキ分解後の行列であり、行列Ｂの偏りをなくすための正規化行列である。行列Ｂが係数ベクトルｕ，ｖの関数であるパラメータ行列である。また、行列Λ１／２は式（９ｃ）を用いて対角行列ｄｉａｇ［√（λ１），√（λ２），…，√（λｄ）］（ｄ’＝ｄの場合）で表される。いま、クロネッカー（ｋｒｏｎｅｃｋｅｒ）演算子を導入すると、モードベクトルは次式で表される。
【０１０２】
［数式３４］
ａ（θｋ，φｋ）＝ａ（μｋ（１），μｋ（２））
【０１０３】
ここで、係数ベクトルｕに着目して上記［性質１］に基づき行列Ｕを構築すると、次式が成立する。
【０１０４】
［数式３５］
Ｕａ（μｋ（２））＝０      （１８）
【０１０５】
従って、次式の関係が求められる。
【０１０６】
【数１２】
2005062040.tif 000015
【０１０７】
ここで、Ｉｉはｉ×ｉの単位行列を表すものとする。同様に係数ベクトルｖに着目して［性質１］と同様に行列Ｖを定義すると［性質２］より次式を得る。
【０１０８】
［数式３６］
Ｖａ（μｋ（２））
＝ｅｘｐ（ｊ２π（ｆΔｘ／ｃ）μｋ（１））Ｉｈａ（μｋ（２））
（２０ａ）
【０１０９】
［数式３７］
Ｉｈ＝［ＩＬｙ−ｄ＋１｜０，…０］
（２０ｂ）
【０１１０】
ここで、Ｉｈは（Ｌｙ−ｄ＋１）×Ｌｙ行列となる。上式の両辺に
［数式３８］
ｅｘｐ（ｊ２π（ｉｆΔｘ／ｃ）μｋ（１））
ｉ＝１，２，…，Ｌｘ−２
を乗算して並べると次式が得られる。
【０１１１】
【数１３】
2005062040.tif 000016
【０１１２】
ここで、ａ（μｋ（２））はＬｙ×１のベクトルであり、ａ（μｋ（１））はＬｘ×１のベクトルである。従って、
【数１４】
2005062040.tif 000017
は（Ｌｘ×Ｌｙ）×１のベクトルである。また、Ｖは（Ｌｙ−ｄ＋１）×Ｌｙの行列である。式（１９）及び式（２１）の各右辺は信号部分空間に対応し、この両式を積み重ねても、その直交関係は成立する。それを次式のごとく２次元ＭＯＤＥ法の行列Ｂとする。
【０１１３】
【数１５】
2005062040.tif 000018
【０１１４】
上記式の右辺の分母は、｛（Ｌｘ−１）×（Ｌｙ−ｄ＋１）｝×（Ｌｘ×Ｌｙ）の行列であり、ｄは入射波の個数ｄである。この行列Ｂに関して式（１７ａ）の尤度関数が最小となるよう係数ベクトルｕ，ｖを決定すればよい。そのような係数ベクトルｕ，ｖが得られれば、上記［性質２］から多項式の根として２次元ＤＯＡの推定が可能となる。さて、式（２２）において係数ベクトルｕ，ｖは複雑な形で組み込まれているが、若干の変形により、次式のような係数ベクトルｕ，ｖに関する連立方程式に帰着可能できる。
【０１１５】
【数１６】
2005062040.tif 000019
【０１１６】
ここで、上記式の導出は以下のようにできる。式（１７ａ）の右辺におけるＢＨＥＳΛ１／２は次式で与えられる。
【０１１７】
【数１７】
2005062040.tif 000020
【０１１８】
上記数１７を用いると次式のように変形することができ、式（２３）を導出できる。
【０１１９】
【数１８】
2005062040.tif 000021
【０１２０】
なお、行列Ｇｉの構造については当該明細書の最後部に添付の付録を参照されたい。この方程式により得られたベクトルｆ、すなわち係数ベクトルｕ，ｖによる多項式の根の関係から２次元到来方向が推定可能となる。ただし、上式には、行列Ｂによる雑音部分空間を正規化するための正規化行列Ｗｎ１／２が含まれている。実際には、正規化行列Ｗｎ１／２は未知であるので、通常、初期値を単位行列とし、得られた推定結果から正規化行列Ｗｎ１／２を構築し推定を繰り返す反復推定が必要となる。ただし、行列Ｗ１／２は行列Ｂの偏りをなくすための正規化行列であり、一般に、
［数式３９］
Ｗｎ１／２＝Ｉ
とした初回の推定結果においても良好な推定結果が得られ、数回の反復で収束する。
【０１２１】
ここでは、方位角係数μｋ（２），μｋ（１）の順に推定するよう定式化したが、データベクトルの順序を変えることにより、方位角係数μｋ（１），μｋ（２）の順に推定するよう多項式の係数を定義することも可能である。現実には、同一の天頂仰角θｋ又は方位角φｋを有する場合も有り得るため、推定順序の異なる２つの推定結果をマージする工夫が必要となる。
【０１２２】
図３は、図１の電波到来方向探知装置の電波到来方向探知コンピュータ２０によって実行されるメインルーチンである電波到来方向推定処理を示すフローチャートである。図３のステップＳ１において、図４の等価ウエイト行列を校正して計算する処理を実行し、ステップＳ２において、図５の固有値分解処理を実行し、ステップＳ３において、図６の到来角推定処理を実行して当該電波到来方向推定処理を終了する。
【０１２３】
図４は、図３のサブルーチンである等価ウエイト行列を校正して計算する処理（ステップＳ１）を示すフローチャートである。なお、リアクタンス値テーブルメモリ１３には、放射パターンを順次設定するために、例えば、表１に示した９組のリアクタンス値セット（ｍ＝０，１，２，…，８）が格納されている。ただし、これは一例であって、少なくとも１組のオムニパターンと複数組のセクタパターンのリアクタンス値セットを含む、互いに異なる複数のリアクタンス値セットであってもよい。もしくは、オムニパターンを含まない、複数組のセクタパターンのリアクタンス値セットを含む、互いに異なる複数のリアクタンス値セットであってもよい。
【０１２４】
図４のステップＳ１１において、電波到来方向設定パラメータｎを１に初期化し、ステップＳ１２において、回転機構３０を用いてアレーアンテナ装置１００を、回転機構３０を用いて機械的に回転させて送信無線信号に対して到来方位角φ＝１５×（ｎ−１）［度］となるように設定する。次いで、ステップＳ１３において、放射パターンパラメータｍを０に初期化し、ステップＳ１４において、リアクタンス値コントローラ１０を制御してリアクタンス値テーブルメモリ１３からｍ番目のリアクタンス値セットを読み出してアレーアンテナ装置１００の各可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）に設定する。そして、ステップＳ１５において、受信信号ｙｍ（ｔ）を受信して電波到来方向探知コンピュータ２０内の一時記憶メモリに格納し、ステップＳ１６において放射パターンパラメータｍ≧８であるか否かが判断され、ＮＯであるときは、ステップＳ１７で放射パターンパラメータｍを１だけインクリメントした後、ステップＳ１４に戻る。一方、ステップＳ１６でＹＥＳであるときは、ステップＳ１８に進む。
【０１２５】
次いで、ステップＳ１８で一時記憶メモリに格納された９個の受信信号ｙｍ（ｔ）を数３のように受信信号ベクトルとして電波到来方向探知コンピュータ２０内の一時記憶メモリに格納し、ステップＳ１９で、受信信号ベクトルに基づいて数４を用いて相関行列Ｒｙを計算し、固有値分解法を用いて受信信号に対する複数の固有値及び複数の固有ベクトルを計算する。さらに、ステップＳ２０において、計算された複数の固有ベクトルのうち信号の固有値（最大の固有値）に対応する固有ベクトルをｅ１（ｎ）とし、ステップＳ２１において電波到来方向設定パラメータｎ≧２４であるか否かが判断され、ＮＯのときはステップＳ２２に進み、電波到来方向設定パラメータｎを１だけインクリメントした後、ステップＳ１２に戻る。一方、ステップＳ２１でＹＥＳであるときは、ステップＳ２３に進み、得られた信号固有ベクトルからなる行列ＥＳと、各放射パターンのステアリングベクトルからなるステアリング行列Ａとに基づいて、数７に基づく数８を用いて校正された等価ウエイト行列Ｗを推定して計算して、元のメインルーチンに戻る。
【０１２６】
図５は、図３のサブルーチンである固有値分解処理（ステップＳ２）を示すフローチャートである。
【０１２７】
図５のステップＳ３１において放射パターンパラメータｍを０に初期化し、ステップＳ３２においてリアクタンス値テーブルメモリ１３からｍ番目のリアクタンス値セットを読み出してアレーアンテナ装置１００の各可変リアクタンス素子１２−１乃至１２−９（１２−５を除く。）８に設定する。このとき、ステップＳ３３において受信信号ｙｍ（ｔ）を受信して電波到来方向コンピュータ２０内の一時記憶メモリに格納し、ステップＳ３４においてｍ≧８であるか否かが判断され、ＮＯであるときは、ステップＳ３５で放射パターンパラメータｍを１だけインクリメントした後、ステップＳ３２に戻る。一方、ステップＳ３４でＹＥＳであるときは、ステップＳ３６において一時記憶メモリに格納された９個の受信信号ｙｍ（ｔ）を数３のように受信信号ベクトルとして一時記憶メモリに格納する。次いで、ステップＳ３７において受信信号ベクトルに基づいて式（８ｂ）を用いて相関行列Ｒｙを計算し、固有値分解法を用いて受信信号に対する９個の固有値及び９個の固有ベクトルを計算し、ステップＳ３８において計算された雑音の固有値に基づいて数９を用いて平均化した雑音電力の推定値σ２を計算する。さらに、ステップＳ３９において相関行列Ｒｙと雑音電力の推定値σ２に基づいて式（８ａ）を用いて素子空間に変換された相関行列Ｒｓを計算し、ステップＳ４０において固有値分解法（式（９ａ）乃至式（９ｃ））を用いてＲｓに対する９個の固有値及び９個の固有ベクトルを計算した後、元のメインルーチンに戻る。
【０１２８】
図６は図３のサブルーチンである到来方向推定処理（ステップＳ３）を示すフローチャートである。
【０１２９】
図６のステップＳ４１において到来信号個数パラメータｋを１に初期化し、ステップＳ４２において式（２４）、（２５ａ）、（２５ｂ）及び（２５ｃ）を用いて行列Ｇｋを計算し、ステップＳ４３においてパラメータｋ≧信号到来信号個数ｄであるか否かが判断され、ＮＯであるときは、ステップＳ４４においてパラメータｋを１だけインクリメントした後、ステップＳ４２に戻る。一方、ステップＳ４３でＹＥＳであるときは、ステップＳ４５において正規化行列Ｗｎ１／２に初期値である単位行列Ｉを設定し、ステップＳ４６において反復回数パラメータｌｏｏｐを１に初期化する。次いで、ステップＳ４７において式（２３）を最小化するベクトルｆ（ｆ＝［ｕ０…ｕｄ｜ｖ０…ｖｄ−１］Ｔ）を公知のＱＲ分解法により計算し、ステップＳ４８において式（２２）により２次元のＭＯＤＥ法の行列ＢＨを求め、式（１７ｂ）により正規化行列Ｗｎを計算し、ステップＳ４９において反復回数パラメータｌｏｏｐ≧予め設定した反復回数パラメータｌｏｏｐｔｈであるか否かを判断し、ＮＯであるときはステップＳ５０において反復回数パラメータｌｏｏｐを１だけインクリメントした後、ステップＳ４７に戻る。一方、ステップＳ４９でＹＥＳであるときは、ステップＳ５１においてベクトルｆを用いて式（１４ａ）、（１４ｂ）によりＵ（ｚ）及びＶ（ｚ）の根を式（１５）を用いて計算する。次いで、ステップＳ５２において、上述された次式
［数式４０］
μｋ（１）＝ｓｉｎθｋｃｏｓφｋ、
［数式４１］
μｋ（２）＝ｓｉｎθｋｓｉｎφｋ
の関係より天頂からの仰角（天頂仰角）θｋ及び方位角φｋ（ｋ＝１，２，…，ｄ）を推定して計算し、ＣＲＴディスプレイ２１に表示して元のメインルーチンに戻る。
【実施例１】
【０１３０】
本発明者らは、本実施形態に係る電波到来方向探知装置を試作してそのシミュレーションを行い、以下のシミュレーション結果を得た。
【０１３１】
以下では、モーメント法による数値解析により、モノポールで構成された９素子方形のアレーアンテナ装置１００を用いたＭＯＤＥ法の基本特性を示す。シミュレーションでは、アンテナ素子長３ｃｍ、素子ワイヤ半径０．５ｍｍ、素子間隔Δｘ＝Δｙ＝３ｃｍとした。また、搬送波周波数ｆｃは２．４５ＧＨｚであり、素子長ほぼλ／４のモノポールアンテナが、λ／４間隔で配列した方形モノポールアレーアンテナに相当する。なお、モーメント法による解析では、無限地板の接地導体１１を仮定した。以下の解析では、電子制御導波器アレーアンテナ装置における特性を正確に把握するため、等価ウエイト行列Ｗの推定（すなわち校正）及び各シミュレーションデータとともに、図１のアレーアンテナ装置１００に平面波が到来する際の受信値をモーメント法により求め、使用している。複数の波を入射する際には、送信電力を一定とした（式（７ｄ）のｓｉは必ずしも一定ではない）。簡単のため、到来波数は既知とし、雑音は付加していない。
【０１３２】
９素子の電子制御導波器アレーアンテナ装置であるため、等価素子空間変換を行うためのリアクタンス行列Ｗのランクは９となる。従って、最低でも９個の独立なリアクタンス値セットが必要となる。ここでは、上述の表１に示す９個のリアクタンス値セットを用いた。
【０１３３】
［数値計算結果：様々な方向での１波到来方向推定］
図７及び図８は、１波入射のＤＯＡ推定結果を重ね書きしたものである。等価素子空間変換を行うためには、表１に示したリアクタンス時の等価ウエイト行列Ｗを正確に推定する必要がある。図７は、天頂仰角θ＝９０°の校正信号を用いて、方位角φ＝０°、３０°、…、３００°、及びφ＝４５°、１３５°、２２５°、３１５°の１６個の校正データから得られた等価ウエイト行列Ｗを用いた推定結果である。φの推定結果は、いずれも正確であるが、天頂仰角θに関しては、校正に用いたθ＝９０°方向（水平面方向）から離れるに従い、大きな偏差（バイアス）を有していることが分かる。図８は校正信号の天頂仰角θを８０°としたものである。天頂仰角θ＝８０°以上の領域において、推定結果が改善されていることが分かる。式（４）から分かるように等価ウエイト行列Ｗの行は、素子間相互結合の関数である。アレーアンテナ装置１００は強い素子間相互結合を利用してビームを形成している。非特許文献１３に示されているように、素子間相互結合行列の等価回路モデルは、微小ダイポールのような素子上の電流分布が一定と見なせる場合には全てのＤＯＡに対して厳密に成立するが、ある程度の素子長になるとシングルモード動作の場合であっても到来方向依存性が表れる。本実施形態で取り扱う電子制御導波器アレーアンテナ装置では方位角φ方向依存性は少ないが、天頂仰角θ方向の影響は無視できないことが分かる。ただし、これらの結果が示しているように、校正信号のＤＯＡを適切に設定すれば、改善可能である。
【０１３４】
［数値計算結果２：インコヒーレント２波推定］
図９及び図１０にインコヒーレント２波推定結果を３例（データ１、データ２、データ３）示す。インコヒーレントであるので、推定誤差のバイアスの傾向は、基本的に１波推定と類似しており、校正波のθ方向に依存した偏差がみられる。
【０１３５】
［数値計算結果３：コヒーレント２波推定］
図１１は図９に対応したコヒーレント２波推定結果である。コヒーレント波の場合、２波が結合した１つの信号部分空間から推定するため、方位角φ方向にも２ないし３°の誤差の増加が認められるが、ＳＳＰを用いることなく分離検出されている。図１２は図１０の校正において方位角φ＝４５°、１３５°、２２５°、３１５°の校正信号のみをθ＝８０°にした推定結果である。データ２に関しては、全体がシフトしているが、データ１、データ３に関しては改善が見られる。全方向に対して偏差をなくすことは困難であるが、校正信号の到来方向を最適化することによる特性の改善が期待できることが分かる。また、ここでは表１のリアクタンス値セットを用いたが、リアクタンス値セットの設定も推定性能に大きな影響を与えることが予想される。
【０１３６】
以上説明したように、等価素子空間変換の概念を用いることにより、電子制御導波器アレーアンテナ装置であるアレーアンテナ装置１００においても様々なアレー信号処理アルゴリズムの適用が可能となる。本実施形態では、一例として、９素子方形モノポールアレーアンテナであるアレーアンテナ装置１００を用いたＭＯＤＥ法のシミュレーション結果を示した。コヒーレント波、インコヒーレント波ともに、水平面方向のＤＯＡ推定は良好な結果が得られた。ただし、仰角が大きく（すなわち天頂仰角θが小さく）なると、偏差が増加する場合がある。
【０１３７】
以上の実施形態においては、９素子の方形の電子制御導波器アレーアンテナ装置を一例として用いているが、本発明はこれに限らず、少なくとも非励振素子Ａ１，Ａ２，Ａ４と、励振素子Ａ５とを備えた最低４素子の方形電子制御導波器アレーアンテナ装置を用いて電波到来方向探知装置を実現できる。すなわち、素子数は４以上の整数であってもよい。
【０１３８】
付録：信号部分空間の固有値、固有ベクトルｅｈｋを、固有値λｋ及び固有ベクトルｅｋを用いて次式のように定義すると、式（１７ａ）の主要部分ＢＨｅｈｋは次式のように変形できる。
【０１３９】
【数１９】
2005062040.tif 000022
【０１４０】
【数２０】
2005062040.tif 000023
【０１４１】
ここで、行列Ｐ，Ｑ，Ｒは以下のような変形により定義される行列である。
【０１４２】
【数２１】
2005062040.tif 000024
【０１４３】
【数２２】
2005062040.tif 000025
【０１４４】
【数２３】
2005062040.tif 000026
【図面の簡単な説明】
【０１４５】
【図１】本発明の実施形態に係る電波到来方向探知装置の構成を示すブロック図である。
【図２】図１のアレーアンテナ装置１００の詳細な構成を示す断面図である。
【図３】図１の電波到来方向探知装置の電波到来方向探知コンピュータ２０によって実行されるメインルーチンである電波到来方向推定処理を示すフローチャートである。
【図４】図３のサブルーチンである等価ウエイト行列を校正して計算する処理（ステップＳ１）を示すフローチャートである。
【図５】図３のサブルーチンである固有値分解処理（ステップＳ２）を示すフローチャートである。
【図６】図３のサブルーチンである到来方向推定処理（ステップＳ３）を示すフローチャートである。
【図７】本実施形態に係る電波到来方向探知装置のシミュレーション結果であって、様々な方向での１波到来方向推定結果（天頂仰角θ＝９０度の校正信号を利用；設定値○及び推定値×を含む。）を示すグラフである。
【図８】本実施形態に係る電波到来方向探知装置のシミュレーション結果であって、様々な方向での１波到来方向推定結果（天頂仰角θ＝８０度の校正信号を利用；設定値○及び推定値×を含む。）を示すグラフである。
【図９】本実施形態に係る電波到来方向探知装置のシミュレーション結果であって、インコヒーレント２波の推定結果（天頂仰角θ＝９０度の校正信号を利用；設定値○及び推定値×を含む）を示すグラフである。
【図１０】本実施形態に係る電波到来方向探知装置のシミュレーション結果であって、インコヒーレント２波の推定結果（天頂仰角θ＝８０度の校正信号を利用；設定値○及び推定値×を含む）を示すグラフである。
【図１１】本実施形態に係る電波到来方向探知装置のシミュレーション結果であって、コヒーレント２波の推定結果（天頂仰角θ＝９０度の校正信号を利用；設定値○及び推定値×を含む）を示すグラフである。
【図１２】本実施形態に係る電波到来方向探知装置のシミュレーション結果であって、コヒーレント２波の推定結果（天頂仰角θ＝８０度の校正信号を利用；設定値○及び推定値×を含む）を示すグラフである。
【符号の説明】
【０１４６】
Ａ５…励振素子、
Ａ１，Ａ２，Ａ３，Ａ４，Ａ６，Ａ７，Ａ８，Ａ９…非励振素子、
１…低雑音増幅器（ＬＮＡ）、
２…ダウンコンバータ、
３…Ａ／Ｄ変換器、
４…無線受信機、
９…同軸ケーブル、
１０…リアクタンス値コントローラ、
１１…接地導体、
１２−１，１２−２，１２−３，１２−４，１２−６，１２−７，１２−８，１２−９…可変リアクタンス素子、
１３…リアクタンス値テーブルメモリ、
２０…電波到来方向探知コンピュータ、
２１…ＣＲＴディスプレイ、
３０…回転機構、
１００…アレーアンテナ装置。
【図１】
2005062040.tif 000027
【図２】
2005062040.tif 000028
【図３】
2005062040.tif 000029
【図４】
2005062040.tif 000030
【図５】
2005062040.tif 000031
【図６】
2005062040.tif 000032
【図７】
2005062040.tif 000033
【図８】
2005062040.tif 000034
【図９】
2005062040.tif 000035
【図１０】
2005062040.tif 000036
【図１１】
2005062040.tif 000037
【図１２】
2005062040.tif 000038