JP5320792B2

JP5320792B2 - 到来方向推定装置、到来方向推定方法および到来方向推定プログラム

Info

Publication number: JP5320792B2
Application number: JP2008088091A
Authority: JP
Inventors: 和雄白川
Original assignee: Fujitsu Ltd
Current assignee: Fujitsu Ltd
Priority date: 2008-03-28
Filing date: 2008-03-28
Publication date: 2013-10-23
Anticipated expiration: 2028-03-28
Also published as: EP2105758A2; US7912680B2; US20090243933A1; EP2105758B1; EP2105758A3; JP2009243947A

Description

本発明は到来方向推定装置、到来方向推定方法および到来方向推定プログラムに関し、特に複数のアンテナを用い目標からの到来信号が到来する方向を推定する到来方向推定装置、到来方向推定方法および到来方向推定プログラムに関する。

複数のセンサを用いた目標（信号送信源または信号反射物体）からの到来信号のＤＯＡ（Direction Of Arrival：到来方向）推定には、ＤＢＦ（デジタルビームフォーミング法）、ＳＳＭ（Sub Space Method：部分空間法）、ＭＬ（Maximum Likelihood Method：最尤推定法）等が用いられている。

特許文献１には、ＤＯＡ推定を高速かつ高精度に行うため、到来信号のベースバンド信号の相関ベクトルを組み合わせた空間平均共役分散行列を用いてＤＯＡを推定する技術が開示されている。
国際公開第２００６／０６７８６９号パンフレット

しかしながら、特許文献１の技術においてもＤＯＡ推定の計算速度・ロバスト性には改善の余地がある。本件開示の到来方向推定装置、到来方向推定方法および到来方向推定プログラムは、到来方向の推定を高速かつ高精度に行うことが可能な到来方向推定装置、到来方向推定方法および到来方向推定プログラムを提供することを目的とする。

本到来方向推定装置は、Ｎ個のセンサが受信した目標からの到来信号よりｖ_１〜ｖ_Ｎで表されるＮ個のベースバンド信号ベクトルを生成する信号ベクトル生成手段と、１≦ＭかつＭ≦（Ｎ−１）／２なる自然数Ｍを行の次元として優先的に定め、前記信号ベクトルの部分ベクトルｖ_１〜ｖ _Ｎから（Ｎ−Ｍ）×Ｍ次の行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を生成するＨａｎｋｅｌ行列生成手段と、前記行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を用いて行列Ｒを生成し、前記行列ＲをＲ＝［Ｒ_１｜Ｒ_２］^Ｔによって２つの部分行列Ｒ_１およびＲ_２に分割し、前記部分行列Ｒ_１とＲ_２とから前記到来信号の到来方向を推定する推定手段と、を具備し、前記行列Ｒ_ｆ１は前記信号ベクトルの部分ベクトルｖ_１〜ｖ_Ｎ−１を構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、前記行列Ｒ_ｆ２は前記信号ベクトルの部分ベクトルｖ_２〜ｖ_Ｎを構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、前記行列Ｒ_ｂ１は、Ｒ_ｆ１ ^＊を前記信号ベクトルの共役複素部分ベクトルｖ_１ ^＊〜ｖ_Ｎ−１ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ１＝Ｊ_Ｎ−ＭＲ_ｆ１ ^＊Ｊ_Ｍであり、前記行列Ｒ_ｆ２ ^＊は前記信号ベクトルの共役複素部分ベクトルｖ_２ ^＊〜ｖ_Ｎ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ２＝Ｊ_Ｎ−ＭＲ_ｆ２ ^＊Ｊ_Ｍであることを特徴とする。

また、本到来方向推定方法は、Ｎ個のセンサが受信した目標からの到来信号よりｖ_１〜ｖ_Ｎで表されるＮ個のベースバンド信号ベクトルを生成するステップと、１≦ＭかつＭ≦（Ｎ−１）／２なる自然数Ｍを行の次元として優先的に定め、前記信号ベクトルの部分ベクトルｖ_１〜ｖ _Ｎから（Ｎ−Ｍ）×Ｍ次の行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を生成するステップと、前記行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を用いて行列Ｒを生成し、前記行列ＲをＲ＝［Ｒ_１｜Ｒ_２］^Ｔによって２つの部分行列Ｒ_１およびＲ_２に分割し、前記部分行列Ｒ_１とＲ_２とから前記到来信号の到来方向を推定するステップと、を有し、前記行列Ｒ_ｆ１は前記信号ベクトルの部分ベクトルｖ_１〜ｖ_Ｎ−１を構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、前記行列Ｒ_ｂ１は、Ｒ_ｆ１ ^＊を前記信号ベクトルの共役複素部分ベクトルｖ_１ ^＊〜ｖ_Ｎ−１ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ１＝Ｊ_Ｎ−ＭＲ_ｆ１ ^＊Ｊ_Ｍであり、前記行列Ｒ_ｆ２ ^＊は前記信号ベクトルの共役複素部分ベクトルｖ_２ ^＊〜ｖ_Ｎ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ２＝Ｊ_Ｎ−ＭＲ_ｆ２ ^＊Ｊ_Ｍであることを特徴とする方法である。さらに、本到来方向推定プログラムは、コンピュータに上記方法を実行させるプログラムである。

本到来方向推定装置、到来方向推定方法および到来方向推定プログラムによれば、到来方向の推定を高速かつ高精度に行うことが可能となる。

まず、特許文献１に係る技術を用いてＤＯＡ推定を行う場合に、計算速度およびロバスト性を改善する余地が残る要因について説明する。なお、以下で、行列やベクトル後に付けた上添え字^Ｈは複素共役転置を、^Ｔは置換を、^＊は各複素共役を取ることを、そして、行列の後に付けた上添え字^−１は逆行列を示す。また、行列Ｊ_Ｐは反対角要素が１でその他の要素が０であるＰ次反対角単位行列を示し、行列Ｉ_ｐはＰ次単位行列である。０_Ｐ×ＱはＰ×Ｑ次のゼロ行列、０_ＰはＰ次正方ゼロ行列を示す。

図１は、本発明の受信用センサ（例としてアンテナ１０とする）を示す図である。図１では、Ｎ個のセンサとして、アンテナＡ_１〜Ａ_ＮがＸ軸方向の直線上に等しい間隔ｄを隔てて配置されている。これにより、複数のアンテナ１０は等間隔リニアアレーアンテナ（ＵＬＡ）を構成する。アンテナ１０から距離ｒ_ｍ、角度θ_ｍ（Ｙ軸の正方向を角度の０度とし、時計回りを正方向として測る）なる位置に、独立したＭ個の目標１２（信号送信源または信号反射物体、図１では煩雑さを避ける為、１個で代表させている）が存在しているものとする。

さて、説明を行う便宜上、具体的な装置として波長λなる高周波信号をキャリア信号とするレーダを考える。レーダは、受信用アンテナ１０と空間的に隔てた位置に設置された送信用アンテナから目標に向けて探知信号を送信し、これが目標で反射される事によって生成される、所謂、エコー信号をアンテナ１０で受信し、エコー信号と送信信号とをミキシングしてベースバンド信号に復調し、更に、適切な信号処理を施して各目標の位置や速度情報を推定する。よって図１の状況では、Ｍ（≦Ｎ−１）個の目標からのエコー信号ｘ_ｍ（ｔ）がアンテナ１０に到来する。なお、ｍ＝１，．．．，Ｍであり、ｔは時間である。

このとき、ｎ番目のアンテナＡ_ｎからの入力信号を復調して得られた信号ｖ_ｎ（ｔ）は式１で表される。

ただし、

である。ここで、ｘ_ｍ（ｔ）はベースバンド信号、ｎ_ｍ（ｔ）は加法性ガウス雑音、λはキャリア信号である。

信号ｖ_ｎ（ｔ）をベースバンド信号ベクトルとして表すと式３となる

ただし、

である。ｘ（ｔ）とｎ（ｔ）との間に相関が無いものとして、式３からｖ（ｔ）の共分散行列Ｒ_ｖｖを計算すると、式７となる。

ここで、Ｒ_ｘｘはベースバンド信号ベクトルｘ（ｔ）の共分散行列であり、式８で表される。

Ｒ_ｖｖがＤＯＡ推定する場合の対象となる。

ここではレーダについて考えているので、ＵＬＡで受信された信号は、同じ信号源からの送信された信号が各目標によって反射されてきものであり、コヒーレント性が高い為、Ｎ×Ｎ次行列であるＲ_ｖｖの階数（ランク）は１に縮退する。したがって、このままではＲ_ｖｖからＤＯＡを推定することは難しい。

そこで、Ｒ_ｖｖの主対角線方向にＬ×Ｌ次部分行列（但し、Ｌ＜Ｎ）をＮ−Ｌ＋１個取り、これらの部分行列を足し合わせて平均することにより、Ｒ_ｖｖのランクを回復する、前方空間平均法（forward spatial smoothing：ｆｓｓ）が用いられる。また、ＵＬＡの基準点を反転させｆｓｓと同様の操作を行う後方空間平均法（backward spatial smoothing：ｂｓｓ）も用いられる。通常は双方を組み合わせ、ｆｂｓｓ（forward backward spatial smoothing）として用いられる。

Ｒ_ｖｖにｆｂｓｓを適用して得られるＲ_ｖｖ ^ｆｂｓｓと表すことにすれば、例えば、ｃａｐｏｎ法によるＤＯＡ推定には、式９が用いられる。ただし、ａ（θ）はＬ次のモードベクトルであり、式４を構成するベクトル要素ａ（θ_ｍ）と同様の構造を持つ。

ｍｕｓｉｃ（Multiple Signal Classification）法では、Ｒ_ｖｖ ^ｆｂｓｓを式１０のように固有値分解して行列Ｅ_Ｎを求め、式１１によってＤＯＡ推定を行う。

すなわち、式９または式１１において、θをパラメータとするモードベクトルa(θ)=[1,exp(j2παsinθ)],…,exp(j2πα(L-1)sinθ)]^Tを定義し、θを走査しながら式１１を計算することで行列に含まれる角度情報を調べ、ＤＯＡを推定する。ただし、α＝ｄ／λである。

通常、信号共分散行列には式７の様に雑音成分σ^２Ｉが不可避的に含まれるが、その影響を最小限に抑えてＤＯＡの推定精度を向上させる為、特許文献１の技術では以下のような処理を行っている。すなわち、ＵＬＡを構成する各アンテナで受信した信号を復調して得られたベースバンド信号から相関ベクトルを計算し、擬似共分散行列を生成する。次いで、この擬似共分散行列からプロパゲータ行列を生成し、ＤＯＡ推定を行う。以下にこの推定方法について詳しく説明する。

ここでは、簡単のため、ｆｓｓのみについて考える。ベースバンド信号の相関ベクトルｒ_ｖ１を式１２のように計算する。

これ以降、Ｍは、Ｎ個のアンテナで構成されたＵＬＡによって、ＤＯＡを推定する事が可能な最大の目標数を表すものとする。例えば、Ｍは（Ｎ−１）／２以下の最大の自然数である。

次に、ｒ_ｖ１の要素を並べ替えて、Ｒ_ｆ１を式１３のように生成する。

ここで、

であるから、Ｒ_ｆ１はＲ_ｆ１＝ＡＸなる形をしている。すなわち、式１５に示す行列Ｒは、空間平均を適用した共分散行列Ｒ_ｖｖと同じ位相情報を含んでいる。

以下、Ｎ＝５、Ｍ＝２の場合を例に説明する。このとき行列Ｒは式１６となる。

次に、式１７のように、行列ＲをＭ×Ｍ次部分行列Ｒ_１および（Ｎ−２Ｍ）×Ｍ次部分行列Ｒ_２とに分割する。

部分行列Ｒ_１およびＲ_２を用いて、Ｍ×（Ｎ−２Ｍ）次行列Γ＝（Ｒ_１Ｒ_１ ^Ｈ）^−１Ｒ_１Ｒ_２ ^Ｈを生成する。

行列ΓとＮ−２Ｍ次単位行列：Ｉ_Ｎ−２Ｍから（Ｎ−Ｍ）×（Ｎ−２Ｍ）次プロパゲータ行列Π＝［Γ｜−Ｉ_Ｎ−２Ｍ］^Ｔを生成する。

さらに、この行列Πを用いて、例えば、（Ｎ−Ｍ）×（Ｎ−Ｍ）次核行列Ωを、Ω＝Π（Π^ＨΠ）^−１Π^Ｈによって定義し、これとアレイモードベクトルa(θ)=[1,exp(j2παsinθ)],…,exp(j2πα(N-M-1)sinθ)]^Tとを用いて、例えば、角度スペクトラムＰ（θ）を式１８のように定義する。ただし、α＝ｄ／λである。

そして、パラメータθを走査しながら角度スペクトラムＰ（θ）を計算し、そのピーク位置から到来信号の到来角度θ_ｍを得ることができる。または、ｚ＝ｅｘｐ（ｊ２παｓｉｎ（θ））、α＝ｄ／λ、に対し、代数方程式ａ（１／ｚ）^ＴΩａ（ｚ）＝０の解から到来信号の到来角度θ_ｍを得ることができる。

以上のように、特許文献１では、式１５および１６のように行列Ｒを求め、行列Ｒからプロパゲータ行列Πを求め、核行列Ωを求め、到来方向を算出している。

説明を簡単にするため、以下式１６に示した様に、Ｎ＝５、Ｍ＝２の場合を例に取る。プロパゲータ行列Πの主要部分であるΓは式１９のようになる。

ＤＯＡ推定を行うためには、Γを構成するベクトルの位相成分があれば十分である。しかしながら、式１９に掛かっている因子１／｜ｖ_５｜^２は、ＤＯＡ推定に直接関与しないばかりか、計算の高速化を阻害さえする。また実際の計算に於いては、桁落ちや丸め誤差等の原因となり、高精度化を阻害する。

そこで、本到来方向推定装置は、上記因子を含まない行列を用いてＤＯＡ推定を行うことを特徴とする。

１≦ＭかつＭ≦（Ｎ−１）／２なる自然数Ｍを行の次元として優先的に定め、Ｎ−Ｍを列の次元とし、ベースバンド信号ベクトルｖ_１〜ｖ_Ｎの部分ベクトルｖ_１〜ｖ_Ｎ−１を構成要素として並べ、式２０のようなＨａｎｋｅｌ行列Ｒ_ｆ１を生成する。

具体例として、Ｎ＝５およびＭ＝２の場合についてＲ_ｆ１を式２１に示す。

次に、例えば、式２２に示す様に行列Ｒを定義し、このＲをＭ×Ｍ次部分行列Ｒ_１および（Ｎ−２Ｍ）×Ｍ次部分行列Ｒ_２に分割する。

上記部分行列Ｒ_１およびＲ_２から、Ｍ×（Ｎ−２Ｍ）次行列Γ＝（Ｒ_１Ｒ_１ ^Ｈ）^−１Ｒ_１Ｒ_２ ^Ｈを計算すると、式２３のようになる。

式２３の行列には、ＤＯＡ推定に用いられる位相情報が保存されている反面、式１９に含まれていた余分な因子：１／｜ｖ_５｜^２は含まれない。

以下、式１７と同様に、行列ΓとＮ−２Ｍ次元単位行列Ｉ_Ｎ−２Ｍとからプロパゲータ行列Π＝［Γ｜−Ｉ_Ｎ−２Ｍ］^Ｔを生成する。

更に、この行列Πを用いて、例えば、（Ｎ−Ｍ）×（Ｎ−Ｍ）次核行列Ωを、Ω＝Π（Π^ＨΠ）^−１Π^Ｈによって定義し、これとアレイモードベクトルa(θ)=[1,exp(j2παsinθ)],…,exp(j2πα(N-M-1)sinθ)]^Tとを用いて、例えば、角度スペクトラムＰ（θ）を式１８のように定義する。但し、α＝ｄ／λである。

そして、パラメータθを走査しながら角度スペクトラムＰ（θ）を計算し、そのピーク位置から到来信号の到来角度θ_ｍを得ることができる。または、ｚ＝ｅｘｐ（ｊ２παｓｉｎ（θ））、α＝ｄ／λ、として、代数方程式ａ（１／ｚ）^ＴΩａ（ｚ）＝０の解から到来信号の到来角度θ_ｍを得ることができる。

本到来方向推定装置によれば、行列Ｒ_ｆ１として、式１３の代わりに式２１のような行列Ｒ_ｆ１を用いる。これにより式２３のΓは、式１９の１／｜ｖ_５｜^２のような因子を含まない。よって、Γの計算、プロパゲータ行列Πの計算を高速高精度に行うことができる。以下、上記原理に基づく本発明の実施例について説明する。

図２は実施例１に係る到来方向推定装置のブロック図である。図３は到来方向推定装置が行う処理のフローチャートである。図２および図３を参照して説明を行う。図２に於いて、到来方向推定装置１００は、複数のアンテナ１０（アンテナＡ_１〜Ａ_Ｎ）、受信部２０および到来方向推定部３０を含んでいる。複数のアンテナ１０は目標からの到来信号を受信する。受信部２０は受信器２２および信号ベクトル生成手段２４を含んでいる。受信器２２は、例えば、キャリア信号の周波数が７６ＧＨｚである到来信号を受信し、変調方式に応じた復調を行い、適切な帯域のベースバンド信号にダウンコンバートする。そして更に、このベースバンド信号をＡ／Ｄ変換してデジタル信号に変換する（以下、単にベースバンド信号と記す）。信号ベクトル生成手段２４は、各アンテナに対応するベースバンド信号をベクトル列として整形する。以上により、信号ベクトル生成手段２４は式３のベースバンド信号ベクトルｖ_１〜ｖ_Ｎを生成する（図３のステップＳ１０）。

到来方向推定部３０は、到来信号数予測／設定手段３１、Ｈａｎｋｅｌ行列生成手段３２および推定手段４５を含んでおり、ベースバンド信号ベクトルからＤＯＡを推定する。推定手段４５は、プロパゲータ行列生成手段３３、核行列生成手段３４および到来方向算出手段３５を含んでいる。到来方向推定部３０はコンピュータのＣＰＵ（中央演算処理装置）等の演算装置、演算装置に実行させるプログラム、或いはその記憶装置上の実体である。また、受信部２０と到来方向推定部３０との全てが、ソフトウェア無線装置と同様の意味で、演算装置から構成されていても良い。

到来信号数予測／設定手段３１は、到来信号数Ｍを予測、若しくは設定して出力する（図３のステップＳ１２）。到来信号数Ｍを予測する場合は、最尤法に基づくＡＩＣ（Akaike Information Criteria）、ＭＤＬ（Minimum Description Length）、若しくは他の適切な指標を用いることができる。到来信号数予測／設定手段３１は、到来信号数Ｍの予測を行わない場合は、アンテナ数Ｎから自動的に決定される値である［（Ｎ−１）／２］をＭの値としてもよい。なお、この数式において記号［ｘ］は、実数ｘに対して、ｘ以下の最大の自然数を戻す事を意味する。

推定手段４５は、先ず、式３の信号ベクトルのＳＮ比が事前に決定しておいた或る閾値以上かを判定し（ステップＳ１４）、ＳＮ比が閾値より小さい場合には、例えば、相関ベクトルｃｏｎｊ（ｖ_Ｎ）×［ｖ_１，…，ｖ_Ｎ］^Ｔを算出し、これを信号ベクトルとして置き直す（ステップＳ１６）。ここでｃｏｎｊ（ｖ_Ｎ）は、信号ｖ_Ｎの共役複素数を意味する。

Ｈａｎｋｅｌ行列生成手段３２は、式２０のようなＨａｎｋｅｌ行列Ｒ_ｆ１を生成する（ステップＳ１６に於いて信号ベクトルを置き換えた場合は、式１６の様な構造となるが、今後は特に触れない）。プロパゲータ行列生成手段３３では、行列Ｒ＝［Ｒ_ｆ１］とする（ステップＳ１８）。

プロパゲータ行列生成手段３３は、行列Ｒを２つの部分行列Ｒ_１およびＲ_２に分解し、前述のように行列Γ＝（Ｒ_１Ｒ_１ ^Ｈ）^−１Ｒ_１Ｒ_２ ^Ｈを生成し、さらに行列Γと単位行列Ｉ_Ｎ−２Ｍとからプロパゲータ行列Π＝［Γ｜−Ｉ_Ｎ−２Ｍ］^Ｔを生成する（ステップＳ２２）。

核行列生成手段３４は、核行列Ω＝Π（Π^ＨΠ）^−１Π^Ｈを生成する（ステップＳ２４）。

到来方向算出手段３５は、パラメータθを含むアレイモードベクトルa(θ)=[1,exp(j2παsinθ)],…,exp(j2πα(N-M-1)sinθ)]^Tと核行列Ωとを用いて、式２４のような角度スペクトラムＰ（θ）を生成する（ステップＳ２６）。但し、α＝ｄ／λである。到来方向算出手段３５は、θを走査しながら角度スペクトラムＰ（θ）を計算し、角度スペクトラムＰ（θ）のピーク位置から到来信号の到来角度θ_ｍを計算する。または、ｚ＝ｅｘｐ (ｊ２παｓｉｎ（θ）)、α＝ｄ／λ、として、代数方程式ａ（１／ｚ）^ＴΩａ（ｚ）＝０の解から到来信号の到来角度θ_ｍを計算する。到来方向算出手段３５は、到来角度θ_ｍを推定した到来方向として出力する（ステップＳ２８）。

Ｈａｎｋｅｌ行列生成手段３２は、信号ベクトルの他の部分ベクトルｖ_２〜ｖ_Ｎを構成要素とし、式２５のようなＨａｎｋｅｌ行列Ｒ_ｆ２をＲ_ｆ１の代わりに生成し、プロパゲータ行列生成手段３３は、このＲ_ｆ２を用いてプロパゲータ行列Πを生成してもよい。

また、Ｈａｎｋｅｌ行列生成手段３２は、共役複素部分ベクトルｖ_１ ^＊〜ｖ_Ｎ−１ ^＊を構成要素とするＨａｎｋｅｌ行列、Ｒ_ｂ１＝Ｊ_Ｎ−ＭＲ^＊ _ｆ１Ｊ_ＭをＲ_ｆ１の代わりに生成し、プロパゲータ行列生成手段３３は、このＲ_ｂ１を用いてプロパゲータ行列Πを生成してもよい。なお、Ｒ^＊ _ｆ１は、Ｒ_ｆ１の各要素について共役複素数をとった行列である。

一例として、Ｎ＝５かつＭ＝２の場合を説明すれば、Ｒ_ｂ１は式２６のように表される。

次に、行列Ｒを式２７のように定義し、更に、行列ＲをＭ×Ｍ次部分行列Ｒ_１および（Ｎ−２Ｍ）×Ｍ次部分行列Ｒ_２に分割する。

上記部分行列Ｒ_１およびＲ_２から行列Γ＝（Ｒ_１Ｒ_１ ^Ｈ）^−１Ｒ_１Ｒ_２ ^Ｈを計算すると、式２８のようになる。

式２８においても式２３と同様に、ＤＯＡ推定に用いられる位相情報が保存され、式１９の１／｜ｖ_５｜^２のような因子は含まない。よって、Γの計算、プロパゲータ行列Πの計算を高速高精度に行うことができる。

さらに、Ｈａｎｋｅｌ行列生成手段３２は、共役複素部分ベクトルｖ_２ ^＊〜ｖ_Ｎ ^＊を用いてＨａｎｋｅｌ行列Ｒ_ｂ２＝Ｊ_Ｎ−ＭＲ^＊ _ｆ２Ｊ_ＭをＲ_ｆ１の代わりに生成し、プロパゲータ行列生成手段３３は、このＲ_ｂ２を用いてプロパゲータ行列Πを生成してもよい。なお、Ｒ^＊ _ｆ２は、Ｒ_ｆ２の各要素について共役複素数をとった行列である。

プロパゲータ行列生成手段３３は、（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１、またはＲ_ｂ２をＫ個組み合わせ、（Ｎ−Ｍ）×（Ｋ×Ｍ）の行列として行列Ｒを構成してもよい。例えば、Ｋ＝２の場合Ｒ＝［Ｒ_ｆ１｜Ｒ_ｆ２］、Ｒ＝［Ｒ_ｂ１｜Ｒ_ｂ２］、Ｒ＝［Ｒ_ｆ１｜Ｒ_ｂ２］、またはＲ＝［Ｒ_ｆ２｜Ｒ_ｆ２］、Ｋ＝４の場合Ｒ＝［Ｒ_ｆ１｜Ｒ_ｆ２｜Ｒ_ｂ１｜Ｒ_ｂ２］とすることができる。この場合、Ｋの値に応じて行列ＲをＭ×（Ｋ×Ｍ）の部分行列Ｒ_１および（Ｎ−２Ｍ）×（Ｋ×Ｍ）の部分行列Ｒ_２に分割し、Ｒ＝［Ｒ_１｜Ｒ_２］^Ｔとする。

以上のように、プロパゲータ行列生成手段３３は、行列Ｒを行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１、およびＲ_ｂ２の少なくとも１つを用いて生成し、この行列Ｒを適切な次元の２つの部分行列Ｒ_１、Ｒ_２に分割する。

実施例１によれば、到来方向推定装置１００は、Ｈａｎｋｅｌ行列生成手段３２がＲ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を生成する。推定手段４５が、行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を用いて行列Ｒを生成し、行列ＲをＲ＝［Ｒ_１｜Ｒ_２］^Ｔなる２つの部分行列Ｒ_１およびＲ_２に分割し、部分行列Ｒ_１とＲ_２とを用いて前記行列Γ、Π、そしてΩを計算し、到来信号の到来方向を推定する。これにより、式１９の１／｜ｖ_５｜^２のような因子を含まないΓ（式２３）を用いる事で、Γの計算、プロパゲータ行列Πの計算を高速高精度に行うことができる。

図４は、実施例２に係る到来方向推定装置のブロック図である。図２を参照して説明を行えば、到来方向推定装置１００ａは、実施例１の図２と比較し、到来方向推定部３０ａがスケーリング行列生成手段３６を含んでいる。同図の他の構成は、実施例１の図２と同じであるから説明を省略する。

スケーリング行列生成手段３６は、プロパゲータ行列生成手段３３が生成した行列Ｒの部分行列Ｒ_１およびＲ_２を用い、（Ｎ−２Ｍ）×（Ｎ−２Ｍ）次スケーリング行列Λ＝Ｒ_２Ｒ_２ ^Ｈ−Ｒ_２Ｒ_１ ^ＨΓを生成する。核行列生成手段３４は、プロパゲータ行列Πとスケーリング行列Λとから核行列Ω＝Π（Λ）^−１Π^Ｈを生成する。または、Π´＝Π（Π^ＨΠ）^−１／２とし、核行列をΩ＝Π´（Λ）^−１（Π´）^Ｈによって生成してもよい。

核行列生成手段３４は、Π（Λ）^−１Π^Ｈの代わりにΠΠ^Ｈを核行列Ωとして生成してもよい。この場合、角度スペクトラムは式２９となり、

代数方程式は式３０となる。

ただし、ｚ＝ｅｘｐ(ｊ２παｓｉｎ（θ）)、α＝ｄ／λ、に対し、

である。

核行列生成手段３４は、核行列をΩ＝Π［α（Π^ＨΠ）＋（１−α）Λ］^−１Π^Ｈとして生成してもよい。ここで、０≦α≦１である。このように、核行列生成手段３４は、Ｈａｎｋｅｌ行列から導かれる複数の行列を複合させたものとして、核行列を生成してもよい。

ここで、α≠０、１とし、Π^ＨΠとΛとが特異点でないとして、Ａ＝Π^ＨΠ、Ｂ＝Λ、１−α＝βとすると、式３２となる。

以上より上の核行列Ωは、公知の核行列Π（Π^ＨΠ）^−１Π^ＨおよびΠ（Λ）^−１Π^Ｈの加重平均となる。よって、公知の核行列を単独で用いたのでは推定困難な目標についても角度推定が可能となる。

核行列生成手段３４は、核行列Ωの生成を、実施例１ではプロパゲータ行列Πから行い、また、実施例２ではプロパゲータ行列Πとスケーリング行列Λとから行う。

到来方向算出手段３５は、プロパゲータ行列Πと、Ｄ、（Ｂ^ＨＡ^−１Ｂ）、（Ｂ^ＨＢ）、（Ｄ−Ｂ^ＨＢ）、Λ（Π^ＨΠ）等、以下に定める行列Ａ，Ｂ，Ｃ、Ｄを適当に組み合わせた代用行列をスケーリング行列Λの代わりに用いて核行列Ωを生成し、上記実施例１と同様にして、角度スペクトラムまたは代数方程式から到来方向を算出してもよい。なお、Ａ＝Ｒ_１Ｒ_１ ^Ｈ、Ｂ＝Ｒ_１Ｒ_２ ^Ｈ、Ｃ＝Ｒ_２Ｒ_１ ^Ｈ、Ｄ＝Ｒ_２Ｒ_２ ^Ｈである。

例としてＤを用いた角度スペクトラムを式３３、代数方程式を式３４に示す。

である。

さらに、到来方向算出手段３５は、スケーリング行列Λ、その代用行列、またはスケーリング行列Λのノルムを用いて到来方向を算出してもよい。

図５から８は、実施例２をＦＭＣＷ（Frequency Modulated Continuous Wave）レーダに適用し、目標定位を実施した例である（これらの図では、距離の単位はｂｉｎである）。

図５（ａ）から（ｄ）は、目標が２つ存在する場合の角度スペクトラムを角度と距離に対し示した図である。設定した目標の距離は４０ｍ、角度は０°および３°である。図５（ａ）はＦＦＴ−ＤＢＦ法、図５（ｂ）はＦＢＳＳ−ＭＵＳＩＣ法、図５（ｃ）は特許文献１の方法（比較例）、図５（ｄ）は実施例２の方法を用いて計算した角度スペクトラムである。角度スペクトラムが仮定した２つの目標位置に集約している場合、つまり、角度スペクトラムが２つの目標付近で点状に表示されている程、高い精度で目標の到来方向が推定されていることを示している。

図６は、図５（ａ）、（ｃ）および（ｄ）の角度スペクトラムを距離４０ｍの位置で角度軸に平行に切り出したものを重ね描きしたものである。図５および図６より、比較例（特許文献１）および実施例２の方法は、ＦＦＴ−ＤＢＦ法、ＦＢＳＳ−ＭＵＳＩＣ法と比較して高い精度で目標の角度推定を実現している事が分かる。

図７および図８は、空間的に近接した位置に配置された複数の目標に対する定位性能を、５個の目標を仮定して（図中×の点）角度スペクトラムによって評価した例である。モノクロ表示で若干判別しづらいが、等高線の込み合っている箇所が角度スペクトラムのピークである。図７は特許文献１の方法を用いて計算した角度スペクトラムを示し、図８は、実施例２において、核行列をΩ＝Π［α（Π）＋（１−α）Λ］^−１Π^Ｈ（α＝０．５）として計算した角度スペクトラムを示している。図８は図７と比較して、×印の周辺の等高線がより込み合っており、仮定した目標付近に角度スペクトラムのピークがある。このように、核行列として複合された核行列を用いることにより、エコー信号のコヒーレンスが高く、推定が難しい多くの目標を正しく検出することができる。

実施例３は、自動車等の車両に搭載された到来方向推定装置の例である。実施例３に係る到来方向推定装置１００ｂは、到来方向推定装置１００ｂを含む、例えば、ミリ波レーダを搭載した自車に対して脅威のある目標の角度（到来方向）を出力する。目標としては、例えば自車の前を走行する自動車等の車両である。図９は、実施例３に係る到来方向推定装置１００ｂのブロック図であるが、到来方向推定装置１００ｂは、受信部２０と到来方向推定部３０の間に脅威度評価手段４０を含んでいる。脅威度評価手段４０は、距離／速度／粗到来方向推定手段３７、目標選別手段３８およびビーム形成手段３９を含んでいる。

図１０は、実施例３に係る到来方向推定装置１００ｂの処理を示すフローチャートである。以下、図１０を参照しながら説明を行う。例えば、ミリ波信号にＦＭ変調を掛ける為に用いられる信号の１周期分に渡り、受信アンテナ毎にエコー信号のサンプリングを行い、信号ベクトル生成手段２４が対応する信号ベクトルのデータセットを抽出する（ステップＳ３０）。即ち、アンテナ数をＮ、サンプリング数をＱとすると、このデータはＮ×Ｑ次行列（データセット）となる。距離／速度／粗到来方向推定手段３７は、時間領域に沿ってＱ次元ベクトルがＮ組並んだデータセットをそれぞれのアンテナ毎にＦＦＴ（高速フーリエ変換）処理する（ステップＳ３２） − Ｑ次元ベクトルのＦＦＴをＮ回行い、周波数領域のデータセットに変換されたＮ×Ｑ行列を得る。距離／速度／粗到来方向推定手段３７は、ＦＦＴの結果として周波数領域に現れたピークに対応するｂｉｎ位置（周波数インデックスの値と考えても良い）から、各目標と自車との相対距離ｒ_ｍおよび相対速度ｓ_ｍを推定し、更に、各目標（ｂｉｎ）位置に於けるＮ次元データベクトルを取り出して、例えばＦＦＴ−ＤＢＦ（デジタルビームフォーミング）処理を行って各目標の大雑把な角度位置を把握する（ステップＳ３４）。また、推定した相対距離、相対速度、角度位置を、各目標固有の識別情報としてペアリングする。

目標選別手段３８は、各目標の脅威度を、目標の相対距離ｒ_ｍと相対速度ｓ_ｍとに基づいて評価する（ステップＳ３６）。脅威度はｒ_ｍが小さい程、そして、ｓ_ｍが大きな負の値となる程（即ち、近接した位置に存在し、自車に高速で接近してくる目標程）、大きくなる。脅威度は、例えば式３６の評価関数を用いる。

ただし、

ここで、ｓ_ＲＥＦ、ｒ_ＲＥＦは、車両メーカーの安全基準等に従ってあらかじめ設定された所定値であるが、例えばＤＳＲＣ（Dedicated Short Range Communications）等のインフラ装置からの情報を受け、交通状況に応じて動的に変更する様にしても良い。

目標選別手段３８は、最大の脅威度を有する目標の脅威度が、基準値以上か否かを判断する（ステップＳ３８）。判定がＹｅｓの場合、目標選別手段３８は、先に粗到来方向推定手段に於いて推定した目標の角度を参照して、最大脅威度を有する目標が自車レーン等、自車の動線と交差する位置を走行するか否かを判断する（ステップＳ４０）。

ステップＳ３８またはＳ４０で判定がＮｏとなった場合、脅威度推定手段４０は、周波数軸上で定位される各目標の相対距離ｒ_ｍおよび相対速度ｓ_ｍのデータ（ｂｉｎ）を、後段の到来方向推定部３０に引き渡す。到来方向推定部３０は、指定された各ｂｉｎ位置のＮ次元データベクトルに対して、例えばデジタルビームフォーミング等、精度は比較的低いが、計算量の小さい通常の到来方向推定処理を適用して各目標の角度推定を行う（ステップＳ４２）。もちろん、粗到来方向推定手段による推定値を用いても、実施例１または実施例２の方法を適用しても良い。到来方向推定部３０は、各目標の識別情報である、相対距離、相対速度、および角度推定値のそれぞれを１つ前の測定期間のデータとして更新保存する（ステップＳ４４）。更新されたデータは、次回の測定期間に於いて得られたデータに基づいて、脅威度評価手段４０が各目標の脅威度を再評価する際に参照される。

ステップＳ４０において判定がＹｅｓとなった場合、目標選別手段３８は、最大の脅威度を持つと判定された目標の相対距離ｒ_ｍおよび相対速度ｓ_ｍのデータ（ｂｉｎ）に基づいて、対応するＮ次元データベクトルを抽出してビーム形成手段３９に引き渡す。ビーム形成手段３９は、距離／速度／粗到来方向推定手段３威度を持つ目標の角度情報を参照し、当該目標が存在する方向に、ビームスペース法によるビーム形成を行う（ステップＳ４６）。この様にして、脅威度推定手段４０は、脅威度の高い目標に関するデータベクトルをビームスペースデータに変換し、後段の到来方向推定部３０に引き渡す。到来方向推定部３０は、詳細な角度推定が可能な実施例１または実施例２の方法で到来方向の推定を行う（ステップＳ４８）。装置が、新たにデータ測定を開始した場合、ステップＳ３０に戻って同様な処理を繰り返す。

実施例３の脅威度評価手段４０は、信号ベクトルを用い、目標の自車に対する脅威を示す脅威度を評価する。到来方向推定部３０の推定手段４５は脅威度に基づき、目標の到来方向を推定する、または推定しない。このように、実施例３においては、脅威度評価手段４０を導入する事により、自車の安全に対する目標の重要性を脅威度という指標で評価し、各目標の脅威度に応じて到来方向推定部３０と連携して適切な到来方向推定手法を適用する事で、装置としてのレスポンスタイムを最適化できる。

実施例４は、信号品質や目標間の相対離角といった追加指標を導入する事で、より、現実的な運用環境に適した処理を実現する装置である。図１１は、実施例４に係る到来方向推定装置１００ｃのブロック図である。同図に於いて、到来方向推定装置１００ｃは、受信部２０と到来方向推定部３０との間に脅威度評価手段５０を含んでおり、ハードウェアとしての基本的な構成は実施例３と殆ど同じであるから説明を省略する。

図１２は、到来方向処理装置１００ｃの処理を示すフローチャートである。脅威度評価手段５０は、信号ベクトルの抽出（ステップＳ５０）から目標毎に相対距離と相対速度の推定（ステップＳ５３）を行うまで、実施例３のそれと同じ処理を行う。但し、実施例３で行っていたＦＦＴ−ＤＢＦによる粗到来方向推定処理は実施しなくても良い。

脅威度評価手段５０の機能として新たに加わったステップは、脅威度の評価（ステップＳ５６）を実施する前の、信号品質の評価（ステップＳ５４）、一つ前の測定期間に推定された目標緒元データの参照（ステップＳ６６）、そして、ステップＳ６６の参照値の相対離角を基準値との比較（ステップＳ６８）、である。

新たに加わった処理の流れは、先ず、ステップＳ５４で、信号ベクトルのＳＮ比が基準値以上か否かを判断し、判定がＹｅｓの場合、ステップＳ５６に進む。一方、判定がＮｏの場合、ステップＳ６６に進み、一つ前の測定期間に推定された目標緒元データのうち到来方向のデータを参照し、ステップＳ６８に於いて参照値の相対離角と基準値との比較を行う。そして、脅威度評価手段５０は、ステップＳ５４の判定がＹｅｓの場合、または、ステップＳ６８の判定がＮｏの場合、各目標について脅威度の評価を行い（ステップＳ５６）、最大の脅威度を有する目標の脅威度が、基準値以上か否かを判断する（ステップＳ５８）。脅威度の評価は、実施例３のステップＳ３６と同じ方法で行うことができる。

脅威度評価手段５０は、ステップＳ５８の判定がＹｅｓである場合、制御を到来方向推定部３０に渡し、詳細な到来方向の推定を行う（ステップＳ６０）。一方、ステップＳ５８の判定がＮｏ、若しくは、ステップＳ６８の判定がＹｅｓである場合、脅威度評価手段５０は、制御を到来方向推定部３０に渡し、粗い到来方向推定を行う（ステップＳ７０）。このように、脅威度評価手段５０は、信号ベクトルを用い目標の自車に対する脅威を示す脅威度を評価する。到来方向推定部３０の推定手段４５は脅威度に基づき、目標の到来方向の推定を粗い角度分解能で行う、または細かい角度分解能で行う。

以下に、ステップＳ６０で行う処理とステップＳ７０で行う処理との違いを示す。行列Ｒの部分行列Ｒ_１、Ｒ_２を用いて、Ｒ_１Ｒ_１ ^Ｈ＝Ａ、Ｒ_１Ｒ_２ ^Ｈ＝Ｂ、Ｒ_２Ｒ_１ ^Ｈ＝Ｂ^Ｈ、Ｄ＝Ｒ_２Ｒ_２ ^Ｈ＝Ｄとし、ＲＲ^Ｈに逆行列の定理を適用すると、式３８となり、

これから直ちに、式３９が得られる。

ステップＳ６０では、核行列としてΩ＝（ＲＲ^Ｈ）^−１−［（Ｒ_１Ｒ_１ ^Ｈ）^−１｜０_{Ｍ×Ｎ−２Ｍ}；０_{Ｎ−２Ｍ×Ｍ}⊥｜０_Ｎ−２Ｍ］を用い、また、ステップＳ７０では、核行列としてΩ＝（ＲＲ^Ｈ）^−１を用いて到来方向の推定を行う。ステップＳ７０の行列は信号部分空間成分の成分を含んでいる為、推定精度が劣る。しかしながら、ステップＳ７０の行列を計算する途中で、（Ｒ_１Ｒ_１ ^Ｈ）^−１の計算をしておく事は比較的容易であるから、実施例４の方法は推定精度の切り替えが簡単に実現できる。

以上、実施例４によれば、エコー信号の品質（ＳＮ比）が基準値以上であり、かつ目標の脅威度が基準値以上である場合、ステップＳ６０のように到来方向の詳細な推定を行う。一方、目標の脅威度が小さく詳細な位置の把握が必要が無い場合や、エコー信号のＳＮ比が小さく時間を掛けても精度の良い推定が望めない場合、ステップＳ７０のように到来方向の粗い推定を行う。即ち、実施例４は運用環境に応じてより適切な動作 − レスポンス重視、または精度重視 − を選択できる装置である。

以上、本実施例について詳述したが、本発明は係る特定の実施例に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形・変更が可能である。

図１は受信用センサを示す図である。図２は実施例１に係る到来方向推定装置のブロック図である。図３は実施例１に係る到来方向推定装置のフローチャートである。図４は実施例２に係る到来方向推定装置のブロック図である。図５（ａ）から図５（ｄ）は各方法で計算された角度スペクトラムを示す図である。図６は距離が一定のときの角度に対する角度スペクトラムを示した図である。図７は実施例２と比較例の角度スペクトラムを示した図である。図８は実施例２の角度スペクトラムを示した図である。図９は実施例３に係る到来方向推定装置のブロック図である。図１０は実施例３に係る到来方向推定装置のフローチャートである。図１１は実施例４に係る到来方向推定装置のブロック図である。図１２は実施例４に係る到来方向推定装置のフローチャートである。

符号の説明

１０アンテナ
１２目標
２０受信部
３０到来方向推定部
３１到来信号予想／設定手段
３２Ｈａｎｋｅｌ行列生成手段
３３プロパゲータ行列生成手段
３４核行列生成手段
３５到来方向算出手段
３６スケーリング行列生成手段
３７距離／速度／到来方向推定手段
３８目標選別手段
３９ビーム形成手段
４５推定手段
４０、５０脅威度評価手段

Claims

Ｎ個のセンサが受信した目標からの到来信号よりｖ_１〜ｖ_Ｎで表されるＮ個のベースバンド信号ベクトルを生成する信号ベクトル生成手段と、
１≦ＭかつＭ≦（Ｎ−１）／２なる自然数Ｍを行の次元として優先的に定め、前記信号ベクトルの部分ベクトルｖ_１〜ｖ _Ｎから（Ｎ−Ｍ）×Ｍ次の行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を生成するＨａｎｋｅｌ行列生成手段と、
前記行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を用いて行列Ｒを生成し、前記行列ＲをＲ＝［Ｒ_１｜Ｒ_２］^Ｔによって２つの部分行列Ｒ_１およびＲ_２に分割し、前記部分行列Ｒ_１とＲ_２とから前記到来信号の到来方向を推定する推定手段と、
を具備し、
前記行列Ｒ_ｆ１は前記信号ベクトルの部分ベクトルｖ_１〜ｖ_Ｎ−１を構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、
前記行列Ｒ_ｆ２は前記信号ベクトルの部分ベクトルｖ_２〜ｖ_Ｎを構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、
前記行列Ｒ_ｂ１は、Ｒ_ｆ１ ^＊を前記信号ベクトルの共役複素部分ベクトルｖ_１ ^＊〜ｖ_Ｎ−１ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ１＝Ｊ_Ｎ−ＭＲ_ｆ１ ^＊Ｊ_Ｍであり、
前記行列Ｒ_ｆ２ ^＊は前記信号ベクトルの共役複素部分ベクトルｖ_２ ^＊〜ｖ_Ｎ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ２＝Ｊ_Ｎ−ＭＲ_ｆ２ ^＊Ｊ_Ｍであることを特徴とする到来方向推定装置。
前記推定手段は、行列Γ＝（Ｒ_１Ｒ_１ ^Ｈ）^−１Ｒ_１Ｒ_２ ^Ｈとし、Ｉ_Ｎ−２ＭをＮ−２Ｍ次単位行列とし、プロパゲータ行列Π＝［Γ｜−Ｉ_Ｎ−２Ｍ］^Ｔを生成する、またはプロパゲータ行列Π´＝Π（Π^ＨΠ）^−１／２を生成し、前記プロパゲータ行列ΠまたはΠ´から前記到来方向を推定することを特徴とする請求項１記載の到来方向推定装置。
前記推定手段は、前記プロパゲータ行列ΠまたはΠ´を用い、核行列Ω＝Π（Π^ＨΠ）^−１Π^ＨまたはΩ＝Π´（Π´^ＨΠ´）^−１Π´^Ｈを生成し、核行列Ωから前記到来方向を推定することを特徴とする請求項２記載の到来方向推定装置。
前記推定手段は、スケーリング行列Λ＝Ｒ_２Ｒ_２ ^Ｈ−Ｒ_２Ｒ_１ ^ＨΓを生成し、プロパゲータ行列ΠまたはΠ´とスケーリング行列Λとから前記到来方向を推定することを特徴とする請求項２記載の到来方向推定装置。
前記推定手段は、前記プロパゲータ行列ΠまたはΠ´とスケーリング行列Λとから核行列Ω＝Π（Λ）^−１Π^ＨまたはΩ＝Π´（Λ）^−１Π´^Ｈを生成し、核行列Ωから前記到来方向を推定することを特徴とする請求項４記載の到来方向推定装置。
核行列生成手段は、複数の核行列を複合させ、複合した核行列から前記到来方向を推定することを特徴とする請求項３または５記載の到来方向推定装置。
前記信号ベクトルを用い、前記目標の前記到来方向推定装置を搭載した自車に対する脅威を示す脅威度を評価する脅威度評価手段を具備し、
前記推定手段は前記脅威度に基づき、前記目標の到来方向を推定する、または推定しないことを特徴とする請求項１から６のいずれか一項記載の到来方向推定装置。
前記信号ベクトルを用い、前記目標の前記到来方向推定装置を搭載した自車に対する脅威を示す脅威度を評価する脅威度評価手段を具備し、
前記推定手段は前記脅威度に基づき、前記目標の到来方向の推定を粗い角度分解能で行う、または細かい角度分解能で行う、ことを特徴とする請求項１から６のいずれか一項記載の到来方向推定装置。
Ｎ個のセンサが受信した目標からの到来信号よりｖ_１〜ｖ_Ｎで表されるＮ個のベースバンド信号ベクトルを生成するステップと、
１≦ＭかつＭ≦（Ｎ−１）／２なる自然数Ｍを行の次元として優先的に定め、前記信号ベクトルの部分ベクトルｖ_１〜ｖ _Ｎから（Ｎ−Ｍ）×Ｍ次の行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を生成するステップと、
前記行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を用いて行列Ｒを生成し、前記行列ＲをＲ＝［Ｒ_１｜Ｒ_２］^Ｔによって２つの部分行列Ｒ_１およびＲ_２に分割し、前記部分行列Ｒ_１とＲ_２とから前記到来信号の到来方向を推定するステップと、
を有し、
前記行列Ｒ_ｆ１は前記信号ベクトルの部分ベクトルｖ_１〜ｖ_Ｎ−１を構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、
前記行列Ｒ_ｆ２は前記信号ベクトルの部分ベクトルｖ_２〜ｖ_Ｎを構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、
前記行列Ｒ_ｂ１は、Ｒ_ｆ１ ^＊を前記信号ベクトルの共役複素部分ベクトルｖ_１ ^＊〜ｖ_Ｎ−１ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ１＝Ｊ_Ｎ−ＭＲ_ｆ１ ^＊Ｊ_Ｍであり、
前記行列Ｒ_ｆ２ ^＊は前記信号ベクトルの共役複素部分ベクトルｖ_２ ^＊〜ｖ_Ｎ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ２＝Ｊ_Ｎ−ＭＲ_ｆ２ ^＊Ｊ_Ｍであることを特徴とする到来方向推定方法。
Ｎ個のセンサが受信した目標からの到来信号よりｖ_１〜ｖ_Ｎで表されるＮ個のベースバンド信号ベクトルを生成するステップと、
１≦ＭかつＭ≦（Ｎ−１）／２なる自然数Ｍを行の次元として優先的に定め、前記信号ベクトルの部分ベクトルｖ_１〜ｖ _Ｎから（Ｎ−Ｍ）×Ｍ次の行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を生成するステップと、
前記行列Ｒ_ｆ１、Ｒ_ｆ２、Ｒ_ｂ１およびＲ_ｂ２のうち少なくとも１つの行列を用いて行列Ｒを生成し、前記行列ＲをＲ＝［Ｒ_１｜Ｒ_２］^Ｔによって２つの部分行列Ｒ_１およびＲ_２に分割し、前記部分行列Ｒ_１とＲ_２とから前記到来信号の到来方向を推定するステップと、
をコンピュータに実行させ、
前記行列Ｒ_ｆ１は前記信号ベクトルの部分ベクトルｖ_１〜ｖ_Ｎ−１を構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、
前記行列Ｒ_ｆ２は前記信号ベクトルの部分ベクトルｖ_２〜ｖ_Ｎを構成要素とし、

で表される（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列であり、
前記行列Ｒ_ｂ１は、Ｒ_ｆ１ ^＊を前記信号ベクトルの共役複素部分ベクトルｖ_１ ^＊〜ｖ_Ｎ−１ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ１＝Ｊ_Ｎ−ＭＲ_ｆ１ ^＊Ｊ_Ｍであり、
前記行列Ｒ_ｆ２ ^＊は前記信号ベクトルの共役複素部分ベクトルｖ_２ ^＊〜ｖ_Ｎ ^＊を構成要素とする（Ｎ−Ｍ）×Ｍ次のＨａｎｋｅｌ行列、Ｊ_Ｎ−ＭおよびＪ_ＭをそれぞれＮ−Ｍ次およびＭ次反対角単位行列としたとき、Ｒ_ｂ２＝Ｊ_Ｎ−ＭＲ_ｆ２ ^＊Ｊ_Ｍであることを特徴とする到来方向推定プログラム。