空時頻域中欠定混合條件下的波達方向估計*

朱立為，汪亞，王翔，黃知濤(.國防科技大學電子科學與工程學院，湖南長沙40073；.電子信息系統復雜電磁環境效應國家重點實驗室，河南洛陽47003)

空時頻域中欠定混合條件下的波達方向估計*

朱立為1，2，汪亞1，王翔2，黃知濤2
(1.國防科技大學電子科學與工程學院，湖南長沙410073；
2.電子信息系統復雜電磁環境效應國家重點實驗室，河南洛陽471003)

波達方向估計是陣列信號處理領域的熱點問題，但經典的波達方向估計方法通常要求陣元數大于源信號個數，即滿足超定條件，而在實際中往往面臨的是源信號個數大于陣元數的欠定條件。基于此，提出了一種基于空間時頻分布的多重信號分類擴展算法，通過將空間時頻分布矩陣進行擴展，實現了欠定條件下的波達方向估計。相比時頻多重信號分類算法，所提算法能同時適應超定和欠定條件；相比已有的欠定波達方向估計方法，其不但保證了波達方向估計的精度，而且放寬了對源信號稀疏性的要求，同時還降低了對快拍數的要求。仿真實驗結果證明了該方法的有效性。

波達方向；時頻分布；欠定

在陣列信號處理領域，波達方向(Direction Of Arrival，DOA)估計是聲吶、雷達、地震學的一個重要研究課題。傳統的DOA估計方法，如多重信號分類(MUltiple SIgnal Classification，MUSIC)方法［1］、改進MUSIC方法［2-3］、時頻MUSIC方法［4］、極大似然方法［5］、空間譜估計方法［6］、波束成形方法［7］等，一般要求信號環境滿足超定條件，即接收天線的陣元個數要大于潛在的源信號個數。但在實際的信號環境下，這一條件并不總能滿足。例如，在機載或星載的非合作通信應用中，由于地面雷達、通信等各種輻射源的大量使用，并且使用的頻段不斷擴展、相互重疊，加上各種自然輻射產生的無意干擾，使得機載或星載接收設備面臨時域高度密集、頻域嚴重混疊、空間相互交錯的復雜電磁環境，而實際上機載或星載設備本身受體積限制，陣元數不能隨意增加，且陣元數目越多，接收設備結構越復雜，價格也越昂貴。因此需要利用有限的接收陣元對盡可能多的源信號的DOA進行估計，即需要研究可以適用于欠定混合條件下的DOA估計方法。

近年來，相繼出現了一些適應欠定混合條件的DOA估計算法。主要可分為兩類:一是基于高階累積量的欠定DOA估計算法［8-9］。文獻［9］利用四階累積量代替協方差矩陣，將MUSIC算法擴展成4-MUSIC算法(4-MUSIC)。為了進一步提高陣列對多個潛在源信號的DOA估計能力，文獻［10］將四階累積量推廣到任意偶數階累積量，提出了2q-MUSIC(q＞1)算法，并詳細分析了采用不同陣列結構等價的虛擬陣元數目。雖然這些高階累積量算法可以通過擴展產生虛擬的陣元以適應欠定混合條件，但是高階累積量的計算過程復雜且需要較多的樣本數目才能保證算法的估計精度。二是其他一些基于時頻域稀疏性的DOA估計算法，如退化分解估計技術(Degenerate Unmixing Estimation Technique，DUET)等也能適應欠定混合條件［11-15］，但這些方法對信號的時頻域稀疏性要求過于嚴格，即任意時頻點只有一個源信號起主導作用，其余源信號取值為0［16］，這在實際應用中并不能完全滿足。

本文提出了基于空間時頻分布(Spatial Time-Frequency Distributions，STFD)的MUSIC算法(STFD-MUSIC)，通過擴展空間時頻混合矩陣［17］來實現欠定條件下的DOA估計。該方法并不需要假設信號是稀疏的，且能夠降低算法對快拍數的要求。

1 算法分析

假設N個窄帶遠場信號S(t)=［s1(t)，…，sN(t)］T入射到由M個陣元組成的天線陣上，觀測信號可以寫為:

其中，X(t)=［x1(t)，…，xM(t)］T和V(t)=［v1(t)，…，vM(t)］T分別表示M個陣元的輸出即觀測信號和噪聲，θ=［θ1，…，θN］代表各源信號的入射角度，混合矩陣A(θ)=［a1(θ1)，…，aN(θN)］表示傳輸矩陣，作為線性矩陣，其K列可以表示為:

本文研究的欠定條件指的是陣元數小于源信號數目，即N＞M；研究目的是估計θ。

1.1 空間時頻分布基本理論

本文算法的核心思路是利用觀測信號的空間時頻分布矩陣實現陣列的虛擬擴展，實現欠定混合條件下的DOA估計。

空間時頻分布矩陣定義如式(3)、式(4):

聯立式(4)和式(5)，可得接收信號的空間時頻分布矩陣與源信號的空間時頻分布矩陣之間的關系，其表示如式(6):

由文獻［18］易知，接收信號的空間時頻分布矩陣的對角元素稱為自源點(auto-terms)，也稱自項，其對應著X(t)中代表各個信號的向量平方；而反對角元素稱為互源點(cross-terms)，對應著多個信號向量的線性組合。

1.2 算法假設條件

為了能夠在時頻域上估計出源信號的DOA，采用了Wigner-Vill分布。同時，假設源信號滿足以下條件:

假設1:在時頻平面上信號的自源時頻點與互源時頻點幾乎是不混疊的；

根據假設2可知本文算法在天線陣元數一定時可接收處理的最大源信號個數。

1.3 算法描述

1.3.1 自源點選擇

對于任意自源點(t，f)∈Ωs(Ω表示整個時頻平面)，接收信號X(t)的空間時頻分布矩陣可以表示為:

式中，DˉSS(t，f)表示對角矩陣，其表示為DˉSS(t，f)=diag［Ds1s1(t，f)，…，DsNsN(t，f)］。

先找出自源點，同時抑制互源點，則可以實現對各個源信號的DOA估計。利用的準則如式(8)所示:

在超定混合條件下，直接對DXX(t，f)進行奇異值分解可以很容易地估計出源信號的DOA。但在欠定混合條件下，不能直接分解來估計源信號的DOA。下面將介紹一種能適用于欠定混合條件下基于空間時頻分布的DOA估計算法。1.3.2基于空間時頻分布的擴展矩陣構造

根據式(8)所述的準則，求出自源點的個數，設在式(8)準則下找出K個自源點，則K個自源點對應的時頻分布矩陣表示如式(9)所示:

定義矩陣D∈K×M，其中(D)km=(Dk)mm(k =1，…，K，m=1，…，M)，則C可以表示成如式(11)所示的形式:

式中，A⊙A*=［a1θ1)，…，aN(θN)?θN)］，其中*，⊙和?分別表示復數的轉置、ri-Rao乘積和Kronecker乘積。因此，要解決的問題就是當N＞M時，如何利用矩陣C來估計角度θ。

一般來說，假設N≤min(M2，K)。根據文獻［19］中的證明可知，當潛在源信號個數與陣元數滿足假設2的關系式時，則A⊙A*和D是列滿秩的。而A⊙A*和D是列滿秩的，意味著C也是列滿秩的，并且C的秩就等于源信號個數N。

1.3.3 子空間求解與角度估計

與傳統的MUSIC方法類似，可以直接對自相關陣RC=CCH進行奇異值分解來獲得噪聲子空間。在本文算法中，雖然信號子空間和噪聲子空間可以通過任意一個RCi進行奇異值分解求得，但是當RCi奇異性較強時，子空間的估計精度會很差，進而影響DOA的估計結果。因此，提出對擴展自相關矩陣集合{RCi|i=1，…，K}進行聯合對角化的方法求解信號子空間和噪聲子空間。相比單一矩陣的處理方法，利用多個矩陣的估計方法可以同時保證算法的估計精度和信噪比適應能力。采用聯合對角化的方法分解自相關陣RC= CCH，其分解后的主要形式如式(12):

式中，RDi=(，Λ為×對角陣。Λ的N個較大對角對應著信號子空間，M2-N個較小對角線元素代表著噪聲子空間。因此，可以利用V矩陣相應的M2-N列向量構成的噪聲子空間來估計源信號的DOA。根據文獻［20］中的研究，可以使用式(13)所示的代價函數進行聯合對角化。

當式(13)的值最小時，可以估計出信號子空間與噪聲子空間。因為RCi是Hermitian陣，所以RCi中組成信號子空間的列向量與組成噪聲子空間的列向量是正交的。詳細來說，Span{VS}= Span{Α⊙Α*}，其中VS是N維信號子空間的基，則Α⊙Α*中所有列向量都與VN中的列向量正交，而且是M2-N維噪聲子空間的基。因此，所有{ai(θi)(θi)，1≤i≤N}向量都與VN中的列向量正交，因此有式(14)所示的關系式:

在實際運用過程中，常常定義空間偽譜的概念來運用上述正交性進行DOA估計，如式(15):

從式(15)可知，只要找出最大“譜線”P(θ)的位置處對應的角度，則這個角度就等于要估計的某一源信號的DOA。

綜上所述，本文所提出的基于空間時頻分布的欠定混合條件下DOA估計的核心算法步驟如下:

(i-1)M+j，k)ij計算拓維后的新矩陣，同時求出C^的協方差矩陣=，然后對協方差矩陣進行聯合對角化，估計出VN。

2 仿真分析

2.1 評價準則

為了驗證算法的性能，仿真時采用平均測向均方根誤差作為指標來衡量DOA的估計性能，定義平均測向均方根誤差(Root-Mean-Square Error，RMSE):

2.2 仿真實驗

仿真分析中，用4個線性調頻信號作為源信號，其頻率分別為:1000MHz，1000.1MHz，1000.4MHz，1000.5MHz，碼速率都為200kbit/s，角度分別為:-43°，-13°，26°，48°。仿真實驗1到4中接收天線是陣元數目為3的均勻線陣，仿真實驗5中接收天線是陣元數目分別為4和5的均勻線陣，相鄰陣元間距為半波長。為降低采樣率和計算量，先將射頻信號的頻率變到中頻，中頻頻率分別為350kHz，450kHz，750kHz和850kHz，采樣率為2MHz。仿真將［-90°，90°］的空域以0.1°間隔均勻采樣，得到離散的假設角度集θ～，進行500次蒙特卡洛試驗。

在仿真實驗中，對本文算法與4-MUSIC算法進行了性能比較，而沒有將本文算法與基于時頻稀疏性的DOA估計算法進行比較，這是因為實驗設置的源信號在時頻域是重疊的，不滿足時頻域稀疏性，基于稀疏性的DOA估計算法無法完成源信號DOA的估計。

在圖1中，仿真實驗給出了信噪比(Signal to Noise Ratio，SNR)為10dB，采樣點數分別為1024和512時，本文算法與4-MUSIC算法的任意10次仿真實驗的空域偽譜估計結果。從圖1(a)可以看出本文算法的空間譜中出現了4個顯著的譜峰，能夠清晰地分辨出四個源信號的DOA；圖1(b)是4-MUSIC算法的空間譜，圖中出現了4個大致的譜峰，勉強能分辨出四個源信號的DOA。當采樣點數降低時，從圖1(c)中可以看出，本文算法仍然可以清晰地分辨出源信號的DOA；而圖1(d)中4-MUSIC算法的空間譜上的譜峰已很不明顯，難以分辨出源信號的DOA。

圖1本文算法與4-MUSIC算法空間偽譜Fig.1 Spatial pseudo-spectrums of the proposed and 4-MUSIC algorithm

圖2 為信噪比為10dB時，本文算法與4-MUSIC算法的平均測向均方根誤差RMSE隨采樣點數變化的曲線(圖中CRLB為克拉美羅下限)。從圖中可看出，隨著采樣點數的增加，平均測向均方根誤差隨之減少；在相同的采樣點數條件下，本文算法的平均測向均方根誤差更小，DOA估計精度更高；而在相同的DOA估計精度下，本文算法需要采樣點數要少，因為在相同的估計性能條件下，相比4-MUSIC的高階累積量的估計，本文算法的空間時頻分布矩陣估計需要的采樣數據更少。

圖2不同采樣點數的RMSE估計結果Fig.2 Estimation results of RMSE for different sample points

圖3 為采樣點數等于1024、信噪比從-5dB到30dB變化時，本文算法與4-MUSIC算法的平均測向均方根誤差RMSE隨信噪比變化曲線。從圖中可以看出，隨著信噪比的變化，本文算法的平均測向均方根誤差更接近CRLB，DOA估計精度更高。當DOA估計精度相同時，本文算法比4-MUSIC算法適應的信噪比更低。

圖3 不同的SNR的估計結果(欠定條件)Fig.3 Estimation results of RMSE for different SNR(underdetermined)

假設陣列誤差矢量ei(1≤i≤2)是相互獨立的高斯隨機變量，即滿足E{eieHj}=σ2δIN，令σ2= 0.000 3，圖4給出了本文算法和4-MUSIC算法的平均測向均方根誤差RMSE隨信噪比變化的曲線。從圖中可以看出，相比4-MUSIC算法，本文算法的陣列誤差適應能力更強。

圖4不同的SNR的估計結果(考慮陣列誤差)Fig.4 Estimation results of RMSE for different SNR(model error)

圖5 為采樣點數設為1024，信噪比從-5dB到30dB變化，陣元數分別為5和6時，本文算法與4-MUSIC算法的平均測向均方根誤差RMSE隨信噪比變化曲線。結合圖3的結果可以看出，本文的DOA估計算法可以有效應用于適定和超定的接收條件，極大地擴展了算法的實際應用范圍，并且隨著陣元數的增加，DOA估計性能也隨之提高。

3 結論

針對欠定混合條件下DOA估計問題，提出了一種基于空間時頻域的估計算法。該方法通過Khatri-Rao乘積和Kronecker乘積對空間時頻混合矩陣進行擴展，以適應欠定混合條件。提出的欠定DOA估計算法中，采用聯合對角化的方法來處理空間時頻矩陣，增強了算法的穩健性，提高了算法在低信噪比條件下的DOA估計能力。仿真結果表明，所提算法的DOA估計性能優于4-MUSIC算法。同時，隨著天線陣元數的增加，所提算法對DOA的估計性能還會大幅度提升。

圖5 不同的SNR的估計結果(接收陣元為5，6)Fig.5 Estimation results of RMSE for different SNR(5 and 6 array elements)

References)

［1］Schmidt R O.Multiple emitter location and signal parameter estimation［J］.IEEE Transactions on Antennas and Propagation，1986，34(3):276-280.

［2］Zhou Q C，Gao H T，Wang F，etal.Modified DOA estimation methods with unknown source number based on projection pretransformation［J］.Progress in Electromagnetics Research B，2012，38:387-403.

［3］Zhou Q C，Gao H T，Wang F.A high resolution DOA estimatingmethod without estimating the number of sources［J］.Progress in Electromagnetics Research C，2012，25:233-247.

［4］Belouchrani A，Amin M G.Time-frequency MUSIC［J］.IEEE Signal Processing，Letters，1999，6(5):109-110.

［5］Ziskind I，Wax M.Maximum likelihood localization of multiple sources by alternating projection［J］.IEEE Transactions on Acoustics Speech＆Signal Processing，1988，36(10):1553-1560.

［6］Malioutov D，Cetin M，Willsky A S.A sparse signal reconstruction perspective for source localization with sensor arrays［J］.IEEE Transactions on Signal Processing，2005，53(8):3010-3022.

［7］Li J，Stoica P.Robust adaptive beamforming［J］.Series in Telecommunications＆Signal Processing，2006，15(7): 2345-2348.

［8］Chevalier P，Albera L，Ferréol A，et al.On the virtual array concept for higher order array processing［J］.IEEE Transactions on Signal Processing，2005，53(4):1254-1271.

［9］Porat B，Friedlander B.Direction finding algorithms based on high-order statistics［J］.IEEE Transactions on Signal Processing，1991，39(9):2016-2024.

［10］Chevalier P，Ferreol A，Albera L.High-resolution direction finding from higher order statistics:the 2-MUSIC algorithm［J］.IEEE Transactions on Signal Processing，2006，54(8):2986 -2997.

［11］Rickard S，Dietrich F.DOA estimation of many W-disjoint orthogonal sources from two mixtures using DUET［C］// Proceedings of the Tenth IEEEWorkshop on Statistical Signal and Array Processing，IEEE，2000:311-314.

［12］Araki S，Sawada H，Mukai R，et al.DOA estimation for multiple sparse sources with normalized observation vector clustering［C］//Proceedings of 2006 IEEE International Conference on Acoustics，Speech and Signal Processing，ICASSP 2006，IEEE，2006.

［13］Zhang W Y，Rao B D.A twomicrophone-based approach for source localization of multiple speech sources［J］.IEEE Transactions on Audio Speech＆Language Processing，2010，18(8):1913-1928.

［14］Zhou Z，Lu S J，Zhang E Y，et al.Underdetermined DOA estimation of LFM signals［C］//Proceedings of 2012 Second International Conference on Instrumentation，Measurement，Computer，Communication and Control(IMCCC)，IEEE，2012:869-872.

［15］Lie JP，Ng B P，See C M S.Multiple UWB emitters DOA estimation employing time hopping spread spectrum［J］.Progress in Electromagnetics Research，2008，78:83-101.

［16］王翔，黃知濤，任嘯天，等.基于時頻單源點檢測和聚類驗證技術的欠定混合盲辨識算法［J］.國防科技大學學報，2013，35(2):69-74.WANG Xiang，HUANG Zhitao，REN Xiaotian，et al.Blind identification of underdetermined mixtures based on detection of time frequency single source point and cluster validation technique［J］.Journal of National University of Defense Technology，2013，35(2):69-74.(in Chinese)

［17］Lu F B，Huang Z T，Jiang W L.Underdetermined blind separation of non-disjoint signals in time-frequency domain based onmatrix diagonalization［J］.Signal Processing，2011，91(7):1568-1577.

［18］Belouchrani A，Amin M G，Thirion-Moreau N，et al.Source separation and localization using time-frequency distributions: an overview［J］.IEEE Signal Processing Magazine，2013，30(6):97-107.

［19］De Lathauwer L.A link between the canonical decomposition in multilinear algebra and simultaneousmatrix diagonalization［J］.Siam Journal on Matrix Analysis and Applications，2006，28 (3):642-666.

［20］Golub G H，Van Loan C F.Matrix computation［M］.2nd ed.USA:the John HopKins University Press，1989.

Underdeterm ined direction of arrival estimation based on spatial time-frequency distributions

ZHU Liwei1，2，WANG Ya1，WANG Xiang2，HUANG Zhitao2
(1.College of Electronic Science and Engineering，National University of Defense Technology，Changsha 410073，China；2.The State Key Laboratory of Complex Electromagnetic Environment Effects on Electronics and Information System，Luoyang 471003，China)

In the field of array signal processing，direction of arrival(DOA)estimation is a hotspot problem.Classical DOA estimation methods usually require the number of sensor should be larger than the source signals’(which the so-called over-determined case is).However，whatwe encounter in practice is always the underdetermined case in which the number of source signal is larger than the sensors’.To solve the problem，amultiple signal classification(MUSIC)extension algorithm based on spatial time-frequency distribution was proposed to achieve the underdetermined DOA estimation by expanding the dimension of the spatial time-frequency distributionsmatrices.Compared with the existing timefrequency MUSIC，the proposed algorithm can be applied to both the over-determined and the underdetermined cases.The proposed algorithm also has advantages over the existing underdetermined DOA estimation methods for it guarantees the estimation precision，relaxes the requirements for source signal sparseness and lowers standards of the number of snapshots.Simulation results confirm the validity and high performance of the proposed algorithm.

direction of arrival；time-frequency distributions；underdetermined

TN911

A

1001-2486(2015)05-149-06

10.11887/j.cn.201505023

http://journal.nudt.edu.cn

2014-09-28

CEMEE國家實驗室開放課題基金資助項目(2014K104B)；國家自然科學基金資助項目(61401490)

朱立為(1985—)，男，湖南郴州人，博士研究生，E-mail:iendwin@163.com；黃知濤(通信作者)，男，教授，博士，博士生導師，E-mail:taldcn@yahoo.com.cn