幾種空間譜估計測向算法的比較

朱肅嫻

（中國船舶重工集團公司第七二三研究所，江蘇 揚州 225001）

0 引 言

陣列信號處理的一個重要應用是在空間譜估計測向中的應用，而空間譜估計測向在雷達、通信等很多方面有著極其廣泛的應用。介紹常用的空間譜估計算法，包括MUSIC算法、ESPRIT算法及改進算法，并進行仿真分析，比較其優缺點，以期為裝備科研生產提供理論依據。

1 空間譜與DOA估計

對于均勻線陣，當空間僅有一個信號源且不考慮信號噪聲時，均勻線陣接收到的信號可表示為：

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

因為：

所以：

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

2 空間譜估計測向

空間譜估計測向系統的原理[1]：空間譜估計測向采用信號來波在每個天線陣元上感應產生電壓信號幅度和其信號相位，與來波方向相關對應特性而實現對空間多個信號同時測向。空間譜估計測向應用多元天線陣，如圖1所示。

圖1 空間譜估計測向系統組成

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

3 經典空間譜估計算法

3.1 MUSIC算法

3.1.1 算法原理

假設空間遠場信號源共有D個，各信號互不相關，各陣元噪聲nm(t),m=1,…,M也互不相關，噪聲與信號Sk(t),k=1,…,D也不相關[3]。

第m個陣元的輸出結果為：

式中：

θk為k個信號的到達方向（線陣的法向為0°），d為陣元間隔，λ為信號工作波長。

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

其中：

求各陣元輸出的空間相關矩陣，有：

式中：

其中，E{·}為求期望值，T為轉置，H為復共軛轉置，σ2為噪聲的方差。

對R進行特征值分解，并將其特征值按降序排列，可以得出下面關系：

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

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

進一步分析，得如下結果：

（1）R的最小特征值等于σ2，重數為M-D，即：

據此空間信號源的數目，可立即估計出來：

最小重數為1。因此，M元陣可同時側向的信號數目的最大值為：

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

顯然，有：

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

式中 ||·||2表示向量的 2范數，PMU(θ)峰值所對應的θ值是信號源方向的估值[4]。

3.1.2 Matlab仿真

輸入信號：輸入三個不同頻率（即就是不相關信號）的復信號，頻率分別為1 100 Hz、1 300 Hz和1 500 Hz，陣元數目為8，采樣頻率為5 000 Hz，信源入射角度分別為20°、30°和45°。輸出結果分別如圖2、圖3、圖4所示。

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

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

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

對于圖2，信噪比SNR=0 dB，快拍數L=2 048，測向結果為20.007 9、30.149 2和45.160 7。

對于圖3，信噪比SNR=20 dB，快拍數L=2048，測向結果為20.065 2、30.034 7和45.046 2。

對于圖4，信噪比SNR=20 dB，快拍數L=128，測向結果為20.179 8、30.149 2和45.046 2。

通過對比可得，對于MUSIC算法，相對于快拍數的影響，信噪比的影響更明顯。由本算法的原理可知，該算法的缺點在于不能夠區分相干信號的信號源方向。

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

（1）通過犧牲一定有效的陣元數量，使信號具有不相干性，即通過預處理使天線接收到的陣列信號失去相干性，然后應用相關算法獲得精確的到達角。相關算法包括前后向預測投影矩陣方法、空間平滑方法和數據矩陣分解方法。

（2）應用移動陣列的方法或應用頻率平滑的方法處理相干信號，但是不損失陣元數，應用結合去相干與空間譜估計方法處理。此算法包含頻率平滑方法和旋轉不變子空間算法。

下面對常用的兩種算法（旋轉不變子空間算法和前向空間平滑算法）進行應用仿真，并比較其優缺點。

3.2 前向空間平滑算法

3.2.1 前向平滑

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

圖5 平滑空間的計算方法的示意圖

對協方差矩陣方法求得的數值做特征值分解處理構造噪聲子空間，根據正交性估計信源方向。可以看出，若m＞D且L≥D時，經過其子陣平滑的m階對應的信源協方差矩陣的秩恢復為D。另外，空間平滑技術的缺點是由于子陣列比原陣列小導致陣列有效孔徑減小。

3.2.2 Matlab仿真

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

輸出結果：分別輸出子陣個數為1、2、3、4、5、6、7、8，對應子陣陣元數為8、7、6、5、4、3、2、1時，測向結果分別如圖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的前向空間平滑算法仿真

結論：當信源數＜子陣陣元數且信源數≤信源數時，可以準確測量經過前后向平滑陣列協方差矩陣的信源數個相干信源。缺點在于，算法原理決定了該算法不能區分特別近情況下的相干信號[6]。

3.3 旋轉不變子空間算法

3.3.1 旋轉不變子空間算法原理

旋轉不變子空間算法包含MUSIC算法和ESPRIT算法，都是對陣列接收數據協方差矩陣進行特征分解。雖然這兩種算法采用的特征分解特性不同，但是兩種算法可以互為補充進行計算。

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

假設兩個子陣的第一個接收數據是X1，其中另外一子陣的收獲數據是X2，按照上述所述陣列的數學模型可以獲知：

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

其中：

根據Φ對角陣，根據公式可以簡化矩陣如下：

可以發現，要能夠獲得到達角信號，只要得到每個子陣列關系的旋轉不變關系Φ。

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

理想狀況下，計算出X協方差矩陣的R：

特征分解，得：

可以得到這些特征值之間的相互關系為：

再根據快速拍攝獲得的數據，對式（27）進行修正得到結果：

于是，有：

根據只有一個非奇異矩陣T，得出：

得出的結構對每個子陣列都是可以應用的，如下：

可以發現，US1=US2=A，即有：

再根據陣列流型在每個子陣列上如下的關系：

采用式（31）計算這兩個子陣間對應的信號子空間有如下關系：

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

然后特征分解：

式中ψ的特征值對應構成其對角陣在數值上和Φ一定等同，且ψ特征矢量是矩陣的各列。因此，可以獲得旋轉不變關系矩陣ψ，然后可以應用公式得信號的來源。

3.3.2 Matlab仿真

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

輸出結果，如表1所示。

表1 旋轉不變子空間算法仿真

4 結 論

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

通過仿真分析可以看出，MUSIC算法和ESPRIT算法主要采用了接收數據協方差矩陣的噪聲子空間的正交特性和旋轉固定性（計算量小，不需要進行譜峰搜索）。實際工程中，可根據對應的信號相關性、測向分辨率的要求及硬件約束來選擇恰當的譜估計算法。