空間譜估計經典算法性能比較

石國德，王明皓，呂朝暉

（1.沈陽航空航天大學 遼寧 沈陽 110136；2.沈陽飛機設計研究所 遼寧 沈陽 110035）

空間譜是陣列信號處理中的一個重要概念，“空間譜”表示信號在空間各個方向上的能量分布。若能夠估計出空間譜就可以得到信號源的波達方向 （direction of arrival，DOA），所以空間譜估計也常稱為DOA估計。

最早的陣列的DOA估計算法可以追述到常規波束形成算法（CBF）法[1]，也叫做Bartlett波束形成算法，該算法是時域傅里葉譜估計在空域的一種簡單擴展，即用空域各陣列數據代替傳統時域處理的時域數據，與時域傅里葉變換一樣空域分辨力受到空域傅里葉限即瑞利限的限制。故CBF不能分辨一個波束寬度內的多個信號源。接著突破瑞利限的DOA估計算法如非線性估計算法如Capon波束形成技術[2]得到了研究。

20世紀70年代起利用子空間分解類的算法開始興起，最具代表性的經典算法有：1）多重信號分類（MUSIC）[3]，即利用接收數據的協方差矩陣進行特征分解得到兩個相互正交的信號子空間和噪聲子空間，信號的子空間的劃分依據是與陣列流形空間構成相同的空間，剩余構成噪聲子空間。利用噪聲子空間與陣列流形正交構造出空間譜（并不代表能量分布，只說明噪聲子空間和陣列流形的正交程度），從而得到DOA估計。2）旋轉不變信號子空間算法（ESPRIT）[4-5]，是以陣列流形中隱含的空間平移不變性，即陣列形式可以劃分成兩個完全相同的子陣，從而利用與陣列形成相同空間的信號子空間得到旋轉不變方程，然后利用各種解旋轉不變方程方法如最小二乘 （LS）[4]、總體最小二乘（TLS）[5]解方程得到DOA估計。

20世紀80年代后期又出現了一批子空間擬合類算法，比較經典的有最大似然（ML）[6]算法、子空間擬合（SF）[7]算法，但由于DOA估計似然函數是非線性的，故求其最優解需要多維搜索，運算量比較大。

1 DOA估計經典算法

1.1 DOA估計信號模型

考慮p個遠場窄帶信號入射到空間M陣元的均勻線陣（ULA），這里假設陣元數等于通道數，即各陣元接收到信號后經各自的傳輸通道送入到處理器，還假設目標為點目標，第一個陣元為參考點，為避免測角模糊設陣元間距為半波長。則信號模型用矢量形式可表示為：

其中 X（t）表示 M×1 維接收快拍數據矢量，S（t）為 P×1 維信號源矢量，N（t）為M×1維噪聲數據矢量，為了方便分析這里假設噪聲為均值為0、方差為σ2的高斯白噪聲。A表示空間 M×N 的陣列流形矩陣，A=[a（θ1）,a（θ2），…，a（θP）]，其中第p（p=1，…，P）個信號陣列流形為 a（θp）=[1，exp（jπsinθp），…，exp（jπ（M-1）sinθp）]。 DOA 估計技術可理解為從模型（1）中估計出P個未知參量θp（p=1，…，P）。下面介紹6種經典算法。

1.2 常規波束形成（CBF）

假設 M 個陣元的加權矢量為 w=[w1，w2，…，wM]T，則整個陣列的輸出為y=wHx（t），那么L次快拍輸出的平均功率為

其中R為接收陣列的協方差矩陣。當加權值為w=a（θ）時，式（2）即為CBF空間譜圖，等于

通過掃描得到得到譜曲線，找到P個譜峰所對應的角度值就得到了DOA估計值。

1.3 Capon波束形成

可以證明式 （2）的波束形成的最優權矢量為wopt=則整個陣列的輸出可以簡化為

由于期待信號是未知的，可通過掃描得到Capon波束形成譜曲線即

亦叫做最小方差法（MVM），找到個譜峰所對應的角度值就得到了DOA估計值。

1.4 多重信號分類（MUSIC）

MUSIC算法的提出開創了信號譜估計算法研究的新時代，促進了特征結構類算法的興起和發展，該算法以成為空間譜估計的標志性算法。該算法不同于上述兩種算法，上述算法是針對接收數據協方差矩陣進行直接處理。而MUSIC算法需要對協方差矩陣進行特征值分解，即

US是由P個大特征值對應的特征矢量張成的信號子空間，UN是由M-P個小特征值對應的特征矢量張成的子空間也即噪聲子空間。需要指出的是這里的P值是未知的，需要用信息論方法或者平滑秩方法等事先估計。信號子空間和噪聲子空間是正交的，故理想條件下陣列流形正交于噪聲子空間，即aH（θ）UN=0。但是實際接收陣列協方差矩陣是有偏差的，可用有限個快拍數據來估計即行特征值分解得到近似噪聲子空間但該噪聲子空間并不與陣列流形完全正交。因此可以通過下面的空域搜索得到MUSIC空間譜：

1.5 旋轉不變子空間算法（ESPRIT）

ESPRIT算法最基本的假設是存在2個完全相同的子陣，且兩個子陣間的間距是Δ已知的。假設兩個相同的子陣分別為：

其中 φp=（2πΔsinθp/λ），由上述可知只要得到兩個子陣間的旋轉不變關系Φ就可求取DOA估計值。將兩個子陣模型合并得：

類似于MUSIC算法，對上述接收協方差矩陣進行特征分解得到信號子空間

則可得到

對式（12）進行解方程可得Ψ=U*S1US2（LS結果），上標*表示pense-moore偽逆。再對Ψ進行特征分解得到Φ，進而可求取DOA估計值。

1.6 最大似然估計（ML）

在對最大似然進行推導之前先限定幾點：1）信號協方差矩陣是正定的；2）陣元數大于信號源數，快拍數大于陣元數；3）不同快拍之間的噪聲是不相關的，且為服從正態分布的加性白噪聲。對于確定性對大似然即是假設信號源是確定未知的最大似然。L次快拍服從的聯合條件概率密度函數：

det{·}表示矩陣的行列式。對式（13）兩邊取負對數得到：

對于式（14）未知σ2，s并不是我們所關心的，故固定其余兩個未知量給出一個量的最大似然估計：=tr{P}，=A*x，其中 PA=A（AHA）-1AH表示矩陣 A 的投影矩陣，A*表示矩陣A的偽逆。把這兩個似然估計值代入式（14）并忽略常數項可得到θ的最大似然估計為：

1.7 子空間擬合（SF）

子空間擬合與最大似然類似，最大似然是接收數據與真實數據的擬合，而子空間擬合是信號子空間或噪聲子空間于真實的子空間的擬合。為了簡要說明這里只給出信號子空間擬合（SSF）。快拍數無限情況下有US=AT，但事實上快拍數有限情況下只是近似，可通過構造一個擬合關系來找出T使得兩者在最小二乘意義下擬合的最好，定義如下擬合關系：

T^作為輔助參量可以用最小二乘求解得：

將式（17）代入到式（16）得到SSF的求解的代價函數：

表1 各算法優缺點綜述Tab.1 Each algorithm and disadvantages summary

2 仿真比較與性能分析

2.1 單信號DOA估計性能隨信噪比變化曲線

考慮ULA的陣元數M=10，快拍數L=100，單個目標信號源DOA參數為θ=30°，圖1給出均方根誤差RMSE比較圖，RMSE定義如下

其中，k表示蒙特卡洛實驗次數，θ^，θ分別為信號源DOA估計值和真值。由圖1可以看出單信號源情況下，ML和CBF方法的均方根誤差最小，其次是MUSIC和SSF，接著為Capon，最差的是ESPRIT。但是文中的6種方法信噪比門限基本一致。

2.2 單信號DOA估計性能隨快拍數變化曲線

圖1 DOA估計性能隨信噪比變化曲線Fig.1 Estimation performance with SNR change curve

該仿真信噪比SNR=10 dB，快拍數變化，其余仿真條件與2.1相同，圖2給出單信號情況下DOA估計性能隨快拍數變化曲線。從圖2可以得到與2.1相同的結論，但需要注意的是在快拍數較小的時候（此處為10個快拍），Capon波束形成算法的性能是最差的。

圖2 DOA估計性能隨快拍數變化曲線Fig.2 Estimation performance with number of snapshots change curve

2.3 單信號DOA估計性能隨陣元數變化曲線

該仿真信噪比SNB=10 dB，其余仿真條件與2.1相同，圖3給出單信號情況下DOA估計性能隨陣元數變化曲線。從圖3得到與2.1相同的結論。

圖3 DOA估計性能隨陣元數變化曲線Fig.3 Estimation performance with number of array elements change curve

2.4 雙信號DOA估計性能隨信噪比變化曲線

考慮ULA的陣元數M=10，快拍數L=100，兩個目標信號源 DOA 參數為 θ=[30°，35°]， 圖 4給出信號比為 SNR=10 dB的空間譜圖，從圖4（a）中可以看出MUSIC算法分辨率高于Capon波束形成，而Capon波束形成高于CBF。

圖4 搜索算法譜圖Fig.4 Search algorithm spectrogram classification

3 結 論

文中對近幾十年的空間譜估計經典算法進行了總結，對比了各種經典算法的性能，列出了各種算法的優缺點，使得工程操作員在DOA估計算法實現時提供了理論依據，總結如下：1）若芯片處理速度不夠快，要求精度不高時，可選擇CBF和ESPRIT算法，這3種算法計算量相對來說不是很大，實時性較容易實現。2）芯片處理速度相對較快，要求精度較高時，可選擇Capon和MUSIC算法。3）芯片處理速度快，要求精度高時，可選擇ML和SF算法。

[1]Krim H，Viberg M.Two decades of array signal processing research[J].IEEE Singal Processing Magazine，1996，13（4）:67－94.

[2]Capon J.High-resolution entropy spectral analysis[J].Proceedings of IEEE，1969，57（8）:1408－1418.

[3]Schmidt R.Multiple emitter location and signal parameter estimation [J].IEEE Transactions on Antennas and Propagation，1986，34（3）:276－280.

[4]Roy R，Kailath T.ESPRIT-a subspace rotation approach to estimation of parameters of cissoids in noise [J].IEEE Transactions on Acoustics Speech and Signal Processing 1986，34（10）:1340－1342.

[5]Roy R，Kailath T.ESPRIT-estimation of signal parameters via rotational invariance techniques [J].IEEE Transactions on Acoustics Speech and Signal Processing，1989，37 （7）:984－995.

[6]Stoica P，Nehorai A.MUSIC，maximum likelihood，and cramerrao bound[J].IEEE Transactions on Acoustics Speech and Signal Processing，1989，37（5）:720－741.

[7]Viberg M，Ottersten B，Kailath T.Detection and estimation in sensor arrays using weighted subspace fitting[J].IEEE Transactions on Acoustics Speech and Signal Processing 1991，39（11）:2436－2449.