(57)【要約】  本発明は、１つまたは複数の電波または音波の到来方向を判断する方法とシステムに関し、より詳しくは、利用できる観察結果の数が少ない、あるいは利用できる観察結果を受信するときの強度が低い状況で特に有利な方法とシステムに関する。この方法は、サンプルサイズが小さい場合に、従来のサブスペース信号処理アルゴリズムの性能を大幅に改善する。アルゴリズムは、観察結果の次元が観察回数と同程度の大きさであっても、一貫した推定が行えるように特に設計されている。これにより、センサまたはアンテナの数とサンプル数が同程度の大きさである有限サンプルサイズシナリオ、あるいは受信した信号の強度が従来の一次検索方式による検出を確実に行うのには不十分なシナリオにおいて、確実に良好な挙動が見られる。
2009529660.tif 000002
【特許請求の範囲】
【請求項１】
Ｍ個（Ｍは自然数であり、Ｍ＞１である）の受信素子（１，２，３，４…）によって少なくとも１つの波形（ｗａ）を受信し、前記Ｍ個の受信素子（１，２，３，４…）の各々の出力において１つの検出された信号（ｄ１，ｄ２，ｄ３，ｄ４…）を取得するステップと、
Ｍ個の前処理ユニット（３０，３１，３２，３３…）において前記Ｍ個の検出された信号（ｄ１，ｄ２，ｄ３，ｄ４…）を処理し、Ｍ個の処理済み信号（ｆ１，ｆ２，ｆ３，ｆ４…）を取得するステップと、
Ｎ回（Ｎは自然数である）の時刻において、信号処理ユニット（４０）で前記Ｍ個の処理済み信号（ｆ１，ｆ２，ｆ３，ｆ４…）をサンプリングし、対応するサンプル値を取得するステップと、
前のステップで取得した前記サンプル値から、Ｎ個の列ベクトルｙ（１），…，ｙ（Ｎ）（ただし、ｙ（ｎ）は、あるサンプリング時刻ｎについて、Ｍ行の列ベクトルを表し、前記Ｍ行の各々はそのサンプリング時刻ｎにおいて対応する受信素子（１，２，３，４…）に関連付けられた前記サンプル値を有する）を作るステップと、
前記Ｎ個の列ベクトルの集合ｙ（１），…，ｙ（Ｎ）から、サンプル空間相関行列Ｒ＾
を計算するステップと、
前記サンプル空間相関行列Ｒ＾から、すべての固有ベクトルｅ＾ｋ（ただし、ｋ＝１…Ｍ）とすべての固有値
【数１】
2009529660.tif 000025
を取得するステップと、を含む、少なくとも１つの波形（ｗａ）の到来方向を検出するための方法であって、
さらに、
前のステップにおいて取得した前記すべての固有ベクトルｅ＾ｋと前記すべての固有値
【数２】
2009529660.tif 000026
から、信号依存パラメータμｋ（ただし、ｋ＝１…Ｍ）の集合を取得するステップと、
すでに取得した前記固有値
【数３】
2009529660.tif 000027
と前記信号依存パラメータμｋから、重みの集合を構成するステップと、
特定のコスト関数から、特定の要求事項を満足するようなＫ個の解を求めるステップと、
を含むことを特徴とする方法。
【請求項２】
請求項１に記載の方法であって、
Ｎ＞Ｍであり、ＮとＭは自然数であることを特徴とする方法。
【請求項３】
請求項に記載の方法であって、
Ｎ＝Ｍであり、ＮとＭは自然数であることを特徴とする方法。
【請求項４】
請求項１に記載の方法であって、
Ｎ＜Ｍであり、ＮとＭは自然数であることを特徴とする方法。
【請求項５】
請求項１から４のいずれかに記載の方法であって、
前記少なくとも１つの波形（ｗａ）を受信する前記ステップは、Ｍ個のセンサまたはアンテナ素子のアレイ（１，２，３，４…）によって行われることを特徴とする方法。
【請求項６】
請求項１から５のいずれかに記載の方法であって、
Ｋ個の波源があり、その各々が波形を発生することを特徴とする方法。
【請求項７】
請求項１に記載の方法であって、さらに、
前記信号処理ユニット（４０）において、列ベクトルａ（θ）の数値を取得し、保存するステップを含み、前記列ベクトルａ（θ）はＭ個の受信素子の各々に対応するＭ本の行を有し、前記Ｍ本の行の各々は、探査対象の到来方向θに応じた関数を有することを特徴とする方法。
【請求項８】
請求項７に記載の方法であって、さらに、
前記Ｍ個の受信素子のアレイを、
前記Ｍ個の受信素子のアレイに、前記方向θから到来する波形をシミュレートする参照信号の集合を印加し、
前記信号処理ユニット（４０）において、前記Ｍ個の素子の各々で取得される前記受信した参照信号の応答を測定してａ（θ）を取得する
ことによって校正するステップを含むことを特徴とする方法。
【請求項９】
請求項１に記載の方法であって、さらに、
前記Ｍ個の受信素子のアレイ（１，２，３，４…）を、
前記Ｍ個の受信素子のアレイ（１，２，３，４…）に参照信号の集合を印加し、
前記信号処理ユニット（４０）において、前記Ｍ個の前処理ユニット（３０，３１，３２，３３，…）の前記出力での前記受信した参照信号の応答を測定し、
前記Ｍ個の前処理ユニット（３０，３１，３２，３３，…）の前記応答を、前記Ｍ個の前処理ユニット（３０，３１，３２，３３，…）の前記出力で測定された前記信号が、前記第一のステップで挿入された前記参照信号にできるだけ近くなるように変調する
ことによって校正するステップを含むことを特徴とする方法。
【請求項１０】
請求項１から９のいずれかに記載の方法であって、
前記信号依存パラメータの集合μｋを取得する前記ステップは、次式
【数４】
2009529660.tif 000028
のμのすべての個別の解μ１≦…不等号μＭを取得することによって行われることを特徴とする方法。
【請求項１１】
請求項１０に記載の方法であって、
Ｎ＜Ｍの場合、μ１＝…＝μＭ−Ｎ＝０であり、μＭ−Ｎ＋１≦…≦μＭは前記式のμの正の解であることを特徴とする方法。
【請求項１２】
請求項１から１１のいずれかに記載の方法であって、
前記重みの集合φｍを構成するステップは、次式、
【数５】
2009529660.tif 000029
を適用することによって行われ、上式において、
Ｋは波源の数、したがって前記Ｍ個の受信ユニット（１，２，３，４…）において受信される波形（ｗａ，ｗｂ，ｗｃ…）の数であり、
【数６】
2009529660.tif 000030
、ただし、
【数７】
2009529660.tif 000031
であることを特徴とする方法。
【請求項１３】
請求項１から１２のいずれかに記載の方法であって、
前記サンプル空間相関行列Ｒ＾は、
【数８】
2009529660.tif 000032
の形態をとることを特徴とする方法。
【請求項１４】
請求項１から１３のいずれかに記載の方法であって、
前記サンプル空間相関行列Ｒ＾は、当業者に周知の他のどの形態でもよいことを特徴とする方法。
【請求項１５】
請求項１２に記載の方法であって、請求項７に従属するとき、
前記コスト関数から前記Ｋ個の解を求めるステップは、以下のコスト関数、
ｆ（η（θ））
（ただし、ｆ（・）は実数値増加関数である）のＫ個の最深極小値を選択するか、あるいは以下のコスト関数、
ｇ（η（θ））
（ただし、ｇ（・）は実数値減少関数である）のＫ個の最高極大値を選択することによって行われ、
【数９】
2009529660.tif 000033
であることを特徴とする方法。
【請求項１６】
請求項１２に記載の方法であって、Ｍ個の受信素子（１，２，３，４…）が等間隔線形アレイを形成するとき、前記方法はさらに、
行列
【数１０】
2009529660.tif 000034
を構成するステップと、
多項式
【数１１】
2009529660.tif 000035
（ただし、
【数１２】
2009529660.tif 000036
であり、（・）＊は複素共役を表し、Γｍ，ｎは行列Γのｍ番目の行とｎ番目の列に対応する成分を表す）を構成するステップと、
複素多項式Ｄ（ｚ）の根から、数量ｆ（｜Ｄ（ｅｘｐ（ｊａｒｇ（ｚｋ）））｜）（ただし、ｆ（・）は実数値増大関数である）を極小化するか、あるいは数量ｇ（｜Ｄ（ｅｘｐ（ｊａｒｇ（ｘｋ）））｜）（ただし、ｇ（・）は実数値減少関数である）を極大化するものとして選択するステップと（両関数においてｊは虚数単位であり、ａｒｇ（・）は複素数の偏角を表し、｜｜は複素数のモジュールを表す）、
ｚ１，…，ｚＫからＫ個の推定到来方向、
【数１３】
2009529660.tif 000037
（ただし、ｄは前記アレイの前記素子の間の距離を表し、λは衝突波形の波長である）を取得するステップと、
を含むことを特徴とする方法。
【請求項１７】
請求項１２に記載の方法であって、請求項７に従属するとき、
前記Ｍ個の受信素子（１，２，３，４…）が等間隔線形アレイを形成するとき、前記方法はさらに、
前記Ｍ個の素子のアレイを、それぞれＬ個の受信素子を含むＭｓｓ個のより小さなサブアレイ（ただし、ＭｓｓとＬは自然数である）に分割するステップと、
前記Ｎ個の列ベクトルｙ（１），…，ｙ（Ｎ）の各々を、それぞれＬ本の行を含むＭｓｓ個のより小さな列ベクトルｙＭｓｓ（１），…，ｙＭｓｓ（Ｎ）に分割するステップと、
前記Ｍｓｓ個のより小さなサブアレイの各々においてＭｓｓ個のサンプル空間相関行列を計算するステップと、
前記Ｍｓｓサンプル空間相関行列を平均し、そこから平均Ｌ×Ｌの空間相関行列Ｒ＾を構成するステップと、
ａ（θ）の最初のＬ行を選択することにより、列ベクトルａｓｓ（θ）を構成するステップと、
コスト関数
ｆ（η（θ））
（ただし、ｆ（・）は実数値増加関数である）のＫ個の最深極小値を選択するか、あるいはコスト関数
ｇ（η（θ））
（ただし、ｇ（・）は実数値減少関数である）のＫ個の最高極大値を選択する（さらに、両関数において、
【数１４】
2009529660.tif 000038
である）ステップと、
を含むことを特徴とする方法。
【請求項１８】
請求項１２に記載の方法であって、請求項７に従属するとき、
特定のコスト関数からＫ個の解を求める前記ステップは、Ｋ回の繰返しの工程によって前記コスト関数のＫ個の極値を選択することによって行われ、ｋ回目の反復、ｋ＝１…Ｋにおいて、ｋ番目の極値は、コスト関数ｆ（θｋ（θ））（ただし、ｆ（・）は実数値増加関数である）の対応する最小値か、コスト関数ｆ（θｋ（θ））（ただし、ｇ（・）は実数値減少関数である）の対応する最大値のいずれかとして選択され、
両関数において、
【数１５】
2009529660.tif 000039
であり、またθ１，…，θｋ−１は前回の反復において選択された極値であり、‖・‖はベクトルのユークリッドノルムを表すことを特徴とする方法。
【請求項１９】
請求項１２に記載の方法であって、さらに、
前記信号処理ユニット（４０）で得られた前記サンプル値にＭｂ×Ｍの行列Ｂ（ただし、Ｍ＜Ｍｂである）を乗じることにより、到来方向の選択された領域を通じてアルゴリズムの分解能を改善するステップを含み、前記行列Ｂの数値は前記受信素子の前記アレイ配置に依存することを特徴とする方法。
【請求項２０】
請求項１９に記載の方法であって、
前記サンプル空間相関行列Ｒ＾は、
【数１６】
2009529660.tif 000040
の形態をとることを特徴とする方法。
【請求項２１】
請求項２０に記載の方法であって、
処理対象となる信号の到来方向は、コスト関数
ｆ（ηＢ（θ））
（ただし、ｆ（・）は実数値増加関数である）のＫ個の最深極小値を選択するか、あるいはコスト関数
ｇ（ηＢ（θ））
（ただし、ｇ（・）は実数値減少関数である）のＫ個の最高極大値を選択することによって推定され、両関数において、
【数１７】
2009529660.tif 000041
であることを特徴とする方法。
【請求項２２】
請求項１から２１のいずれかに記載の方法であって、
前記少なくとも１つの受信波形（ｗａ）は無線周波数信号またはマイクロ波信号であることを特徴とする方法。
【請求項２３】
請求項１に記載の方法の前記ステップを実行するように構成された手段を備えることを特徴とするシステム。
【請求項２４】
請求項２３に記載のシステムであって、
前記受信ユニット（１，２，３，４…）は、Ｍ個のセンサまたはアンテナ素子のアレイであることを特徴とするシステム。
【請求項２５】
請求項２３または２４のいずれかに記載のシステムであって、
前記Ｍ個の前処理ユニット（３０，３１，３２，３３…）は無線周波数または音響フロントエンド（３０，３１，３２，３３…）であることを特徴とするシステム。
【請求項２６】
請求項２３から２５のいずれかに記載のシステムであって、
前記少なくとも１つの受信された波形（ｗａ）は、少なくとも１つの波源によって生成されることを特徴とするシステム。
【請求項２７】
請求項２６に記載のシステムであって、
前記少なくとも１つの波源は送信機であることを特徴とするシステム。
【請求項２８】
請求項２６に記載のシステムであって、
前記少なくとも１つの波源は、他の波源により生成された音波または電波の反射であることを特徴とするシステム。
【請求項２９】
請求項２３から２８のいずれかに記載のシステムであって、
前記少なくとも１つの受信波形（ｗａ）は無線周波数信号またはマイクロ波信号であることを特徴とするシステム。
【請求項３０】
コンピュータプログラムであって、
コンピュータ、デジタル信号プロセッサ、フィールドプログラマブルゲートアレイ、特定用途集積回路、マイクロプロセッサあるいはその他の形態のプログラム可能なハードウェア上で実行されたときに請求項１に記載の方法の前記ステップを実行するように構成されたコンピュータプログラムコード手段でなることを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は1つまたは複数の電波または音波の到来方向を判断するための方法とシステムに関し、より詳しくは、利用できる観察結果が少ない、あるいは利用できる観察結果の受信時の強度が低い状況において特に有利な上記の方法とシステムに関する。
【背景技術】
【０００２】
多くの分野において、特定の受信機から見た、1つまたは複数の電波または音波の到来方向（ＤｏＡ）あるいはまた到来角度を判断することに対して高い関心が寄せられている。衝突角度の推定は、特定の領域に複数のセンサまたはアンテナ素子を空間的に分散させたアレイを用いることで可能になる。センサの各々は、無線周波数／音響フロントエンドおよびこれに対応する信号処理技術を備える。アレイ素子の各々が受信する波形または信号は、中央処理ユニットによって一緒に処理され、中央処理ユニットはそのシナリオでの異なる波源(source)の到来方向を判断する。波源は、それ自体が送信機かもしれないし、他の信号源から発せられた音波または電波の反射かもしれない。
【０００３】
たとえば、レーダの分野においては、ある基準点から見た、異なる標的の正確な位置を推定することが有益である。これは、特定の無線周波数波形を送信し、シナリオ内の異なる物体から跳ね返された後の遅延と信号が受信される到来方向を推定することによって行うことができる。アンテナアレイを用いれば、三角化法に頼らずに、標的の正確な位置を推定できる。
【０００４】
その一方で、移動通信等の他の利用分野の中には、移動局等の信号源の位置を推定することに関心が寄せられるものもある。前述のように、これは、アンテナアレイを受信機として使用することによって有効に行われる。基地局に向けて同時に送信している異なる移動局の到来方向は、本願で説明する技術を用いて推定できる。
【０００５】
個別の波形発生源の到来方向（ＤｏＡ）を判断するための方法は非常に多くあり、中で最も重要な貢献を果たしているものとしては、ケイポン推定法［Ｊ．ケイポン: ”High resolution frequency-wavenumber spectrum analysis”, Proceedings of the IEEE, vol. 57, pp. 1408-1418, Aug. 1969］、最尤（ＭＬ）推定法［Ｐ．ストイカおよびＡ．ネホライ: ”Performance study of conditional and unconditional direction-of-arrival estimation”, IEEE Transactions on ASSP, vol. 38, pp. 1783-1795, Oct. 1990］、ＭＵＳＩＣ（“ＭＵｌｔｉｐｌｅ  ＳＩｇｎａｌ  Ｃｌａｓｓｉｆｉｃａｉｏｎ”）と呼ばれる超解像法［Ｒ．シュミット: ”Multiple emitter localization and signal parameter estimation”, Proceedings of the RACD, Spectral Estimation Workshop, Rome (NY), pp. 243-258, 1979. “Modern Spect5ral Analysis II”, S.B.ケスラー（編集者）, IEEE Press, New York, 1986に転載］がある。
【０００６】
サンプル数が十分に多いと、ＭＬ法では最良の推定性能を発揮するが、それと引き換えに、計算は格段に複雑となる（ＭＬの解は、複数の可変値の中のコスト関数を最適化することによって計算され、この解は多次元検索から得られる）。ケイポン法とＭＵＳＩＣスペクトルのどちらによる推定方式も、解を一次元検索から得ることから、計算上の負荷の点でははるかに効率的である。これら２つの方法の違いは、ＭＵＳＩＣ法が、ケイポン法による解決策と異なり、観察結果の空間相関行列の固有構造を利用することである。現実には、理論的にも実験的にも、固有構造方式に基づくＤｏＡ推定法は、ケイポン推定方式等、他の従来の一次元検索に基づく方法より優れていることが示されている。
【０００７】
これに対し、Ｖ．ギルコはＧ推定(G-estimation)の理論を発表した［Ｖ．ギルコ: ”An Introduction to Statistical Analysis of Random Arrays”, The Netherlands: VSP, 1998］。この理論は、サンプルサイズが小さい場合に優れた検出特性を有する推定量を得るための体系的方法を提供する。
【０００８】
ＭＵＳＩＣ等の部分空間に基づく方法の主な欠点は、相互の間隔の狭い標的を検出する能力が、利用できるサンプルの数か受信する波形の強度のいずれかが特定の基準値以下になったときに大幅に低下する点である。強度が低い、またはサンプルサイズが小さい領域におけるＭＵＳＩＣ方式の検出性能の低下を緩和するための方法がいくつか提案されているが、その多くはある特定のアレイ形状にしか適用できない。たとえば、多項式の根(polynomial-root)を利用するＭＵＳＩＣ方式または状態空間方式は、等間隔の線形アレイアーキテクチャ（つまり、すべてのセンサ／アンテナが直線に沿って等間隔で配置されているアレイ上）以外では運用できない。これに対し、文献では、ＭＵＳＩＣ方式に重み付けを行った方法が検出性能低下の閾値を超えて有効性を改善する上で役立つ可能性があると指摘されている。しかしながら、現在のところ、一般的な状況におけるＭＵＳＩＣアルゴリズムの優れた特性を失うことなく、強度が低い、またはサンプルサイズの小さい状況でこの方法の有効性を最適化するためにこれらの重みをどのように選択するかについては、明確な指標がない。
【発明の開示】
【課題を解決するための手段】
【０００９】
本発明の主目的は、ＭＵＳＩＣの検出性能低下効果を克服し、利用できるサンプル数が少ないか、あるいは受信した波形の強度が、一次元検索方式によるＤｏＡ検出を確実に行うのに十分な大きさではない状況で、近接させて配置された標的の到来方向（以下、ＤｏＡという）を適正に検出することを可能にする方法とシステムを提供することである。本発明はＧ推定の理論の教示を利用しているため、本願による方法を以下、ＧＭＵＳＩＣ方式と呼ぶ。本願による方法はさらに、多項式の根、ビームスペースプロジェクション(beamspace projection)、逐次推定または空間平均等、従来のＭＵＳＩＣ方式の一定の改良版をさらに補完することができる。
【００１０】
本発明の目的物は、マルチセンサまたはマルチアンテナアレイ受信機のＤｏＡ推定のための新規な方法とシステムである。本発明の主な特徴は、固有構造による従来の方法の閾値効果を克服し、利用できるサンプルの数または受信した波形の強度が、到来角度の一次探査に基づく従来の方法によるＤｏＡ検出を実行するのに不十分である状況で、複数のＤｏＡを正しく検出することである。
【００１１】
本願で提案するＤｏＡ特定方法は、探査角度に依存するコスト関数の最低および最高ピークを選択するステップを含む。この関数は、観察結果の推定共分散行列の固有値と固有ベクトルから構成される。
【００１２】
ＭＵＳＩＣアルゴリズムの周知の改良版と比較した本発明による方法の別の利点は、線形アレイ、非線形アレイ、等間隔アレイ、円形アレイ、三角アレイ等、どのようなタイプの受信素子のアレイであっても有効である点である。
【００１３】
本願の提案による、少なくとも１つの波形の到来方向を判断する方法は、Ｍ個の受信素子で少なくとも１つの波形を受信するステップと、Ｍ個の受信素子の各々の出力における少なくとも１つの検出信号を取得するステップと、上記Ｍ個の検出信号を処理して、Ｍ個の処理済信号を得るステップと、Ｎ回の時刻に上記Ｍ個の処理済信号をサンプリングするステップと、前のステップで取得したサンプルから、Ｎ個の列ベクトルｙ（１），…，ｙ（Ｎ）を作るステップ（ただし、ｙ（ｎ）は、あるサンプリング時刻ｎについて、Ｍ行の列ベクトルを表し、Ｍ行の各々はそのサンプリング時刻ｎにおける対応する受信素子に関連付けられるサンプル値を含む）と、Ｎ個の列ベクトルの集合から、サンプル空間相関行列Ｒ＾を計算するステップと、サンプル空間相関行列Ｒ＾からすべての固有ベクトルｅ＾ｉとすべての固有値λ＾ｉを取得するステップと、すべての固有ベクトルｅ＾ｉとすべての固有値λ＾ｉから、信号依存パラメータμκの集合を取得するステップと、固有値λ＾ｉ
と信号依存パラメータμκから、重みの集合を構成するステップと、特定のコスト関数から、特定の要求事項を満たすようなＫ個の解を求めるスップと、を含む。
【００１４】
上記の特定のコスト関数からＫ個の解を求めるステップは、特定のコスト関数のＫ個の最深極小値(deepest local minima)を選択するか、あるいは特定のコスト関数のＫ個の最高極大値(highest local maxima)を選択することによって行われる。
【００１５】
本発明はまた、上記の方法の各ステップを実行するように適応された手段を備えるシステムならびに、コンピュータ、デジタル信号プロセッサ、フィールドプログラマブルゲートアレイ、特定用途集積回路、マイクロプロセッサ、マイクロコントローラ、あるいはその他の形態のプログラム可能なハードウェアの上で実行された時に上記の方法の各ステップを実行するように適応されたコンピュータプログラムコード手段でなるコンピュータプログラムを提供する。
【００１６】
本発明の利点は、以下の説明から明らかとなる。
【００１７】
説明を十分なものとし、本発明をよりよく理解できるように、一連の図面を添付する。この図面は明細書の一部であり、本発明の好ましい実施例を示しており、これは本発明の範囲を限定するものではなく、本発明の実施方法の例として解釈されるべきである。図面は後述の図で構成される。
【発明を実施するための最良の形態】
【００１８】
本発明の実現は、図１に示す装置によって達成される。１つまたはいくつかの波形（ｗａ，ｗｂ，ｗｃ…）は、少なくとも１つの受信素子（１，２，３，４…）によって受信される。これらの波形は好ましくは無線周波数またはマイクロ波信号である。受信素子は、好ましくは、複数のセンサまたはアンテナで構成されるセンサアレイである。１つまたはいくつかの波形（ｗａ，ｗｂ，ｗｃ…）がそれぞれ１つまたは複数の波源から供給され、この波源はそれ自体が送信機であってもよく、あるいは供給される波形は他の波源から発生される音波または電波が反射されたものであってもよい。受信素子（１，２，３，４…）の各々で捕捉または検出される波形または信号（ｄ１，ｄ２，ｄ３，ｄ４…）は、前処理ユニットまたはフロントエンド（３０，３１，３２，３３…）が受信し、これらが信号を前処理し（ｆ１，ｆ２，ｆ３，ｆ４…）、信号処理ユニット（４０）によるその後の処理に備える。信号処理ユニット（４０）は、特定の同期された時刻においてフロントエンド（３０，３１，３２，３３…）から供給された信号（ｆ１，ｆ２，ｆ３，ｆ４…）のサンプルをとり、個別の時間の多次元観察結果またはスナップショットの集合［（ｆ１１，ｆ１２，ｆ１３，ｆ１４…），（ｆ２１，ｆ２２，ｆ２３，ｆ２４…），…，（ｆＮ１，ｆＮ２，ｆＮ３，ｆＮ４…）］（ただし、Ｎは（ｆ１，ｆ２，ｆ３，ｆ４…）のサンプルがとられる時刻の数を表す自然数である）を生成する。
【００１９】
Ｍ個の受信素子、センサまたはアンテナ（１，２，…，Ｍ）が利用可能である、つまり、センサまたはアンテナアレイがＭ個の受信素子で形成されるとすると、各サンプリング時刻ｎにおいて、信号処理ユニット（４０）は、Ｍ個の異なるフロントエンド（３０，３１，３２，３３，…）から供給されるＭ個の数値の集合（ｆｎ１，ｆｎ２，…，ｆｎＭ）を取得する。これらの数値（ｆｎ１，ｆｎ２，…，ｆｎＭ）は、受信機の特徴に応じて実数値または複素数値のいずれでもよい。Ｍ個の数値はＭ×１の列ベクトルｙ（ｎ）（ただし、ｎはその観察結果と関連するサンプリング時刻を示す自然数である）に蓄積される。したがって、受信された信号または波形がＮ個の異なる時刻において同調してサンプリングされた場合、信号処理ユニット（４０）はＮ個の列ベクトルまたはスナップショットの集合｛ｙ（１），…，ｙ（Ｎ）｝を処理する。利用できるスナップショットの数Ｎは、センサの数Ｍより少なくとも、これと同じでも、これより多くてもよい。
【００２０】
Κ個の、たとえば送信機、標的または他の波源から発生された音波あるいは電波の反射等の波源があるとすると、この状況において、上記の列ベクトルの各々は、次のように表すことができる。
ｙ（ｎ）＝Ａｓ（ｎ）＋ｎ（ｎ）                      （１）
ただし、
・ｓ（ｎ）は、Κ個の異なる波源または反射体により送信または反射された信号のサンプルを含むΚ×１の列ベクトルである。
・ｎ（ｎ）は、ｎ番目の時刻においてフロントエンド（３０，３１，３２，３３…）により生成された（熱）騒音サンプルを含むＭ×１の列ベクトルである。
・Ａは列ａ（θ１），…，ａ（θΚ）（ただし、ａ（θ）は角度θから送られる信号に関連するアレイ応答であり、θ１，…，θΚは、その状況に存在するΚ個の指向性波源のＤｏＡである）のＭ×Ｋ行列である。アレイ応答ａ（θ）は、Ｍ×１の列ベクトルであり、その成分(entry)は到来角度の関数である。各波形はＭ個の受信素子のセンサまたはアンテナアレイにおいて、異なる到来角度θ（つまり到来方向）で到来するため、式（１）により表される信号においては、Κ個の波源の各々が角度θ１，…，θΚでの関数ａ（θ）を評価することによって得られる列ベクトルＭ×１とともに寄与する。
【００２１】
これらの関数は既知であるか、あるいは信号処理ユニット（４０）において他の手段により取得されると仮定される。この列ベクトルａ（θ）の構成方法を説明するために、Ｍ個の素子の等間隔の線形アレイを考え、平面波形がアレイのブロードサイド(broadside)に関してθの角度でアレイ上に衝突すると仮定する。第一のアンテナで受信される信号をｓ（ｎ）とする（つまり、送信された信号であり、説明を簡単にするために、遅延は考慮しない）。すると、この信号の帯域幅が十分に低いと、第二のアンテナが受信する信号はｓ（ｎ）ｅｘｐ（ｊ２πｄ／λ・ｓｉｎ（θ））（ただし、ｊは虚数単位、ｄはアレイ素子間分離、λは衝突する波形の波長である）となる。さらに、ｍ番目のアンテナが受信する信号は、ｓ（ｎ）ｅｘｐ（ｊ２π（ｍ−１）ｄ／λ・ｓｉｎ（θ））となる。したがって、Ｍ×１の列ベクトルの各アンテナで受信される信号を集めると、この列ベクトルをｓ（ｎ）・ａ（θ）（ただし、前述のように、ｓ（ｎ）は送信された信号、ａ（θ）はｍ番目の成分がｅｘｐ（ｊ２π（ｍ−１）ｄ／λ・ｓｉｎ（θ））に等しいＭ×１である）として表現できることがわかる。一般に、列ベクトルａ（θ）は、アレイの配置に合わせて構成でき、到来方向θの複合関数としてアレイ応答を行わせる。
【００２２】
すでに述べたように、これらの関数ａ（θ１）…ａ（θΚ）は既知の、あるいは信号処理ユニット（４）において他の手段により得られるものと仮定される。
【００２３】
原則として、受信連鎖(receiving chain)が理想的なら、受信素子のアレイの配置がわかっていれば、ａ（θ）がわかる。しかしながら、実際には受信連鎖は理想的ではなく、特低の量の歪みが生じる。このため、このタイプの受信素子のアレイの構成においては、関数ａ（θ）の形態を求めるステップを実施するか、あるいは受信連鎖の応答を補正して、所望の関数ａ（θ）を得ることが普通である。この工程は時々、「アレイ校正」と呼ばれる。関数ａ（θ）の形態を求めるために、通常、参照信号が通常アンテナまたはセンサに導入され、これらの参照信号が方向θから発せられる波形をシミュレートする。次に、受信連鎖の応答が信号処理ユニット（４０）で測定され、この測定の結果がａ（θ）と呼ばれる。基本的に、この校正工程は、探査マージン(exploration margin)内のどのθの数値についても実行する必要があるため、θの各数値について考えられるすべての関数ａ（θ）の集合を求める。しかしながら、校正工程は探査そのものの前に行われるため、現実には、θのどの数値についてもａ（θ）はわからない（θの無限の数値について実験を繰り返さなければならないであろう）。このため、現実には、θのできるだけ多くの、またできるだけ相互に近接した個別の数値の集合を選択し、校正工程中に関数ａ（θ）を評価し、把握する。
【００２４】
後に、到来方向の検出アルゴリズムの実行中に、コスト関数η（θ）は通常、θの数値のうち関数ａ（θ）がわかっているものについて、θの残りの数値を補間することによって評価される。
【００２５】
平面アレイに関するかぎり、つまり受信素子が平面上で分散されるアレイであるかぎり、関数ａ（θ）は１つの角度（線形アレイ、つまり受信素子が直線上に分散される場合、これに当てはまる）ではなく、２の角度（仰角と方位角）に依存する。
【００２６】
すなわち、θはパラメータかパラメータのベクトルのいずれかとすることができ、その数値は、到来方向を決定する角度である。
【００２７】
Ｍ個の受信素子のアレイ（１，２，３，４…）は次の方法でも校正できる。参照信号の集合をＭ個の受信素子のアレイ（１，２，３，４…）に印加するし、信号処理ユニット（４０）においてＭ個の前処理ユニット（３０，３１，３２，３３，…）の出力での受信参照信号の応答を測定し、Ｍ個の前処理ユニット（３０，３１，３２，３３，…）の応答を、Ｍ個の前処理ユニット（３０，３１，３２，３３，…）の出力で測定された信号が最初のステップで挿入された参照信号とできるだけ近くなるように変調する。
【００２８】
図２は、観察結果｛ｙ（１），…，ｙ（Ｎ）｝からの到来方向θ１，…，θΚを取得するための好ましい方法を示している。提案された方法は、「背景技術」に記載のＧ推定量の理論の教示を利用することから、ＧＭＵＳＩＣと呼ばれる。従来のＭＵＳＩＣアルゴリズムと反対に、ＧＭＵＳＩＣはサンプルサイズとセンサまたはアンテナの数の両方が同程度の大きさであるときに、一貫したＤｏＡ推定が行われるように設計されている。これにより、強度が低い、またはサンプルサイズが小さい場合の推定性能を大きく改善することができる。
【００２９】
ＭＵＳＩＣアルゴリズムでは、Ｎが大きいがＭが一定のとき、つまり、ＮがＭと相対して大きい場合に到来方向または到来角度の推定結果は一貫している。これに対し、本発明によるアルゴリズムによれば、ＮとＭの両方が高い、つまり、ＭとＮが大きさにおいて相互に同等であるときに、一貫した信頼性の高い結果が得られる。現実において、これは、本発明によるアルゴリズムが、サンプルの数Ｎがセンサまたはアンテナの数Ｍよりはるかに多いという条件下でなくても、きわめて高い信頼性で挙動することを意味する。これは、もともとのＭＵＳＩＣアルゴリズムでは実現できない。したがって、本発明によるアルゴリズムは、「サンプルサイズが小さい場合」において、ＭＵＳＩＣアルゴリズムよりずっと有効に機能する。
【００３０】
上記と同じ論理は、強度の弱い場合にも当てはまる。到来方向検出またはスペクトル推定のどのアルゴリズムにおいても、検出対象の信号の強度が大きければ大きいほど、信号のＤｏＡの推定を高い信頼度で実現するのに必要なサンプル数は少なくてすむ。本発明によるアルゴリズムを用いれば、ＭＵＳＩＣアルゴリズムと比較して、少ないサンプル数でも良好な検出が可能となるため、同じサンプル数の場合に、ＭＵＳＩＣアルゴリズムがその信号検出のために必要とするであろう信号強度より弱い強度の信号を検出することができる。
【００３１】
この方法は、次のステップを含む。
１．Ｎ個の列ベクトルまたはスナップショットの集合｛ｙ（１），…，ｙ（Ｎ）｝から、サンプル空間相関行列またはサンプル共分散行列Ｒ＾を計算する。
好ましい実施例において、このサンプル空間相関行列は次の形となる。
【数１】
2009529660.tif 000003
ただし、（・）Ｈは、実数値の信号を扱っている場合は移項、複素数値信号を扱っている場合は共役転置を意味する。しかしながら、このサンプル空間相関行列は、（２）の形態には限定されず、サンプル空間相関行列Ｒ＾は、当業界で知られている他のいずれの形態でもよい。
【００３２】
２．行列Ｒ＾の固有値と固有ベクトルを得る。これらはそれぞれ、｛ｅ＾ｋ，ｋ＝１…Ｍ｝と
【数２】
2009529660.tif 000004
と表される。
【００３３】
３．次式のμの個別の解を得る。
【数３】
2009529660.tif 000005
これらの解は、μ１≦…≦μＭの順序となる。
Ｍ≦Ｎの場合、これらＭ個の解は基本的にすべて異なる。λ＾ｋの数値のいくつかが異ならない場合、つまり、これらがリスト｛λ＾１，…，λ＾Ｍ｝の中で繰り返される場合、μｋもリスト｛μ１，…，μＭ｝の中で、λ＾ｋがリスト｛λ＾１，…，λ＾Ｍ｝で繰り返されるのと同じ回数だけ繰り返される。
Ｍ＞Ｎの場合、μ１＝…＝μＭ−Ｎ＝０と設定し、μＭ−Ｎ＋１≦…≦μＭが上記の式のμに正の解を含むようにする。これらの解は、原則としてすべて異なる。前述のように、λ＾ｋの数値のいくつかが異ならない場合、つまり、リスト｛λ＾Ｍ−Ｎ＋１，…，λ＾Ｍ｝において繰り返される場合、数値μｋもまた、λ＾ｋがリスト｛λ＾Ｍ−Ｎ＋１，…，λ＾Ｍ｝の中で繰り返されるのと同じ回数だけ｛μＭ−Ｎ＋１，…，μＭ｝において繰り返される。
【００３４】
４．次の重みを構成する。
【数４】
2009529660.tif 000006
【数５】
2009529660.tif 000007
【数６】
2009529660.tif 000008
ただし、
【数７】
2009529660.tif 000009
である。
【００３５】
５．推定されたＤｏＡは、コスト関数ｆ（η（θ））（ただし、ｆ（・）は実際値の増加関数である）のＫ個の最深極小値を選択するか、あるいはコスト関数ｇ（η（θ））（ただし、ｇ（・）は実際値の減少関数である）のＫ個の極大値を選択することによって選ばれ、
【数８】
2009529660.tif 000010
である。
【００３６】
関数、
【数９】
2009529660.tif 000011
のＫ個の最深極小値として推定されたＤｏＡを選択する、従来のＭＵＳＩＣ方式との相似性に注意されたい。
【００３７】
本発明は、（従来のＭＵＳＩＣ方式の場合のように）最小固有値に関連するものだけではく、すべての固有ベクトルからの情報を利用する。次に、強度が低い、またはサンプルサイズが小さい状況でのＭＵＳＩＣ方式の検出性能の低下の問題を本発明が克服する点について説明する。
【００３８】
図３は、２つの点波源が方位角３５°と４０°に位置づけられ、その信号がガウス型ホワイトノイズで受信される場合の、上記（７），（８）の関数の特定の例を到来角度の関数として示したグラフである。波源は、Ｍ＝１０個の素子が衝突波形の波長の半分に等しい素子間分離で線形に配置された、センサまたはアンテナの等間隔線形アレイを照明する。この例では、合計Ｎ＝７０個のスナップショットが受信機において利用でき、信号は、背景騒音の強度より１０ｄＢ高い強度で受信された。この特定の実験において、従来のＭＵＳＩＣ方式では２つの指向性波源を区別できず、最深極小値は３９．５１°と−３２．８２°の位置と特定された。本願の提案によるＧＭＵＳＩＣの発明を用いると、システムは３４．９２°と３９．５８°で２つの角度波源を認識でき、これは波源の位置の実際の数値に非常に近い。
【００３９】
図４は、ＭＵＳＩＣ方式の提案中のものと従来のものを５００回実行した結果を平均することにより得られた平均二乗誤差を比較したものであり、そのアレイについてＮ＝７とＮ＝７０のスナップショットが使用できた。シナリオは図３と同じである。結果は、異なる信号対騒音比の数値について示され、ＳＮＲは受信した信号の強度（２つの波源で同じとする）と背景全方向性騒音の強度との商である。参考として、無条件のクラーメル・ラオの不等式の下限(Cramer-Rao Bound)（ＣＲＢ）、つまり、到来方向の不偏推定量で実現可能な最低分散も示している。ＧＭＵＳＩＣ方式が示すＳＮＲの閾値は、従来のＭＵＳＩＣよりはるかに低い点に注意されたい。さらに、提案された方式は、高いＳＮＲ領域でもよりよい推定性能を示している。
【００４０】
図５は、従来のＭＵＳＩＣとＧＭＵＳＩＣシステムの検出確率を、１０４回の実施結果の平均から得たものとして示す。検出確率は、
【数１０】
2009529660.tif 000012
（ただし、θ１とθ２はシミュレートされたシナリオでの２つの点波源のＤｏＡである）として定義される。特に、θ１＝０°に固定し、θ２は０°から１０°の範囲で可変とした。結果は、２つのサンプルサイズ、Ｎ＝７，Ｎ＝７０と２つのＳＮＲの数値、つまりＳＮＲ＝１０ｄＢ，ＳＮＲ＝４０ｄＢについて示されている。ＧＭＵＳＩＣ方式はまた、２つの間隔の狭い波源の検出確率の点で、従来のＭＵＳＩＣより優れていることが実験的に実証された。
【００４１】
本明細書ですでに述べたように、本発明による方法はさらに、多項式根計算、ビームスペースプロジェクション、逐次推定または空間平均等、従来のＭＵＳＩＣ方式を改良したものにおいても実現でき、これらについて以下の好ましい実施例に示す。
【００４２】
上記のように、本願の方法が等間隔線形アレイまたは受信素子に適用された場合、ＤｏＡは、前述の本発明による方法を若干改変したものを使って以下のように推定できる。この手順は通常、多項式根計算と呼ばれる。
【００４３】
Ｍ×Ｍ行列を、
【数１１】
2009529660.tif 000013
と定義し、Γｍ，ｎでｍ番目の行とｎ番目の列に対応する成分を表すものとする。
【００４４】
複素数多項式、
【数１２】
2009529660.tif 000014
で、
【数１３】
2009529660.tif 000015
（ただし、（・）＊は複素共役である）
の場合を考える。
【００４５】
ＤｏＡは、次の数値として推定される。
【数１４】
2009529660.tif 000016
上式において、ｄはアレイの素子間距離を表し、λは衝突波形の波長であり、ａｒｇ（・）は複素数の偏角を表し、ｚ１…ｚΚは、｛ｚ（ｋ），ｋ＝１…２（Ｍ−１）｝として表される複素数多項式Ｄ（ｚ）の２（Ｍ−１）累乗根から、数量ｆ（｜Ｄ（ｅｘｐ（ｊａｒｇ（ｚｋ）））｜）（ただし、ｆ（・）は実数値増加関数である）を最小化するか、あるいは数量ｇ（｜Ｄ（ｅｘｐ（ｊａｒｇ（ｚｋ）））｜）（ただし、ｇ（・）は実数値減少関数である）を最大化するものとして選択され、またｊは虚数単位、ａｒｇ（・）は複素数の偏角を表し、｜｜は複素数のモジュールである。
【００４６】
本明細書において先に提案した方法を受信素子の等間隔線形アレイに適用した場合に取り入れることのできる別の手順は、空間平均法である。この方法は、空間自由度のすべてを空間共分散行列の推定の改善に使用するものである。この手法において、Ｍ個の素子のアレイはＭｓｓ個のより小さなサブレアレイに分割され、その各々がＬ個のセンサまたはアンテナを含む。Ｍｓｓ個のより小さなサブアレイの間に重複がなければ、Ｍｓｓ×Ｌ＝Ｍである。反対に、Ｍｓｓ個のより小さなサブアレイの間に重複があれば、Ｍｓｓ×Ｌ＞Ｍである。Ｎ個の列ベクトルｙ（１），…，ｙ（Ｎ）の各々は、Ｍｓｓ個のより小さな列ベクトルｙＭｓｓ（１），…，ｙＭｓｓ（Ｎ）に分割され、この小さな列ベクトルの各々がＬ本の行を有する。
【００４７】
共分散行列は次のように推定される。まず、Ｍｓｓ個のサブアレイの各々の共分散またはサンプル空間相関行列を（２）における対応するサンプル相関行列として、あるいは当業界で周知のその他の空間相関行列推定値を使って推定する。つまり、Ｍｓｓ個の小さなサブアレイの各々におけるＭｓｓ個のサンプル空間相関行列を計算する。次に、Ｌ×Ｌの共分散行列推定値であるＲ＾を、前のステップで得られたＭｓｓ個の共分散行列を平均することによって構成する。ａｓｓ（θ）を、ａ（θ）の最初のＬ行を選択することによって形成される列ベクトルを表すものとする。この方法は、ａ（θ）の代わりにａｓｓ（θ）と、上記の推定によるサンプル共分散行列Ｒ＾に置き換えるステップを含む。その後、次のコスト関数ｆ（η（θ））のＫ個の最深極小値（ただし、ｆ（・）は実数値増加関数である）か、あるいは次のコスト関数ｇ（η（θ））のＫ個の最高極大値（ただし、ｇ（・）は実数値減少関数である）のいずれかが次式により計算される。
【数１５】
2009529660.tif 000017
【００４８】
前述の本発明によるアルゴリズムを、さらによい結果を導くために改変したまた別の方法は、逐次検出方法に基づく。（７）におけるコスト関数、つまり
【数１６】
2009529660.tif 000018
のＫ個の最深極小値（またはＫ個の最高極大値）を選択する代わりに、Ｋ回の反復工程の後にＤｏＡを推定するものであり、この反復工程では、ｋ回目（ｋ＝１…Ｋ）の反復時において、ｋ番目のＤｏＡが、コスト関数ｆ（θｋ（θ））（ただし、ｆ（・）は実数値増加関数）の対応する最小値またはコスト関数ｇ（θｋ（θ））（ただし、ｇ（・）は実数値減少関数）の対応する最大値のいずれかとして選択され、
【数１７】
2009529660.tif 000019
であり、上式において、
【数１８】
2009529660.tif 000020
Ａｋ＝［ａ（θ１），…，ａ（θｋ−１）］
ここで、θ１，…，θｋ−１はそれ以前の反復において検出されたＤｏＡ、‖・‖はベクトルのユークリッドノルムを表す。
【００４９】
本願で提案する方法はまた、ＤｏＡ推定手順の特性を改善するために、観察結果またはサンプル共分散行列に一次変換が適用された場合にも利用できる。
【００５０】
たとえば、この種の利用分野においては、受信するスナップショットにＭｂ×Ｍの行列Ｂ（ただし、Ｍ＜Ｍｂ）を乗じることが一般的である。この手順は通常、ビームスペースプロジェクションと呼ばれる。行列Ｂは、求めようとするＤｏＡの特定領域を通じたアルゴリズムの分解能を高めるために、アレイ配置を考慮に入れて設計される。つまり、行列Ｂの数値は、受信素子のアレイ配置によって異なる。本発明による方法は、この状況でも運用できるように変更できることは自明である。（２）においてサンプル共分散行列Ｒ＾
を使う代わりに、これをＭｂ×Ｍｂ変換した、
【数１９】
2009529660.tif 000021
を使用しなければならない。
【００５１】
この点から、対象となる信号のＤｏＡは、コスト関数ｆ（ηＢ（θ））（ただし、ｆ（・）は実数値増加関数）のＫ個の最深極小値を選択するか、またはコスト関数ｇ（ηＢ（θ））（ただし、ｇ（・）は実数値減少関数）のＫ個の最高極大値のいずれかを選択するとによって推定され、
【数２０】
2009529660.tif 000022
であり、上式において、
【数２１】
2009529660.tif 000023
はＲ＾Ｂの固有ベクトルを表し、重み
【数２２】
2009529660.tif 000024
はＲ＾をＲ＾Ｂに置き換えたφｍとして構成される。
【００５２】
本文において、「〜を備える(comprise)」とその派生形（「〜を備えている(comprising)」その他)は、排他的な意味ではないと理解すべきであり、つまり、これらの用語は、ここに記載され、定義されたものがさらに別の素子、ステップその他も含む可能性を排除すると解釈すべきではない。
【００５３】
本発明は明らかに、本明細書に記載の具体的な実施例に限定されず、当業者であれば考えられるすべての変更（たとえば、構成要素、構成等の選択に関するものなど）が、付属の特許請求範囲に定義された本発明の一般的範囲の中に包含される。
【図面の簡単な説明】
【００５４】
【図１】本発明の好ましい実施例による、１つまたは複数の波形の到来方向を推定するための受信システムの概略図である。
【図２】本発明の好ましい方法により実行されるステップの機能ブロック図である。
【図３】本願に記載の方法を実践する際に得られる到来方向「ＤｏＡ」抽出に使用されるコスト関数の特定の例のグラフである。
【図４】この方法、あるいは本発明による好ましい実施例を実践することから得られた、それぞれ平均二乗誤差および検出確率の結果を示すグラフである。
【図５】この方法、あるいは本発明による好ましい実施例を実践することから得られた、それぞれ平均二乗誤差および検出確率の結果を示すグラフである。
【図１】
2009529660.tif 000042
【図２】
2009529660.tif 000043
【図３】
2009529660.tif 000044
【図４】
2009529660.tif 000045
【図５】
2009529660.tif 000046
【国際調査報告】
2009529660.tif 000047
2009529660.tif 000048