黃麗薇，陳慧琴，陳玉林

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

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

1 經(jīng)典MUSIC算法

1.1 經(jīng)典MUSIC原理

考慮M元均勻線陣，有N個(gè)窄帶信源平面波入射，信源方向?yàn)?θ1，θ2，…，θN。S(k)=［s1(k)，…，sN(k)］T，si(k)為第i個(gè)信源的復(fù)振幅。陣列的導(dǎo)向矢量 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個(gè)陣元輸出。由于信號(hào)與噪聲獨(dú)立，陣列數(shù)據(jù)X(k)的協(xié)方差矩陣Rx=E［X(k)XH(k)］=ARsAH+σ2I，Rs=E［SSH］為信號(hào)的相關(guān)矩陣，σ2I為噪聲的相關(guān)矩陣。

對Rs進(jìn)行特征分解，得到N個(gè)特征值λ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個(gè)特征值，vi是與 λi對應(yīng)的特征向量，則有 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，即噪聲特征值對應(yīng)的特征向量與信號(hào)源方向?qū)?yīng)的A矩陣的列向量正交。用噪聲特征向量為列，構(gòu)造噪聲矩陣 En= ［vN+1，vN+2，…，vM］，定義空間譜 Pmusic(θ)=。當(dāng)a(θ)和En的各列正交時(shí)，分母為0，由于噪聲存在，分母實(shí)際為一小值，對應(yīng)Pmusic有尖峰。使θ改變，通過尋找尖峰來估計(jì)信源到達(dá)角。

1.2 非相干信號(hào)到達(dá)角估計(jì)MUSIC仿真

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

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

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

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為非相干信號(hào)DOA估計(jì)對應(yīng)仿真圖，空間譜峰值對應(yīng)的角度為－30°和45°，經(jīng)典MUSIC算法可以準(zhǔn)確估計(jì)信號(hào)到達(dá)角。

圖1 經(jīng)典MUSIC算法對非相干信號(hào)到達(dá)角估計(jì)

1.3 相干信號(hào)到達(dá)角估計(jì)MUSIC算法仿真

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

圖2為相干信號(hào)DOA估計(jì)對應(yīng)仿真圖?？梢钥闯觯媒?jīng)典MUSIC算法進(jìn)行相干信源到達(dá)角估計(jì)，空間譜峰偏離實(shí)際值很多，估計(jì)失敗。因?yàn)镸USIC算法要求到達(dá)天線陣元的信號(hào)彼此不相關(guān)，因?yàn)橹挥羞@樣，信源的協(xié)方差矩陣才為滿秩。對于高度相關(guān)或相干的入射信號(hào)，由于接收數(shù)據(jù)的協(xié)方差矩陣的秩降為1，導(dǎo)致信號(hào)子空間的維數(shù)小于信號(hào)源數(shù)，即信號(hào)子空間擴(kuò)散到了噪聲子空間，MUSIC算法的性能會(huì)急劇下降，無法正確估計(jì)信源方向。

圖2 經(jīng)典MUSIC算法對相干信號(hào)到達(dá)角估計(jì)

2 相干信號(hào)到達(dá)角估計(jì)的改進(jìn)MUSIC算法

2.1 前向、后向、前后向空間平滑估計(jì)

2.1.1 空間平滑MUSIC算法原理

空間平滑MUSIC的思想是，改變經(jīng)典MUSIC算法中只有一個(gè)陣列的情況，將窄帶信源下的均勻線陣分成若個(gè)相互重疊的子陣列，各子陣列的協(xié)方差矩陣進(jìn)行平均運(yùn)算，實(shí)現(xiàn)去相關(guān)。前向空間平滑法把M個(gè)陣元分成p個(gè)子陣，每個(gè)子陣m個(gè)陣元，如圖3所示，取左側(cè)第一個(gè)子陣為參考子陣，對于第k個(gè)子陣的數(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é)方差矩陣的恢復(fù)，通過對各子陣協(xié)方差矩陣的均值實(shí)現(xiàn)，即前向平滑的協(xié)方差矩陣為 Rf=后向空間平滑法的子陣劃分如圖4所示，模型類似于前向平滑法，后向平滑的協(xié)方差矩陣為。前后向平滑法為兩個(gè)矩陣求平均，即

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

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

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

用前后向空間平滑MUSIC算法實(shí)現(xiàn)實(shí)例二的DOA估計(jì)，Matlab中關(guān)鍵實(shí)現(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算法對相干信號(hào)DOA估計(jì)仿真圖，對相干信號(hào)的到達(dá)角－30°和45°能進(jìn)行準(zhǔn)確估計(jì)。采用空間平滑技術(shù)，預(yù)處理輸入信號(hào)的協(xié)方差，能有效估計(jì)相干信號(hào)的來波方向，精度較高?？臻g平滑的代價(jià)是，可估計(jì)信源數(shù)變?yōu)?M/3(M為陣元數(shù))，這是由于將接受陣列分成若干多個(gè)子陣，使陣面孔徑和陣元數(shù)減小。而由于孔徑變小，使得該算法對非相關(guān)信號(hào)源的到達(dá)角估計(jì)性能減小，一般只適合于均勻線陣。

圖5 前后向平滑MUSIC算法對相干信號(hào)到達(dá)角估計(jì)

2.2 改進(jìn)MUSIC(MMUSIC)算法

2.2.1 改進(jìn)MUSIC算法原理

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

圖6 改進(jìn)MUSIC算法對相干信號(hào)到達(dá)角估計(jì)

可得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算法進(jìn)行到達(dá)角估計(jì)［6］。本質(zhì)上，MMUSIC就是前后相平滑技術(shù)中，取子陣陣元數(shù)等于總陣元數(shù)的特殊情況。

2.2.2 相干信號(hào)到達(dá)角估計(jì)MMUSIC算法仿真

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

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

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

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

D=diag(D);

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

圖6為改進(jìn)MUSIC算法對相干信號(hào)DOA估計(jì)對應(yīng)的仿真圖，對相干信號(hào)的到達(dá)角－30°和45°度能進(jìn)行準(zhǔn)確估計(jì)。改進(jìn)MUSIC算法用修正后的R進(jìn)行信號(hào)到達(dá)角估計(jì)，具有平均意義，可以提高估計(jì)性能，也不會(huì)影響對非相干信號(hào)的到達(dá)角估計(jì)。

3 結(jié)束語

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

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

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

［2］張小飛，汪飛.陣列信號(hào)處理的理論與應(yīng)用［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)信號(hào)源的 DOA 估計(jì)性能［J］.通信學(xué)報(bào)，2000，21(10):14－17.