MUSIC及其改進算法的研究與實現(xiàn)

黃麗薇，陳慧琴，陳玉林

(東南大學成賢學院，江蘇南京 210088)

陣列信號處理是信號處理的重要分支，信號到達角(DOA)估計是陣列信號的重要部分［1－2］。經(jīng)典MUSIC算法［3］基于接收信號的協(xié)方差矩陣分解，對于非相干信源的到達角估計具有良好的性能。前/后向空間平滑算法［4］和改進 MUSIC(MMUSIC)算法［5］可以實現(xiàn)相干信源的到達角估計。

1 經(jīng)典MUSIC算法

1.1 經(jīng)典MUSIC原理

考慮M元均勻線陣，有N個窄帶信源平面波入射，信源方向為 θ1，θ2，…，θN。S(k)=［s1(k)，…，sN(k)］T，si(k)為第i個信源的復振幅。陣列的導向矢量 a(θi)=［1，e－jωi，…，e－j(M－1)ωi］T，i=1，…，N，A=，d 為陣元間距，λ 為載波波長。n(k)=［n1(k)，…，nM(k)］T，ni(k)為零均值、方差σ2的白噪聲，與信源不相關(guān)。第k次快拍，得到的數(shù)據(jù)為X(k)=AS(k)+n(k)，k=1，2，…，K，K 為快拍次數(shù)。X(k)=［x1(k)，…，xM(k)］T為M個陣元輸出。由于信號與噪聲獨立，陣列數(shù)據(jù)X(k)的協(xié)方差矩陣Rx=E［X(k)XH(k)］=ARsAH+σ2I，Rs=E［SSH］為信號的相關(guān)矩陣，σ2I為噪聲的相關(guān)矩陣。

對Rs進行特征分解，得到N個特征值λ1，λ2，…，λN，定義對角陣∑S=diag(λ1，λ2，…，λN)，∑N=diag(λN+1，λN+2，…，λM)，其中噪聲相關(guān)矩陣的特征值λN+1= λN+2= … = λM= σ2。設(shè) λi為Rx的第i個特征值，vi是與 λi對應的特征向量，則有 Rxvi= λivi。設(shè)λi=σ2是最小特征值，則Rxvi= σ2vi(i=N+1，N+2，…，M)。將 Rx=ARsAH+ σ2I代入得，σ2vi=(ARsAH+ σ2I)vi。右式展開得 ARsAHvi=0，由于 A－1和 R－1存在，所以有 AHvi=0，i=N+1，N+2，…，M，即噪聲特征值對應的特征向量與信號源方向?qū)腁矩陣的列向量正交。用噪聲特征向量為列，構(gòu)造噪聲矩陣 En= ［vN+1，vN+2，…，vM］，定義空間譜 Pmusic(θ)=。當a(θ)和En的各列正交時，分母為0，由于噪聲存在，分母實際為一小值，對應Pmusic有尖峰。使θ改變，通過尋找尖峰來估計信源到達角。

1.2 非相干信號到達角估計MUSIC仿真

實例1 均勻直線陣陣元數(shù)為8，兩個獨立的窄帶信號源，載頻頻率分別為 π/3和 π/4，信噪比都為30 dB，從－30°和45°方向入射到直線陣，陣元間距取半波長，采樣數(shù)為1 024。

Matlab中關(guān)鍵實現(xiàn)程序:

x=A*s+N;%加了高斯白噪聲后的陣列接收信號

R=x*x'/n;

［V，D］=eig(R);%求R特征值和特征向量，D為特征值對角陣，V為特征向量陣

D=diag(D);

Un=V(:，1:M －P);%取所有行，1～6列

sdoa=－90:0.1:90;%搜索范圍為－90°～90°

for i=1:length(sdoa)

a_theta=exp(－j*(0:M－1)'*2* π*d*sin(π*sdoa(i)/180)/l);

Pmusic(i)=1./abs((a_theta)'*Un*Un'*a_theta);

end

圖1為非相干信號DOA估計對應仿真圖，空間譜峰值對應的角度為－30°和45°，經(jīng)典MUSIC算法可以準確估計信號到達角。

圖1 經(jīng)典MUSIC算法對非相干信號到達角估計

1.3 相干信號到達角估計MUSIC算法仿真

實例2 均勻直線陣陣元數(shù)為8，兩個相干窄帶信號源，載頻頻率均為π/6，信噪比為30 dB，從－30°和45°方向入射到直線陣，陣元間距取半波長，采樣數(shù)為1 024。

圖2為相干信號DOA估計對應仿真圖。可以看出，利用經(jīng)典MUSIC算法進行相干信源到達角估計，空間譜峰偏離實際值很多，估計失敗。因為MUSIC算法要求到達天線陣元的信號彼此不相關(guān)，因為只有這樣，信源的協(xié)方差矩陣才為滿秩。對于高度相關(guān)或相干的入射信號，由于接收數(shù)據(jù)的協(xié)方差矩陣的秩降為1，導致信號子空間的維數(shù)小于信號源數(shù)，即信號子空間擴散到了噪聲子空間，MUSIC算法的性能會急劇下降，無法正確估計信源方向。

圖2 經(jīng)典MUSIC算法對相干信號到達角估計

2 相干信號到達角估計的改進MUSIC算法

2.1 前向、后向、前后向空間平滑估計

2.1.1 空間平滑MUSIC算法原理

空間平滑MUSIC的思想是，改變經(jīng)典MUSIC算法中只有一個陣列的情況，將窄帶信源下的均勻線陣分成若個相互重疊的子陣列，各子陣列的協(xié)方差矩陣進行平均運算，實現(xiàn)去相關(guān)。前向空間平滑法把M個陣元分成p個子陣，每個子陣m個陣元，如圖3所示，取左側(cè)第一個子陣為參考子陣，對于第k個子陣的數(shù)據(jù)模型 xk(t)= ［xk，xk+1，…，xk+m－1］=AD(k－1)s(t)+nk(t)， 該子陣數(shù)據(jù)協(xié)方差矩陣 Rk=AD(k－1)Rs(D(k－1))HAH+ σ2I，前向空間平滑 MUSIC 方法對滿秩協(xié)方差矩陣的恢復，通過對各子陣協(xié)方差矩陣的均值實現(xiàn)，即前向平滑的協(xié)方差矩陣為 Rf=后向空間平滑法的子陣劃分如圖4所示，模型類似于前向平滑法，后向平滑的協(xié)方差矩陣為。前后向平滑法為兩個矩陣求平均，即

圖3 前向空間平滑算法子陣劃分圖

圖4 后向空間平滑算法子陣劃分圖

2.1.2 前后向空間平滑MUSIC算法仿真

用前后向空間平滑MUSIC算法實現(xiàn)實例二的DOA估計，Matlab中關(guān)鍵實現(xiàn)程序:

%前向平滑

xf1=x(［1:6］，:);Rf1=xf1*xf1'/n;

xf2=x(［2:7］，:);Rf2=xf2*xf2'/n;

xf3=x(［3:8］，:);Rf3=xf3*xf3'/n;

Rf=(Rf1+Rf2+Rf3)/3;

%后向平滑

xb1=conj(x(［8:－ 1:3］，:));Rb1=xb1*xb1'/n;

xb2=conj(x(［7:－ 1:2］，:));Rb2=xb2*xb2'/n;

xb3=conj(x(［6:－ 1:1］，:));Rb3=xb3*xb3'/n;

Rb=(Rb1+Rb2+Rb3)/3;

%前后向平滑

Rbf=(Rf+Rb)/2;

［U，S，V］=svd(Rbf);

Un=U(:，P+1:M);

圖5為前后向平滑MUSIC算法對相干信號DOA估計仿真圖，對相干信號的到達角－30°和45°能進行準確估計。采用空間平滑技術(shù)，預處理輸入信號的協(xié)方差，能有效估計相干信號的來波方向，精度較高。空間平滑的代價是，可估計信源數(shù)變?yōu)?M/3(M為陣元數(shù))，這是由于將接受陣列分成若干多個子陣，使陣面孔徑和陣元數(shù)減小。而由于孔徑變小，使得該算法對非相關(guān)信號源的到達角估計性能減小，一般只適合于均勻線陣。

圖5 前后向平滑MUSIC算法對相干信號到達角估計

2.2 改進MUSIC(MMUSIC)算法

2.2.1 改進MUSIC算法原理

在經(jīng)典MUSIC算法基礎(chǔ)上，在得到Rx=ARsAH+σ2I后，不直接對信號的相關(guān)矩陣Rs=E［SSH］進行特征分解，而令Y(k)=JMX*(k)，X*(k)為X(k)的復共軛，JM為M階交換矩陣，副對角線上元素為1，其余元素為0，可將向量倒排，圖6為交換矩陣示例。

圖6 改進MUSIC算法對相干信號到達角估計

可得Y(k)的相關(guān)矩陣 Ry=E［Y(k)YH(k)］=JMA*R*sATJM+σ2IM=JMR*xJM。令 R=Rx+Ry=Rx+JMR*xJM，對R特征分解，再用經(jīng)典MUSIC算法進行到達角估計［6］。本質(zhì)上，MMUSIC就是前后相平滑技術(shù)中，取子陣陣元數(shù)等于總陣元數(shù)的特殊情況。

2.2.2 相干信號到達角估計MMUSIC算法仿真

用MMUSIC算法實現(xiàn)實例二的DOA估計，Matlab中關(guān)鍵實現(xiàn)程序:

J=fliplr(eye(M));%產(chǎn)生交換矩陣

R=R+J*conj(R)*J;%conj求復數(shù)的共軛

［V，D］=eig(R)%;%D為對角陣，求R特征值和特征向量，V為8×8陣

D=diag(D);

Un=V(:，1:M －P);

圖6為改進MUSIC算法對相干信號DOA估計對應的仿真圖，對相干信號的到達角－30°和45°度能進行準確估計。改進MUSIC算法用修正后的R進行信號到達角估計，具有平均意義，可以提高估計性能，也不會影響對非相干信號的到達角估計。

3 結(jié)束語

對3種MUSIC算法進行了原理分析和實驗仿真，對各自的特點進行了闡述。算法優(yōu)缺點比較如表1所示，使用時可以根據(jù)實際情況和性能需要，加以選擇。

表1 3種算法優(yōu)缺點比較

［1］張賢達，保錚.通信信號處理［M］.北京:國防工業(yè)出版社，2000.

［2］張小飛，汪飛.陣列信號處理的理論與應用［M］.北京:國防工業(yè)出版社，2013.

［3］Schmidt R O.Multiple emitter location and signal parameter estimation［J］.IEEE Transactions on Antenna and Propagation，1986，34(3):276 －280.

［4］Pillai S U，Kwon B H.Forward/Backward spatial smoothing technique for coherent signal identification［J］.IEEE Transactions on ASSP，1989，37(1):8 －15.

［5］Kunda D.Modified music algorithm for estimating DOA of signals［J］.Signal Processing，1996，48(7):85 －89.

［6］何子述，黃振興，向敬成.修正MUSIC算法對相關(guān)信號源的 DOA 估計性能［J］.通信學報，2000，21(10):14－17.