幾種空間譜估計(jì)測(cè)向算法的比較

朱肅嫻

（中國船舶重工集團(tuán)公司第七二三研究所，江蘇 揚(yáng)州 225001）

0 引 言

陣列信號(hào)處理的一個(gè)重要應(yīng)用是在空間譜估計(jì)測(cè)向中的應(yīng)用，而空間譜估計(jì)測(cè)向在雷達(dá)、通信等很多方面有著極其廣泛的應(yīng)用。介紹常用的空間譜估計(jì)算法，包括MUSIC算法、ESPRIT算法及改進(jìn)算法，并進(jìn)行仿真分析，比較其優(yōu)缺點(diǎn)，以期為裝備科研生產(chǎn)提供理論依據(jù)。

1 空間譜與DOA估計(jì)

對(duì)于均勻線陣，當(dāng)空間僅有一個(gè)信號(hào)源且不考慮信號(hào)噪聲時(shí)，均勻線陣接收到的信號(hào)可表示為：

定義xm(n)的空間傅里葉變換為：

因?yàn)椋?/p>

所以：

|x(φ)|2常稱為空間譜。當(dāng)φ=φ1時(shí)，式（4）取得最大值。因此，可以通過搜索|x(φ)|2峰值位置，利用φ1=2πdsinθ1/λ來確定信號(hào)源的方向。

2 空間譜估計(jì)測(cè)向

空間譜估計(jì)測(cè)向系統(tǒng)的原理[1]：空間譜估計(jì)測(cè)向采用信號(hào)來波在每個(gè)天線陣元上感應(yīng)產(chǎn)生電壓信號(hào)幅度和其信號(hào)相位，與來波方向相關(guān)對(duì)應(yīng)特性而實(shí)現(xiàn)對(duì)空間多個(gè)信號(hào)同時(shí)測(cè)向。空間譜估計(jì)測(cè)向應(yīng)用多元天線陣，如圖1所示。

圖1 空間譜估計(jì)測(cè)向系統(tǒng)組成

空間譜估計(jì)測(cè)向系統(tǒng)設(shè)計(jì)[2]：實(shí)現(xiàn)空間譜估計(jì)測(cè)向系統(tǒng)應(yīng)具備物理支持（天線陣列和數(shù)字接收機(jī)）和軟件系統(tǒng)支持。硬件和軟件兩者應(yīng)相輔相成，硬件的高能性、一致性可以減小采樣數(shù)據(jù)的誤差，可以很好地表現(xiàn)軟件的超分辨性能。

3 經(jīng)典空間譜估計(jì)算法

3.1 MUSIC算法

3.1.1 算法原理

假設(shè)空間遠(yuǎn)場(chǎng)信號(hào)源共有D個(gè)，各信號(hào)互不相關(guān)，各陣元噪聲nm(t),m=1,…,M也互不相關(guān)，噪聲與信號(hào)Sk(t),k=1,…,D也不相關(guān)[3]。

第m個(gè)陣元的輸出結(jié)果為：

式中：

θk為k個(gè)信號(hào)的到達(dá)方向（線陣的法向?yàn)?°），d為陣元間隔，λ為信號(hào)工作波長(zhǎng)。

將式（5）（共M個(gè)方程）寫成矩陣方程，有：

其中：

求各陣元輸出的空間相關(guān)矩陣，有：

式中：

其中，E{·}為求期望值，T為轉(zhuǎn)置，H為復(fù)共軛轉(zhuǎn)置，σ2為噪聲的方差。

對(duì)R進(jìn)行特征值分解，并將其特征值按降序排列，可以得出下面關(guān)系：

（1）特征值 λ1≥ λ2…λD＞ λD+1=…=λM；

（2）特征向量 V1、V2…VD、VD+1、VM。

進(jìn)一步分析，得如下結(jié)果：

（1）R的最小特征值等于σ2，重?cái)?shù)為M-D，即：

據(jù)此空間信號(hào)源的數(shù)目，可立即估計(jì)出來：

最小重?cái)?shù)為1。因此，M元陣可同時(shí)側(cè)向的信號(hào)數(shù)目的最大值為：

（2）矩陣R列空間的基是各相互正交的向量，ΩN為其中噪聲特征的子列空間項(xiàng)，Ωs為其中信號(hào)序列的子列空間，是由除最小特征值以外的特征值對(duì)應(yīng)的特征向量形成的子空間。在信號(hào)源所在方向上，諸方向向量a(θk),k=1,…D均處在信號(hào)子空間Ωs中。因此，構(gòu)造矩陣：

顯然，有：

MUSIC算法就是根據(jù)式（19）求空間譜PMU(θ)，有：

式中 ||·||2表示向量的 2范數(shù)，PMU(θ)峰值所對(duì)應(yīng)的θ值是信號(hào)源方向的估值[4]。

3.1.2 Matlab仿真

輸入信號(hào)：輸入三個(gè)不同頻率（即就是不相關(guān)信號(hào)）的復(fù)信號(hào)，頻率分別為1 100 Hz、1 300 Hz和1 500 Hz，陣元數(shù)目為8，采樣頻率為5 000 Hz，信源入射角度分別為20°、30°和45°。輸出結(jié)果分別如圖2、圖3、圖4所示。

圖2 信噪比SNR=0 dB，快拍數(shù)L=2 048的MUSIC算法仿真

圖3 信噪比SNR=20 dB，快拍數(shù)L=2 048的MUSIC算法仿真

圖4 信噪比SNR=20 dB，快拍數(shù)L=128的MUSIC算法仿真

對(duì)于圖2，信噪比SNR=0 dB，快拍數(shù)L=2 048，測(cè)向結(jié)果為20.007 9、30.149 2和45.160 7。

對(duì)于圖3，信噪比SNR=20 dB，快拍數(shù)L=2048，測(cè)向結(jié)果為20.065 2、30.034 7和45.046 2。

對(duì)于圖4，信噪比SNR=20 dB，快拍數(shù)L=128，測(cè)向結(jié)果為20.179 8、30.149 2和45.046 2。

通過對(duì)比可得，對(duì)于MUSIC算法，相對(duì)于快拍數(shù)的影響，信噪比的影響更明顯。由本算法的原理可知，該算法的缺點(diǎn)在于不能夠區(qū)分相干信號(hào)的信號(hào)源方向。

因此，為解決相干信號(hào)的處理問題，通常采用兩種算法：

（1）通過犧牲一定有效的陣元數(shù)量，使信號(hào)具有不相干性，即通過預(yù)處理使天線接收到的陣列信號(hào)失去相干性，然后應(yīng)用相關(guān)算法獲得精確的到達(dá)角。相關(guān)算法包括前后向預(yù)測(cè)投影矩陣方法、空間平滑方法和數(shù)據(jù)矩陣分解方法。

（2）應(yīng)用移動(dòng)陣列的方法或應(yīng)用頻率平滑的方法處理相干信號(hào)，但是不損失陣元數(shù)，應(yīng)用結(jié)合去相干與空間譜估計(jì)方法處理。此算法包含頻率平滑方法和旋轉(zhuǎn)不變子空間算法。

下面對(duì)常用的兩種算法（旋轉(zhuǎn)不變子空間算法和前向空間平滑算法）進(jìn)行應(yīng)用仿真，并比較其優(yōu)缺點(diǎn)。

3.2 前向空間平滑算法

3.2.1 前向平滑

眾所周知，均勻線陣（ULA）算法（如圖5所示）具有的特性是平移的固定性，所以將這個(gè)陣列采用相關(guān)特性可分化為與其互相重疊的L個(gè)分子陣列，然后它們相互關(guān)聯(lián)的每個(gè)子列陣中的陣元數(shù)量為m=M-L+1，且每個(gè)子列陣包含同一個(gè)信源方位的信息。下面按照多重信號(hào)分類算法計(jì)算這L個(gè)子列陣的自協(xié)方差矩陣，計(jì)算它們的算術(shù)平均值，然后形成一個(gè)相對(duì)等效的m階子列陣列的協(xié)方差矩陣Rf[5]。

圖5 平滑空間的計(jì)算方法的示意圖

對(duì)協(xié)方差矩陣方法求得的數(shù)值做特征值分解處理構(gòu)造噪聲子空間，根據(jù)正交性估計(jì)信源方向?？梢钥闯觯鬽＞D且L≥D時(shí)，經(jīng)過其子陣平滑的m階對(duì)應(yīng)的信源協(xié)方差矩陣的秩恢復(fù)為D。另外，空間平滑技術(shù)的缺點(diǎn)是由于子陣列比原陣列小導(dǎo)致陣列有效孔徑減小。

3.2.2 Matlab仿真

輸入信號(hào)：輸入為兩個(gè)頻率為1 000 Hz的相干信號(hào)，采樣頻率為2 000 Hz，陣元數(shù)目為8，快拍數(shù)為1 024，信源入射角度分別為-30°、45°，信噪比為20 dB。

輸出結(jié)果：分別輸出子陣個(gè)數(shù)為1、2、3、4、5、6、7、8，對(duì)應(yīng)子陣陣元數(shù)為8、7、6、5、4、3、2、1時(shí)，測(cè)向結(jié)果分別如圖6、圖7、圖8、圖9、圖10、圖11、圖12、圖13所示。

圖6 L=1，m=8的前向空間平滑算法仿真

圖7 L=2，m=7的前向空間平滑算法仿真

圖8 L=3，m=6的前向空間平滑算法仿真

圖9 L=4，m=5的前向空間平滑算法仿真

圖10 L=5，m=4的前向空間平滑算法仿真

圖11 L=6，m=3的前向空間平滑算法仿真

圖12 L=7，m=2的前向空間平滑算法仿真

圖13 L=8，m=1的前向空間平滑算法仿真

結(jié)論：當(dāng)信源數(shù)＜子陣陣元數(shù)且信源數(shù)≤信源數(shù)時(shí)，可以準(zhǔn)確測(cè)量經(jīng)過前后向平滑陣列協(xié)方差矩陣的信源數(shù)個(gè)相干信源。缺點(diǎn)在于，算法原理決定了該算法不能區(qū)分特別近情況下的相干信號(hào)[6]。

3.3 旋轉(zhuǎn)不變子空間算法

3.3.1 旋轉(zhuǎn)不變子空間算法原理

旋轉(zhuǎn)不變子空間算法包含MUSIC算法和ESPRIT算法，都是對(duì)陣列接收數(shù)據(jù)協(xié)方差矩陣進(jìn)行特征分解。雖然這兩種算法采用的特征分解特性不同，但是兩種算法可以互為補(bǔ)充進(jìn)行計(jì)算。

假設(shè)存在兩個(gè)一樣的子陣，且每個(gè)子列陣Δ確定，有著相同的信號(hào)，每個(gè)子陣列釋放的信號(hào)相位差Φ(i=1,2,…N)。

假設(shè)兩個(gè)子陣的第一個(gè)接收數(shù)據(jù)是X1，其中另外一子陣的收獲數(shù)據(jù)是X2，按照上述所述陣列的數(shù)學(xué)模型可以獲知：

式中，N1、N2為M-1維度的噪聲矢量代號(hào)；S為空間信號(hào)的N×1維度的信號(hào)矢量代號(hào)；A為陣列的M×N維度的流型矩陣。

其中：

根據(jù)Φ對(duì)角陣，根據(jù)公式可以簡(jiǎn)化矩陣如下：

可以發(fā)現(xiàn)，要能夠獲得到達(dá)角信號(hào)，只要得到每個(gè)子陣列關(guān)系的旋轉(zhuǎn)不變關(guān)系Φ。

下面只需從式（21）和式（22）中得到每個(gè)子陣列之間的關(guān)系。具體地，需要先將這兩個(gè)子陣列模型合為一，如下：

理想狀況下，計(jì)算出X協(xié)方差矩陣的R：

特征分解，得：

可以得到這些特征值之間的相互關(guān)系為：

再根據(jù)快速拍攝獲得的數(shù)據(jù)，對(duì)式（27）進(jìn)行修正得到結(jié)果：

于是，有：

根據(jù)只有一個(gè)非奇異矩陣T，得出：

得出的結(jié)構(gòu)對(duì)每個(gè)子陣列都是可以應(yīng)用的，如下：

可以發(fā)現(xiàn)，US1=US2=A，即有：

再根據(jù)陣列流型在每個(gè)子陣列上如下的關(guān)系：

采用式（31）計(jì)算這兩個(gè)子陣間對(duì)應(yīng)的信號(hào)子空間有如下關(guān)系：

若US1滿秩，則具有唯一的解釋，可以得出：

然后特征分解：

式中ψ的特征值對(duì)應(yīng)構(gòu)成其對(duì)角陣在數(shù)值上和Φ一定等同，且ψ特征矢量是矩陣的各列。因此，可以獲得旋轉(zhuǎn)不變關(guān)系矩陣ψ，然后可以應(yīng)用公式得信號(hào)的來源。

3.3.2 Matlab仿真

輸入信號(hào)：輸入為兩個(gè)頻率為1 000 Hz的相干信號(hào)，采樣頻率為2 000 Hz，陣元數(shù)目為8，快拍數(shù)為1 024，子陣數(shù)為4，信源入射角度分別為-45°、30°、60°，信噪比為20 dB。

輸出結(jié)果，如表1所示。

表1 旋轉(zhuǎn)不變子空間算法仿真

4 結(jié) 論

空間譜測(cè)向中兩種典型的算法——MUSIC算法和旋轉(zhuǎn)不變子空間（ESPRIT）算法。對(duì)陣列接收數(shù)據(jù)協(xié)方差矩陣進(jìn)行特征值分解，是這兩種算法的相同點(diǎn)。MUSIC算法和旋轉(zhuǎn)不變子空間（ESPRIT）算法差異是采用特制分解的特性不同，且它們之間不是相互獨(dú)立的。MUSIC算法的改進(jìn)方法是使用前向空間平滑算法，使其可以適用于相干信號(hào)源。以上算法通過仿真驗(yàn)證可以看出，采樣率、快拍數(shù)、信號(hào)的相干性都對(duì)測(cè)向有一定影響。

通過仿真分析可以看出，MUSIC算法和ESPRIT算法主要采用了接收數(shù)據(jù)協(xié)方差矩陣的噪聲子空間的正交特性和旋轉(zhuǎn)固定性（計(jì)算量小，不需要進(jìn)行譜峰搜索）。實(shí)際工程中，可根據(jù)對(duì)應(yīng)的信號(hào)相關(guān)性、測(cè)向分辨率的要求及硬件約束來選擇恰當(dāng)?shù)淖V估計(jì)算法。