馬 騰，杜 江*，喬理揚，黃天賜

(1.成都信息工程大學 通信工程學院，四川 成都 610225;2.氣象信息與信號處理四川省高校重點實驗室，四川 成都 610225;3.四川大學 電氣工程學院，四川 成都 610065)

0 引言

波達方向(Direction of Arrival,DOA)估計在雷達及通信等領域有著廣泛的應用且是陣列信號處理的重要研究分支及組成部分[1-2]。經典的基于子空間的多重信號分類算法(MUSIC)[3]開創(chuàng)了超分辨空間譜參數估計的先河，該算法測向精度高且理論推導不復雜，但運算量較大，在有些場合不太適合信號的實時處理。其運算量主要集中在陣列輸出的協(xié)方差矩陣的計算、協(xié)方差特征分解得到信號和噪聲子空間及譜峰搜索，且譜峰搜索隨著搜索精度的提升，運算量急劇增加。為了避免或減少譜峰搜索從而降低運算量，有學者提出了經典的求根MUSIC(MUSIC-root)[3]、旋轉不變子空間算法(ESPRIT)[4]，以及閆鋒剛等[5-8]提出的一整套壓縮譜理論，但仍需計算陣列輸出的協(xié)方差及對其特征分解。此外，為了減少整體算法中復矩陣間的運算量，酉變換被引入并衍生了一系列算法[9-10]，其使復矩陣間運算變?yōu)閷嵕仃囬g運算，大大降低了計算復雜度。

為了避免直接對所有陣列輸出計算協(xié)方差及對其特征分解得到信號或噪聲子空間，降低運算復雜度，一些快速參數估計算法被提出。如傳播算子算法(PM)[4]及其衍生算法[11],通過線性運算代替協(xié)方差矩陣的特征分解，間接得出信號及噪聲子空間。應用多級維納濾波(MSWF)[11-12]技術進行快速DOA估計,通過多級維納濾波前向遞推間接獲得信號子空間和噪聲子空間，避免了直接對陣列輸出的協(xié)方差矩陣的特征分解。但基于多級維納濾波技術的方法需要知道原始信號作為參考，實際應用性不強。近些年，一種Nystr?m方法被用于DOA估計中[13-15]。其最初應用于機器學習中的核機器學習領域[16]。該方法在離散空間中也可以用來逼近對應數據的特征矢量，即Nystr?m方法的矩陣表示形式。目前Nystr?m方法是一種較流行的“低秩矩陣逼近算法”之一。在DOA估計中應用該方法，無需計算所有陣列輸出數據的協(xié)方差矩陣及對其直接進行特征值分解得到信號子空間，而是通過陣列輸出的部分塊矩陣直接得到一種等效的逼近信號子空間，極大地減少了計算復雜度。

本文將結合Nystr?m方法提出一種新穎的DOA估計算法，在得到逼近的信號子空間后，無需像大多數現有的DOA估計算法通過矩陣間的運算或譜峰搜索得到DOA估計，而是通過逼近的信號子空間構建低階的特征多項式方程，此多項式方程的階數僅與信號源數目相同(信號源數通常遠少于陣列陣元數)，不同于現有的DOA估計求根類算法所建立的多項式，其階數大約為2倍的陣列陣元數[11，17-18]，對此低階多項式求根相比高階多項式方程求根得到DOA估計更加提升了計算效率。通過理論分析及仿真實驗對比其他算法說明了所提算法的性能。

1 陣列信號數學模型

天線陣列模型如圖1所示，假設有K個遠場非相干窄帶信號入射到共有M個天線陣元的一均勻線性陣列(ULA)，各陣元之間間距d為傳播波長λ的一半。

圖1 天線陣列模型

陣列在t時刻對應第m個陣元的輸出可表示為：

(1)

式中，sk(t)為第k個入射信號；nm(t)為第m個陣元受到的高斯白噪聲干擾；θk為第k個信號的DOA。定義θk方向向量為：

(2)

陣列流型矩陣A定義為A=[a(θ1),a(θ2),…,a(θk)]，所有M個陣元的輸出為x(t)=[x1(t),x2(t),…,xM(t)]T。當快拍數為L時，陣列輸出可寫成如下形式：

X(t)=[x(t1),x(t2),…,x(tL)]=AS(t)+N(t)，

(3)

式中，X(t)為M×L維陣列輸出數據；S(t)=[s(t1),s(t2),…,s(tL)]為K×L維入射到陣列的信號；N(t)=[n(t1),n(t2),…,n(tL)]為M×L維加性高斯白噪聲。則可定義陣列輸出的協(xié)方差矩陣為：

R=E[X(t)XH(t)]=ARSAH+σ2IM，

(4)

式中，RS=E[S(t)SH(t)]為信號源S(t)的信號協(xié)方差矩陣；IM為M×M維單位矩陣。通常在實際應用中，式(4)的協(xié)方差矩陣R通過下式估計得到：

(5)

在得到陣列輸出的協(xié)方差矩陣R后，其特征分解可表示為：

(6)

式中，US為信號子空間，由與信號源數目相同的K個較大特征值所對應的特征向量構成；UN為其余的M-K個特征向量組成的噪聲子空間。得到信號及噪聲子空間后，可通過MUSIC,ESPRIT等算法及其衍生算法進行DOA估計。下一節(jié)將對提出的新穎的DOA估計算法進行闡述。

本文中右上標“*”表示復共軛，“T”表示轉置，“H”表示共軛轉置，“?”表示矩陣的求逆運算。

2 所提算法

2.1 通過Nystr?m方法得到逼近的信號子空間

首先通過Nystr?m方法得到等效逼近的信號子空間，Nystr?m方法的本質是一種快速的“低秩矩陣逼近”算法之一，只用計算部分陣列輸出數據所構成的協(xié)方差矩陣及對其特征分解就可間接獲得逼近的信號子空間，避免了直接計算所有陣列輸出數據所構成的協(xié)方差矩陣及對其直接進行特征分解，從而減少了運算量。由文獻[13-15]可知，逼近的信號子空間可通過如下的Nystr?m方法得到。

將陣列輸出X(t)=[x(t1),x(t2),…,x(tL)]分割成2個子矩陣：

(7)

式中，X1為由X的前q行組成的q×L維矩陣；X2為由X的后M-q行組成的(M-q)×L維矩陣。這里的q為一個定義的參數，需滿足q∈{1,2,…,M}。2個子矩陣各自的協(xié)方差矩陣為：

(8)

(9)

式中，A1為由陣列流型矩陣A的前q行組成的矩陣；A2為其后M-q行組成的矩陣；RS為信號源的協(xié)方差矩陣。為了確保R11必須是一個滿秩矩陣，則q應該大于信源數K且小于min(M,L)。通常，q的值無需隨著陣元數M的增加而增加，當陣元數增加到較大時，依然可以選擇一個較小的q值來分割矩陣，足以確保較準確的DOA估計以及同時減少計算量。在實際應用中，式(8)與式(9)的協(xié)方差矩陣R11，R21由下式估計得到：

(10)

(11)

R21U11=U21Λ11，

(12)

(13)

(14)

很容易證明上式特征向量矩陣U滿足列矢量的相互正交特性，即UHU=I。且由協(xié)方差矩陣R特征分解得到的信號子空間US和近似特征向量矩陣U的前K列組成的矩陣張成同一空間，即span{US}=span{U(:,1:K)}[13]。

2.2 構建低階特征方程進行DOA估計

(15)

定義系數向量b=[b0,b1,…,bK-1]T及其構建的移位循環(huán)矩陣為：

(16)

對式(15)兩邊同時乘以ejJψk，可以得到如下齊次線性方程組：

(17)

結合式(16)、式(17)及陣列流型A,有如下關系式成立：

BA=0(M-K)×1。

(18)

因為陣列流型A與信號子空間US張成同一空間，所以可以用信號子空間代替陣列流型，而本文用Nystr?m方法得到逼近的信號子空間UNS=U(:,1:K)，則用UNS代替陣列流型可得BUNS=0(M-K)×1。

同理，對式(15)兩邊同時乘以e-j(J+K)ψk,可以得到關于式(17)、式(18)的另一種共軛顛倒的齊次線性方程組及關系式：

(19)

(20)

(21)

式中，vec表示矩陣的拉直運算，即把矩陣按照列的順序，一列接一列地組成一個長向量。為了方便表示，用Ψ表示關于系數b的所有拼接組成的長向量構成的矩陣,即：

b=-Ψ?γ。

(22)

得到系數b后再代入式(15),通過求解此特征多項式的K個根αk(k=1,2，…,K)得到DOA估計如下：

(23)

式中，angle表示復數取相位角運算?？梢钥闯觯岬男路f的DOA估計算法不同于現有的大多數DOA估計算法及求根DOA估計算法?，F有的求根DOA估計算法利用了Pisarenko分解，利用MUSIC算法的噪聲子空間與陣列方向向量正交的思想構建2(M-1)階的多項式求根[17-18]，即：

(24)

式中，p(z)=[1,z,…,zM-1]T，z為信號傳播到相鄰兩陣元時的相移，即z=ejψ。上式需要求解2(M-1)個關于單位圓呈鏡像對稱的根并最后選取K個在單位圓內離單位圓最近的根的相位給出DOA估計。而本文所提算法直接通過逼近的信號子空間構建一個階數為K的低階特征多項式方程并直接對求解得到的K個根提取其相位，得到DOA估計。相較于高階多項式方程求解其計算效率更高。

3 計算復雜度分析

計算復雜度的衡量一般通過矩陣間的復乘次數來衡量。通過分析可知，MUSIC及ESPRIT算法需對所有陣列輸出進行協(xié)方差的計算及對其特征分解得到信號或噪聲子空間，所需的計算復雜度分別為O(M2L)及O(M3)。當MUSIC算法進行譜峰搜索得到DOA時，計算復雜度會隨著搜索精度的提升而急劇增加，其譜峰搜索的計算復雜度為O(MKJ)，其中J代表搜索精度。而用求根MUSIC代替譜峰搜索時，相比譜峰搜索其計算量會大大下降，求根進行DOA估計時的計算復雜度為O((2(M-1))3)(若采用常規(guī)的多項式伴隨矩陣特征值分解的方法來求解多項式的根[17])。而PM算法也需進行所有陣列輸出的協(xié)方差計算O(M2L)，通過線性運算求得信號子空間的計算復雜度為O(M2K+MK2+K3)。本文所提算法無需對所有陣列輸出計算協(xié)方差矩陣及特征值分解，經分析得知通過Nystr?m方法得到的逼近信號子空間的計算復雜度主要集中在塊矩陣R11及R12的協(xié)方差矩陣計算上，其計算復雜度分別為O(Lq2)及O(MLq-Lq2)(q為定義的參數)。然后用基于Nystr?m方法構建信號子空間的計算復雜度為O(Mq2)[13]?？梢姡ㄟ^Nystr?m方法得到逼近的信號子空間總共的計算復雜度為O(MLq+Mq2)。通常有M，L?K，則獲取信號子空間這一步的計算復雜度就遠小于上述算法。而獲取低階特征方程系數b所需的計算復雜度為O(2K(K+1)(M-K)+K3)[19]，最后求根進行DOA估計的計算復雜度為O(K3)，相較于求根MUSIC，其計算復雜度在求根環(huán)節(jié)要更低。若M?K，則所提算法的整體計算復雜度主要集中在含有M的高階項，即O(MLq)，而需計算陣列輸出協(xié)方差及其特征分解算法的計算復雜度則主要為O(M3+M2L)。綜上，所提算法具有較低的計算復雜度。幾種算法的計算復雜度對比如表1所示。

表1 不同算法的計算復雜度分析

4 仿真實驗及分析

這一節(jié)進行Matlab數值仿真實驗來說明所提算法的有效性及進行性能對比分析。首先，說明驗證了所提算法(Nystr?m-LD-root)對DOA估計的有效性。而后，與文獻[11，19]所提的2種由PM算法衍生的快速估計算法PRM算法及PM-LD-rMUSIC算法進行了不同信噪比(SNR)及不同快拍數L下對算法估計性能影響的對比討論分析。

4.1 DOA估計有效性驗證

假設仿真實驗中為一陣元數M=16的均勻線陣，陣元間距λ為半波長，SNR=10 dB，快拍數L=100，并設置參數q=8，有3個遠場窄帶非相干信號分別以波達角度θ為-30°,10°,20°入射到此陣列。進行50次獨立實驗得到的DOA估計如圖2所示，可見在此仿真參數設置下所提算法能穩(wěn)定地估計出DOA且每個角度的估計基本呈一條直線，則所提算法具有一定的有效性、可行性以及穩(wěn)定性。

圖2 DOA估計

4.2 信噪比變化對估計性能的影響

用均方根誤差(RMSE)來衡量不同算法的估計性能，并做對比討論。定義RMSE如下：

(25)

假設仿真實驗中為一陣元數為M=12的均勻線陣，陣元間距d為半波長，快拍數L=500，設置所提算法的參數q取值8,10。有4個遠場窄帶非相干信號分別以波達角度θ為5°,15°,25°，40°入射到此陣列。SNR從5 dB以1 dB為步長增加至25 dB，對每一SNR進行500次蒙特卡羅仿真實驗，在此仿真參數設置下對每一算法進行仿真并統(tǒng)計所得不同SNR下的RMSE對比如圖3所示。

圖3 SNR-RMSE對比

由圖3可以看出，所有算法均隨著SNR的升高其RMSE逐漸降低，但所提算法的RMSE要始終低于PM-LD-rMUSIC算法的RMSE，說明所提算法的估計誤差要始終低于PM-LD-rMUSIC算法，估計精度較好。在較低SNR(10 dB以下)時，略遜于PRM算法的估計精度，因PRM算法結合了PM算法的快速估計信號子空間以及MUSIC-root算法的DOA估計精度高的特點，所以在低SNR時PRM算法的DOA估計誤差要稍低一些。但在高SNR(13 dB以上)時所提算法的估計精度要好于其他2種算法。在高SNR時，對所提算法的Nystr?m方法中的q參數選取不是很敏感，在q取8或10時，其RMSE在高SNR時基本一致。因此，由上節(jié)計算復雜度分析可知，當陣元數M及快拍數L固定后在高SNR情況下，可選取一個適當較小的參數q來降低算法的計算復雜度。當在較低的SNR(10 dB以下)時，較大的q值可以提升一點估計性能。

圖4 不同SNR下的估計成功概率

4.3 快拍數變化對估計性能的影響

假設仿真實驗中的陣元數M、陣元間距d、所提算法的參數q以及入射到陣列的信號波達角度θ依然與上一小節(jié)保持一致?？炫臄礚從100以100為步長增加至1 000，對每一快拍數L進行500次蒙特卡羅仿真實驗，在此仿真參數設置下對每一算法進行仿真并統(tǒng)計所得分別在SNR=10 dB以及SNR=15 dB時的不同快拍數L的RMSE曲線如圖5和圖6所示，估計成功概率如圖7和圖8所示。

圖5 信噪比為10 dB時不同快拍數的RMSE

圖6 信噪比為15 dB時不同快拍數的RMSE

圖7 信噪比為10 dB時的估計成功概率

圖8 信噪比為15 dB時的估計成功概率

由圖5和圖6可以看出，無論SNR如何變化，PM-LD-rMUSIC算法的估計誤差隨著快拍數的遞增始終要大于所提算法。在較低SNR(10 dB)時，所提算法在快拍數400以上時與PRM算法的RMSE基本保持一致。而隨著SNR升高到15 dB,所提算法在快拍數200以上其RMSE最低，估計誤差要小于其余2種算法。說明在高SNR時只需較少的快拍數就可達到較好的估計精度。從圖7和圖8的估計成功概率圖也可看出，所提算法在較低SNR(10 dB)時，在快拍數400多就可達到完全成功估計，而PM-LD-rMUSIC在快拍數1 000時還未達到完全估計成功。隨著SNR升高到15 dB,所提算法只需200快拍數就可達到完全估計成功。從圖5～圖8也可看出，所提算法基本對參數q值的選取不敏感，在q為8，10的情況下，其RMSE曲線和估計成功概率曲線基本一致，估計性能漸進一致。

5 結束語

本文提出了一種低計算復雜度的新穎的DOA估計算法。首先通過Nystr?m方法得到逼近的信號子空間，無需計算所有陣列輸出的協(xié)方差矩陣及對其特征分解，降低了計算復雜度。而后通過逼近的信號子空間直接構建關于DOA的低階特征多項式方程，此低階多項式方程的階數僅與信號源數目相同，只需求解此方程的K個根就可得到DOA估計，提升了計算效率。仿真結果及理論分析表明，所提算法可對DOA進行有效估計且計算復雜度較低，在較高SNR時只需較少的快拍數就有較低的估計誤差。