石國(guó)德，王明皓，呂朝暉

（1.沈陽航空航天大學(xué) 遼寧 沈陽 110136；2.沈陽飛機(jī)設(shè)計(jì)研究所 遼寧 沈陽 110035）

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

最早的陣列的DOA估計(jì)算法可以追述到常規(guī)波束形成算法（CBF）法[1]，也叫做Bartlett波束形成算法，該算法是時(shí)域傅里葉譜估計(jì)在空域的一種簡(jiǎn)單擴(kuò)展，即用空域各陣列數(shù)據(jù)代替?zhèn)鹘y(tǒng)時(shí)域處理的時(shí)域數(shù)據(jù)，與時(shí)域傅里葉變換一樣空域分辨力受到空域傅里葉限即瑞利限的限制。故CBF不能分辨一個(gè)波束寬度內(nèi)的多個(gè)信號(hào)源。接著突破瑞利限的DOA估計(jì)算法如非線性估計(jì)算法如Capon波束形成技術(shù)[2]得到了研究。

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

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

1 DOA估計(jì)經(jīng)典算法

1.1 DOA估計(jì)信號(hào)模型

考慮p個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)入射到空間M陣元的均勻線陣（ULA），這里假設(shè)陣元數(shù)等于通道數(shù)，即各陣元接收到信號(hào)后經(jīng)各自的傳輸通道送入到處理器，還假設(shè)目標(biāo)為點(diǎn)目標(biāo)，第一個(gè)陣元為參考點(diǎn)，為避免測(cè)角模糊設(shè)陣元間距為半波長(zhǎng)。則信號(hào)模型用矢量形式可表示為：

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

1.2 常規(guī)波束形成（CBF）

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

其中R為接收陣列的協(xié)方差矩陣。當(dāng)加權(quán)值為w=a（θ）時(shí)，式（2）即為CBF空間譜圖，等于

通過掃描得到得到譜曲線，找到P個(gè)譜峰所對(duì)應(yīng)的角度值就得到了DOA估計(jì)值。

1.3 Capon波束形成

可以證明式 （2）的波束形成的最優(yōu)權(quán)矢量為wopt=則整個(gè)陣列的輸出可以簡(jiǎn)化為

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

亦叫做最小方差法（MVM），找到個(gè)譜峰所對(duì)應(yīng)的角度值就得到了DOA估計(jì)值。

1.4 多重信號(hào)分類（MUSIC）

MUSIC算法的提出開創(chuàng)了信號(hào)譜估計(jì)算法研究的新時(shí)代，促進(jìn)了特征結(jié)構(gòu)類算法的興起和發(fā)展，該算法以成為空間譜估計(jì)的標(biāo)志性算法。該算法不同于上述兩種算法，上述算法是針對(duì)接收數(shù)據(jù)協(xié)方差矩陣進(jìn)行直接處理。而MUSIC算法需要對(duì)協(xié)方差矩陣進(jìn)行特征值分解，即

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

1.5 旋轉(zhuǎn)不變子空間算法（ESPRIT）

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

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

類似于MUSIC算法，對(duì)上述接收協(xié)方差矩陣進(jìn)行特征分解得到信號(hào)子空間

則可得到

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

1.6 最大似然估計(jì)（ML）

在對(duì)最大似然進(jìn)行推導(dǎo)之前先限定幾點(diǎn)：1）信號(hào)協(xié)方差矩陣是正定的；2）陣元數(shù)大于信號(hào)源數(shù)，快拍數(shù)大于陣元數(shù)；3）不同快拍之間的噪聲是不相關(guān)的，且為服從正態(tài)分布的加性白噪聲。對(duì)于確定性對(duì)大似然即是假設(shè)信號(hào)源是確定未知的最大似然。L次快拍服從的聯(lián)合條件概率密度函數(shù)：

det{·}表示矩陣的行列式。對(duì)式（13）兩邊取負(fù)對(duì)數(shù)得到：

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

1.7 子空間擬合（SF）

子空間擬合與最大似然類似，最大似然是接收數(shù)據(jù)與真實(shí)數(shù)據(jù)的擬合，而子空間擬合是信號(hào)子空間或噪聲子空間于真實(shí)的子空間的擬合。為了簡(jiǎn)要說明這里只給出信號(hào)子空間擬合（SSF）?？炫臄?shù)無限情況下有US=AT，但事實(shí)上快拍數(shù)有限情況下只是近似，可通過構(gòu)造一個(gè)擬合關(guān)系來找出T使得兩者在最小二乘意義下擬合的最好，定義如下擬合關(guān)系：

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

將式（17）代入到式（16）得到SSF的求解的代價(jià)函數(shù)：

表1 各算法優(yōu)缺點(diǎn)綜述Tab.1 Each algorithm and disadvantages summary

2 仿真比較與性能分析

2.1 單信號(hào)DOA估計(jì)性能隨信噪比變化曲線

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

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

2.2 單信號(hào)DOA估計(jì)性能隨快拍數(shù)變化曲線

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

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

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

2.3 單信號(hào)DOA估計(jì)性能隨陣元數(shù)變化曲線

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

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

2.4 雙信號(hào)DOA估計(jì)性能隨信噪比變化曲線

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

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

3 結(jié) 論

文中對(duì)近幾十年的空間譜估計(jì)經(jīng)典算法進(jìn)行了總結(jié)，對(duì)比了各種經(jīng)典算法的性能，列出了各種算法的優(yōu)缺點(diǎn)，使得工程操作員在DOA估計(jì)算法實(shí)現(xiàn)時(shí)提供了理論依據(jù)，總結(jié)如下：1）若芯片處理速度不夠快，要求精度不高時(shí)，可選擇CBF和ESPRIT算法，這3種算法計(jì)算量相對(duì)來說不是很大，實(shí)時(shí)性較容易實(shí)現(xiàn)。2）芯片處理速度相對(duì)較快，要求精度較高時(shí)，可選擇Capon和MUSIC算法。3）芯片處理速度快，要求精度高時(shí)，可選擇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.