秦 陽，劉洪彪，胥勇軍，趙 樺，周陳超，房 敏，鄧朝平

（解放軍95806部隊，北京 100076）

0 引 言

在實際應(yīng)用中，常常需要對空間中存在的多個信號源進行分解，以便跟蹤或檢測感興趣的空間信號，抑制那些被認為是干擾的空間信號。例如對天線陣列接收的空間信號所進行的分析與處理稱為陣列信號處理。而空間譜估計技術(shù)是在波束形成技術(shù)、零點技術(shù)和時域譜估計技術(shù)的基礎(chǔ)上發(fā)展起來的一種技術(shù)。與頻譜表示信號在各個頻率上的能量分布相對應(yīng)，空間譜則可解釋為信號在空間各個方向上的能量分布，空間譜估計技術(shù)的目標是研究提高在處理帶寬內(nèi)空間信號角度的估計精度、角度分辨率和提高運算速度的各種算法。經(jīng)過多年的發(fā)展，已經(jīng)產(chǎn)生了大量性能優(yōu)異的測向算法可資利用，典型的有 MUSIC［1］、ESPRIT［2］等。

MUSIC算法是基于特征結(jié)構(gòu)分析的空間譜估計方法，是空間譜估計技術(shù)的典型代表。其測向原理是根據(jù)矩陣特征分解的理論［3］對陣列輸出協(xié)方差矩陣進行特征分解，將信號空間分解為噪聲子空間D和信號子空間C，利用噪聲子空間D與陣列的方向矩陣A的列矢量正交的性質(zhì)，構(gòu)造空間譜函數(shù)P（w）并進行譜峰搜索，從而估計出信號方向信息。本文通過仿真驗證了MUSIC算法可以很好地對信號參數(shù)進行估計。

1 MUSIC算法基本原理

1.1 信號子空間和噪聲子空間［4］

一個系統(tǒng)的輸入信號可以表示為（有用）信號S及噪聲N之和，即：

式中：X＝ ［x1，x2…xp］T，為輸入信號向量；S＝［s1，s2…sp］T， 為 有 用 信 號 向 量；N＝［n1，n2…np］T，為噪聲信號向量。

于是，可以定義輸入信號的相關(guān)函數(shù)矩陣Rxx為：

式中：“＊”表示復數(shù)共軛；由定義可知，信號矩陣Rxx為Hermite矩陣，即滿足Rxx＊＝RxxT。

因此，若輸入信號和噪聲信號互不相關(guān)，則相關(guān)矩陣Rxx可寫為：

式中：Rs和Rn分別為信號與噪聲的相關(guān)矩陣。

假如有p個互不相關(guān)的信號源同時作用于系統(tǒng)的m個輸入端，系統(tǒng)的輸入向量X構(gòu)成了p維酉向量空間Up，則第p個源在第m個輸入端上的輸入信號可以寫為：

式中：ap（t）為輸入信號的復包絡(luò)；ω為信號載頻；φpm為信號的相位滯后。

若第p個信號源的復包絡(luò)ap（t）為隨機的，則它的有用信號向量S也為一隨機向量?？梢杂盟趍個輸入端的相位滯后來定義一個信號特征向量Sp，即：

則由p個信號方向向量所張成的空間Span｛Sp，p＝1，2，…，P｝稱為在全空間UM中的p維信號子空間C（其中p≤m），而將它在全空間UM中的正交補空間C⊥稱為M-P維噪聲子空間D。

信號子空間和噪聲子空間的物理意義是明顯的，若一個理想的系統(tǒng)輸入無噪聲，則經(jīng)過統(tǒng)計平均所得的輸入向量必落在信號子空間C內(nèi)，但是若系統(tǒng)輸入有噪聲時，它就不能完全落在信號子空間內(nèi)，而它中間的噪聲分量則應(yīng)落在噪聲子空間C⊥內(nèi)。

1.2 MUSIC算法［5］

MUSIC算法的基本思想是將陣列輸出數(shù)據(jù)的協(xié)方差矩陣進行特征分解，從而得到與信號分量相對應(yīng)的信號子空間和與信號分量相正交的噪聲子空間，利用2個子空間的正交特性構(gòu)造出“針狀”空間譜峰來估計信號的方向信息。

設(shè)空間有p個互不相關(guān)的信號，以方位角θ1，θ2，…，θp入射到具有m個接收陣元的接收陣元陣列中，入射信號的數(shù)目p小于陣列的陣元數(shù)m，則此陣列系統(tǒng)的信號模型為：

對由上式描述的陣列信號觀測模型做以下假設(shè)：

假設(shè)1：對于不同的wi值，向量a（wi）相互線性獨立；

假設(shè)2：加性噪聲向量N（n）的每個原色都是零均值的復白噪聲，它們不相關(guān)，并且具有相同的方差σn2；

假設(shè)3：矩陣P＝E｛S（n）SH（n）｝非奇異，即rank（P）＝p。

上述3個假設(shè)都只是一般的假設(shè)，在實際中容易得到滿足?？梢姡?/p>

式中：AH為A的共軛轉(zhuǎn)置；I為單位矩陣。

于是，Rxx是個對稱矩陣，令其特征值分解為：

由于A列滿秩，故rank（APAH）＝rank（P）＝p，且這里假定p＜m，于是有：

式中：a12，…，ap2為無加性噪聲時的觀測信號AX（n）的自相關(guān)矩陣APAH的特征值。

同時用UH左乘和U右乘上式，得：

這表明，自相關(guān)矩陣Rxx的特征值為：

根據(jù)信號特征值和噪聲特征值，將特征矩陣U的列向量分為兩部分，即：

式中：信號子空間C＝ ［S1，…，SP］＝ ［U1，…，UP］，由信號特征向量組成；噪聲子空間D＝ ［g1，…，

又由Rxx＝APAH＋σn2I，有RxxD＝APAHD＋σn2D，利用上式的結(jié)果，得到：

上式成立的充分必要條件是AHD＝O。

將A＝ ［a（w1），…，a（wp）］代入式（15）得到：

顯然，當w≠w1，…，wp時，aH（w）D≠OT。

將上式改成標量形式，可以定義一種類似于功率譜的函數(shù)：gm-p］＝ ［Up＋1，…，Um］，由噪聲特征向量組成。

考察：

式（17）取峰值的p個w值給出p個信號的波達方向θ1，θ2，…，θp。

這樣定義的函數(shù)描述了空間參數(shù)（即波達方向）的分布，因此成為空間譜。它能夠?qū)Χ鄠€信號進行識別?？梢钥闯觯怂惴ǖ年P(guān)鍵是要對自相關(guān)矩陣進行特征分解。

2 計算機仿真研究

2.1 仿真數(shù)據(jù)模型

通過對前面MUSIC算法基本原理的分析，N個正弦信號加白噪聲的自相關(guān)函數(shù)可以表示為：

式中：rxx（n）為輸入信號的自相關(guān)函數(shù)；rss（n）為有用信號的自相關(guān)函數(shù)；ruu（n）為白噪聲信號的自相關(guān)函數(shù)；wi為正弦波的頻率；Ai為正弦波的幅度；σn2為白噪聲的方差；δ（k）為沖擊函數(shù)。

本文對這種正弦波加白噪聲的情形進行仿真，假定仿真的觀測數(shù)據(jù)分別由下面兩種情況產(chǎn)生：

（1）單個正弦信號檢測情況

式中：u（n）為一高斯白噪聲，其均值為0，方差為1，采用混合同余法產(chǎn)生。

用MUSIC方法估計觀測數(shù)據(jù)中正弦波的頻率，并給出白噪聲方差σn2與復正弦波的振幅A的估計值。

2.2 仿真圖形結(jié)果

westimate1＝0.250 0π（估計頻率），S＝1.405 8（噪聲估計功率），A＝3.983 5（信號振幅估計），仿真結(jié)果如圖1所示。

（2）多個正弦信號情況

圖1 單個正弦信號短數(shù)據(jù)（N＝64）時的估計結(jié)果

westimate1＝0.250 0π（估計頻率），S＝0.964 6（噪聲估計功率），A＝4.008 3（信號振幅估計），仿真結(jié)果如圖2所示。

圖2 單個正弦信號長數(shù)據(jù)（N＝1 024）時的估計結(jié)果

S＝2.1686（噪聲估計功率），A1＝3.887 2，A2＝2.187 4，A3＝2.910 4（信號振幅估計），仿真結(jié)果如圖3所示。

圖3 多個正弦信號短數(shù)據(jù)（N＝64）時的估計結(jié)果

S＝1.038 5（噪聲估計功率），A1＝4.009 8，A2＝2.021 8，A3＝2.991 7（信號振幅估計），仿真結(jié)果如圖4所示。

圖4 多個正弦信號長數(shù)據(jù)（N＝1 024）時的估計結(jié)果

2.3 仿真結(jié)果分析

（2）在信號數(shù)據(jù)較少，即信號長度較短的情況下（N＝64），噪聲功率和信號幅度的估計值有一定的偏差，這主要是因為信號短，加窗效應(yīng)明顯，自相關(guān)矩陣的估計誤差比較大，因此特征值分解存在誤差。增加信號長度可以減少信號自相關(guān)矩陣的估計誤差，能夠得到比較理想的結(jié)果，這可以從長數(shù)據(jù)（N＝1 024）中看出。

（3）在MUSIC方法中，前面幾個噪聲特征值對應(yīng)的特征矢量都是相同的，于是嘗試用單個特征矢量對單個正弦信號進行頻率估計，結(jié)果出現(xiàn)了虛假峰，即在不等于信號頻率的地方多出來一些峰值。研究特征矢量的正交性可知，特征矢量是兩兩正交的。對矩陣特征分解得到的特征矢量組成的矩陣U進行處理，程序如下：

計算得到的結(jié)果為：

圖5 出現(xiàn)虛假峰的結(jié)果

3 結(jié)束語

通過對MUSIC算法的仿真研究討論，可以看出該算法具有高分辨特性及較小的估計方差，有著優(yōu)異的性能。在信號處理中幾個重要問題，如本文討論的噪聲中幾個疊加信號的參數(shù)估計，還有窄帶傳感器陣列的方向搜索，以及重疊回波的分辨率等，都可以采取此算法進行處理。但采樣協(xié)方差矩陣的特征分解涉及大量不規(guī)則矩陣運算，運算量較大，限制了算法的實時運用。而且這里假定的信號源都是不相關(guān)的，對相關(guān)信號源運用MUSIC算法還需進行去相關(guān)處理［6-7］。因此，對算法的改進研究成為近年來許多學者研究的課題［8-10］。

［1］ Schmid T R O.Multiple emitter location and signal parameter estimation［J］.IEEE Transactions on AP，1986，34（3）：276-280.

［2］ ROY R，Kailath T.ESPRIT-subspace rotation approach to estimation of parameters of cissoids in noise［J］.IEEE Transactions on ASSP，1986，34（10）：1340-1343.

［3］ 程云鵬，張凱院，徐仲.矩陣論［M］.第3版.西安：西北工業(yè)大學出版社，2006.

［4］ 沈福民.自適應(yīng)信號處理［M］.西安：西安電子科技大學出版社，2001.

［5］ 王永良，陳輝，彭應(yīng)寧，等.空間譜估計理論與算法［M］.北京：清華大學出版社，2004.

［6］ Lin J D，F(xiàn)ang W H，Wang Y Y.FSF MUSIC for joint DOA and frequency estimation and its performance analysis［J］.IEEE Transactions on Signal Processing，2006，52（12）：4529-4542.

［7］ Mestr X，Langunas M A.Modified subspace algorithms for DOA estimation with large arrays［J］.IEEE Transactions on Signal Processing，2008，56（2）：598-614.

［8］ 黃磊，吳順君，張林讓，等.快速子空間分解方法及其維數(shù)的快速估計［J］.電子學報，2005，33（6）：977-981.

［9］ 張立杰，黃建國，史文濤，等.正交投影陣列信號子空間估計方法［J］.系統(tǒng)工程與電子技術(shù)學報，2009（9）：2063-2066.

［10］余繼周，陳定昌.一種DOA估計的快速子空間算法［J］.現(xiàn)代電子技術(shù)，2005，28（12）：90-92.