均勻圓陣主特征矢量分析解相干算法*

倪柳柳，陳 輝，倪萌鈺，王曉戈

(空軍預警學院， 湖北 武漢 430019)

波達方向(Direction-Of-Arrival, DOA)估計一直是無線電通信、電子對抗偵察、聲吶以及雷達系統等諸多領域的一個重要問題，在進行DOA估計之前需要確定陣列形狀，常用陣型有均勻線陣(Uniform Linear Array，ULA)、L型陣列及均勻圓陣(Uniform Circular Array，UCA)等，其中，均勻圓陣因其可提供360°的全方位無模糊定位以及各向同性性質而越來越得到重視。同時在實際波達方向估計中，由于反射或折射等引起的傳輸多徑現象，或敵方設置的有源電磁干擾等，相干信號廣泛存在，因此對于相干信號的檢測與估計也是DOA估計中的重要研究方向。因此本文將研究適用于均勻圓陣的解相干算法。

目前，適用于均勻圓陣的解相干算法主要有三大類。一是虛擬陣變換，其思想就是通過陣列變換將圓陣轉換為虛擬均勻線陣[1-3]或者多個相同間距的均勻圓陣[4-6]，然后再利用空間平滑等方法進行解相干處理。如模式空間平滑(Method Of Direction Estimation-Space Smoothing, MODE-SS)算法[1]，就是將均勻圓陣通過模式空間變換轉換成一個同孔徑等距均勻線陣，再進行空間平滑處理。模式空間托普利茲(Method Of Direction Estimation-TOEPlitz, MODE-TOEP)[2]是在模式空間變換后重構一個Toeplitz矩陣，這個矩陣的秩只和DOA信息有關，因此可以解相干。而模式空間差分(Method Of Direction Estimation-DIFFerence, MODE-DIFF)[3]算法將MODE-TOEP算法和差分算法相結合，達到同時估計非相關和相干信源的目的。模式空間中心對稱(Method Of Direction Estimation-Central Symmetry, MODE-CS)[5]利用均勻圓陣的中心對稱性對模式空間數據進行共軛平均，通過相關性構造Hermitian Toeplitz矩陣，實現圓陣相干信源的方位估計。這類算法的優點是解相干的信號源數較多；但其缺點為陣列的變換會帶來誤差，導致DOA估計精度不高。二是采用多維搜索類算法[7-11]進行估計。這類算法不需要進行解相干預處理，就可以直接進行DOA估計，如基于輪換投影的最大似然(Alternative Projection-Maximum Likelihood, AP-ML)算法[10]、加權信號子空間擬合(Weighted Signal Subspace Fitting, WSSF)算法[8,11]等。這類算法的優點就是算法的估計精度高；缺點是運算過程中要用到多維搜索，運算量巨大。三是稀疏重構類DOA算法[12-18]，這類算法利用恢復稀疏信號的思想進行算法設計。Malioutov等將DOA估計算法轉變為基于L1范數的優化問題，提出了基于L1范數的奇異值分解(L1-Singular Value Decomposition, L1-SVD)算法[12]；文獻[13]中的基于L1范數的矩陣協方差稀疏重構(L1-Sparse Representation of Array Covariance Vectors, L1-SRACV)算法利用了協方差矩陣的稀疏特性，并通過求解二階錐優化問題來實現DOA估計；Yang等則提出了離格稀疏貝葉斯(Off Grid-Sparse Bayesian Inference, OG-SBI)DOA估計方法[14-15]，也就是首先設置較粗的網格，而在估計過程中通過不停迭代將真實信源附近的網格不斷逼近真實角度，最終實現DOA估計。這類算法的優點是無須信源數信息，對相干信源不敏感，估計精度較高。缺點為算法都是有偏估計，其DOA估計值相對真實值的偏差取決于網格的劃分，而網格劃分太小的話復雜度將非常高。

近年來，一種基于均勻線陣的新的解相干DOA算法[19-20](Principal-eigenvector Utiliztion for Modal Analysis, PUMA)被提出，該算法利用均勻線陣具備的范德蒙德結構來設計DOA算法，能夠較好地估計相干信源的DOA。受此啟發，本文提出了一種基于虛擬均勻線陣的主特征矢量分析算法(Virtual uniform Linear array Principal Eigenvector Analysis，VLPEA)，算法在模式空間變換的基礎上，分析虛擬均勻線陣主特征矢量的特性，通過求解線性預測算子最終成功估計出信源的方向。下節的分析和實驗證明了本文算法估計精度高，復雜度低。

1 均勻圓陣陣列模型及預處理

本文采用均勻圓陣為陣列模型，其結構如圖1所示。

圖1 均勻圓陣陣列模型Fig 1 Array model of UCA

M個各向同性的陣元均勻分布在半徑為r的圓陣上，遠場窄帶相干信號以方位角度θ和俯仰角度φ入射到各個陣元，θ為信源和陣列中心的連線在xoy面上的投影與x軸的夾角，φ為信源和陣列中心連線與xoy平面之間的夾角，其中θ∈[0,2π]，φ∈[0,π]。

假設空間有K個波長為λ的遠場窄帶信號入射到UCA上，以陣列中心為參考點，第m個陣元位置為pm=(xm,ym,0)，xm=rcosψm，ym=rsinψm，ψm=2π(m-1)/M，ψm為第m個陣元和圓心連線與x軸的夾角，第k個信號坐標為γk=(cosφkcosθk,cosφksinθk,sinφk)，因此陣列模型導向矢量為

a(θk,φk)=[e-jω0τ1k，e-jω0τmk，…,e-jω0τMk]T

(1)

式中：τmk=cosφkcos(θk-ψm)r/c為第k個信號在第m個陣元相對于參考陣元的時延k∈[1,K]；ω0=2πc/λ，c為光速，λ為信號波長。因此流型矩陣為

A=[a(θ1,φ1),a(θk,φk),…,a(θK,φK)]

(2)

第t次快拍時的接收數據矢量為

x(t)=As(t)+n(t),t=1,…,N

(3)

式中，x(t)∈CM，u(t)∈CK和n(t)∈CM分別表示第t快拍數下的陣列接收數據、信源矢量和噪聲矢量，N為總快拍數。

本文只討論所有信源與xoy平面共面的情況，此時φk=0，接收數據變為

x(t)=ACs(t)+n(t)

(4)

式中，AC(θ)=[aC(θ1),aC(θ2),…,aC(θK)]，導向矢量aC(θk)=a(θk,0)。

綜上可知，均勻圓陣的陣列流形矩陣不符合范德蒙德結構，大部分解相干DOA估計方法并不適用，因此考慮利用模式空間變換法將陣列流型進行變換，得到符合范德蒙德結構的陣列流型，為下一步的解相干奠定基礎。構造模式空間變換矩陣為

(5)

式中：J=diag{j-QJ-Q(-β),…,jQJQ(-β)}，Jq(-β)為q階第一類貝塞爾函數，β=2πr/λ，Q=?β」為均勻圓陣可以激發的最大相位模式，這里?·」表示向下取整；F=[w-Q,w-Q+1，…,wQ]，wq=[1,e-j2πq/M,…,e-j2πq(M-1)/M]，q=-Q,…,Q。

當取M>2h+1時，可以用構造的模式空間變換矩陣得到模式空間中虛擬均勻線陣的接收數據矢量：

(6)

(7)

從式(7)可以看出，變換后的陣列流型具有范德蒙德結構，即等效為陣元數為D=2Q+1的虛擬均勻線陣，變換后的陣列數據協方差矩陣為

(8)

(9)

由式(9)可知，模式變換后的各通道噪聲功率不一致，因此必須使用廣義特征分解求解子空間。設λd(d=1,2,…,D)為矩陣束(RL,TTH)廣義特征分解后得到的所有特征值,按從大到小排列，ud為對應的特征矢量，前P個大特征值對應的主特征矢量組成矩陣ULS=[u1,u2,…,uP]，假設有K1個相干信源，則P=K-K1+1，但是其張成的子空間并不是虛擬線陣的信號子空間，對其進行進一步的處理，得到修正的主特征矢量矩陣：

VS=TTHULS

=[v1,v2,…,vP]

=[TTHu1,TTHu2,…,TTHuP]

(10)

同上可得修正的噪聲子空間：

Vn=[vP+1,vP+2,…,vD]

=[TTHuP+1,TTHuP+2,…,TTHuD]

(11)

2 虛擬均勻線陣主特征矢量分析算法

由文獻[21]證明：

span(VS)=span(AL)

(12)

式(12)表明信號子空間VS和虛擬線陣的流行矩陣張成同一空間，而AL具有范德蒙德結構，因此考慮利用線性預測(Linear Prediction, LP)算法的原理，即vk(k=K+1,K+2,…,P)中任一元素和此元素之前K個元素線性相加為0，來求解線性系數對應的角度。

2.1 算法推導

根據LP原理，得到

(13)

式中，[vk]l為vk中第l個元素，ci為之前第i個元素對應的LP系數。因此第k個角度θk包含在如下多項式中：

(14)

式中，zk=ejθk。式(13)用矩陣形式可以表示為

Gkc=hk

(15)

式中

(16)

c=[c1,c2,…,cK]T

(17)

hk=-[[vk]K+1,[vk]K+2,…,[vk]D]T

(18)

Gc=h

(19)

式中

(20)

(21)

令εk=Gkc-hk，可以得到

(22)

顯而易見，理想情況下，式(19)是經典線性系統，可以用最小二乘(Least Squares, LS)法求解c，其閉式解為

cLS=G?h

(23)

(24)

(25)

同理得

(26)

(27)

(28)

式中

(29)

(30)

在無誤差條件下

ε=vec([ε1,ε2,…,εP])=vec(B(c)VS)=0(D-K)P

(31)

式中

(32)

其中，toep(a,b)表示由向量a和b構成的Toeplitz矩陣，然后將誤差考慮進去，式(31)變為

(33)

(34)

文獻[19]指出主特征矢量誤差之間具有如下特性：

(35)

式中，δij為沖激函數。由式(35)可得

(36)

(37)

由式(34)、式(36)和式(37)得

(38)

(39)

(40)

(41)

(42)

(43)

2.2 算法流程

綜上所述本文提出VLPEA算法，算法步驟如下所示。

(44)

式中，“∠”為求角度符號

2.3 算法分析

MODE-TOEP算法利用協方差矩陣的某一列重構具有Toeplitz結構的矩陣，因此可以估計相干信源，然而也是因為上述重構方法，其信息利用并不充分，因此其估計精度無法達到最佳。MODE-DIFF算法使用了差分算法分別估計相干信源和獨立信源，然而差分的同時也去掉了分相干信源的部分信息，因此同樣無法將估計精度提升到最優。改進特征分解(Excitation Mode EigenValue Decomposition, EMEVD)算法[21]是在上述兩種算法的基礎之上再加上一次特征分解，同樣無法解決信息利用不全的問題。而本文的VLPEA算法通過對式(16)進行WLS求解，該矩陣包含了信源的所有信息，因此本文算法精度要比EMEVD算法、MODE-TOEP算法和MODE-DIFF算法高，同時文獻[19]證明了式(19)的WLS代價式在信噪比或者快拍數趨近于無窮大時可以等價于ML算法，因此VLPEA算法估計精度很高。

3 仿真與分析

為驗證VLPEA算法的有效性，本節使用MATLAB進行如下仿真實驗，所有實驗均采用15陣元的均勻圓陣陣列結構(即M=15)，陣元半徑為0.7倍的波長。在仿真實驗中，將本文算法與MODE-TOEP算法、MODE-DIFF算法、MODE-CS算法，ML算法以及OG-SBI算法做比較，其中ML算法采用文獻[10]中的AP-ML算法，OG-SBI算法中的終止迭代條件為ξ≤1×10-6，最大重復迭代次數為1 000，初始網格為2°。在單次實驗中每個真實角度與估計角度的偏差小于0.4倍的相鄰兩角度間隔時，則認定為算法分辨成功，完全均方根誤差(Completed Root Mean Squared Error, CRMSE)定義為成功分辨DOA信號的均方根誤差(Root Mean Square Error, RMSE)，如式(45)所示。

(45)

3.1 實驗1： 比較算法在信噪比變化時的性能

實驗1比較了不同算法的CRMSE和分辨概率隨信噪比變化的性能改變情況。假設有4個等功率信號入射到陣列上，DOA分別為θ1=10°，θ2=25°，θ3=50°，θ4=70°，前兩個信號源相干，快拍數取200，信噪比(Signal to Noise Ratio, SNR)的變化范圍為-5～20 dB，步長為2 dB，進行2 000次蒙特卡洛獨立重復實驗，實驗結果如圖2所示。從圖2(a)可以看出，ML算法的信源CRMSE是最小的，本文提出的VLPEA算法的性能要比EMEVD算法、MODE-TOEP算法、MODE-CS算法和MODE-DIFF算法好，這與上文的分析是一致的，同時在信噪比小于15 dB時，VLPEA算法的CRMSE比OG-SBI算法小，當信噪比更大時OG-SBI算法性能則有所提升。由圖2(b)可以看出，當SNR大于1時，VLEPA算法的成功概率達到100%，而OG-SBI、EMEVD算法、MODE-TOEP算法和MODE-DIFF算法的成功概率在SNR大于10時才能達到100%。

(a) 不同DOA估計算法CRMSE隨信噪比SNR變化(a) CRMSE of different DOA estimationalgorithms versus SNR

(b) 不同算法的DOA估計成功概率隨信噪比SNR變化(b) Success probability of different DOA estimation algorithms versus SNR圖2 不同信噪比情況下的算法性能分析Fig 2 Performance analysis of the algorithms under different SNR

3.2 實驗2： 比較算法在快拍數變化時的性能

圖3顯示的是不同算法的CRMSE和成功概率隨快拍數變化的情況，SNAP即為快拍數，此實驗中SNR固定為20 dB，快拍數變化范圍為20～300，步長為20次快拍數，其他參數的設定和實驗1一致。可以看出，隨著快拍數的增加， VLPEA算法和其他幾種算法的CRMSE都呈現降低的趨勢，其中依然是ML算法的CRMSE始終最小，OG-SBI算法其次，VLPEA算法的CRMSE十分接近上述兩種算法，而EMEVD算法、MODE-TOEP算法、MODE-CS算法和MODE-DIFF算法的CRMSE則比上述三種算法大接近一個數量級，在圖3(b)中，MODE-TOEP算法和MODE-DIFF算法同樣需要較多的快拍數才能保持性能穩定。

(a) 不同算法的DOA估計CRMSE隨快拍數SNAP的變化(a) CRMSE of different DOA estimation algorithms versus SNAP

3.3 實驗3： 比較算法在角度間隔變化時的性能

本次實驗比較各個算法的分辨性能。設定信噪比為20，快拍數為200，兩個相干信源入射到陣列上，第一個角度θ1固定為10°，第二個角度θ2從12°變化到30°，其他參數和前面的實驗一致，結果如圖4所示。圖4的橫坐標為兩個角度的間隔變化，由圖4(b)可以看出，在角度間隔大于5°時，ML算法和VLPEA算法的成功概率能到100%，而EMEVD算法、MODE-TOEP算法、MODE-CS算法和MODE-DIFF算法在信噪比大于10時成功概率才能到100%，因此圖4(a)觀察所有算法成功概率等于100%時的CRMSE，發現ML算法、OG-SBI算法以及VLPEA算法的CRMSE都維持在一個比較低的水平。

(a)不同算法的DOA估計CRMSE隨角度間隔的變化(a) CRMSE of different DOA estimationalgorithms versus spacing

(b) 不同算法的DOA估計成功概率隨角度間隔的變化(b) Success probability of different DOA estimation algorithms versus spacing圖4 不同角度間隔下的算法性能分析Fig.4 Performance analysis of the algorithms under different spacing

3.4 實驗4： 算法的時間復雜度分析

本次實驗比較各個算法的時間復雜度。快拍數固定為200，兩個入射角度分別為10°和25°，信噪比固定為20，陣元數不變，表1所示為不同算法在陣元數變化時的單次平均運行時間，可以看出OG-SBI算法的耗時相當長，這是因為為達到較高的精度，其算法中的迭代門限取得很小，導致重復次數非常高， ML算法也需要大約1 s才能完成DOA估計，而EMEVD算法、MODE-TOEP算法、MODE-DIFF算法，MODE-CS算法以及VLPEA算法都能在極短時間內分辨出信源，且VLPEA算法運算時間最短，達到了毫秒級，說明了VLPEA算法的計算量很低，與前文的分析一致。

表1 不同算法單次平均運行時間比較

4 結論

本文提出了一種均勻圓陣的相干信源波達方向估計算法——VLPEA算法。該算法首先通過模式空間變換技術將均勻圓陣轉換為模式空間中的虛擬均勻線陣，在此基礎上對其數據協方差矩陣廣義特征分解，得到新的主特征矢量矩陣；其次根據均勻線陣的范德蒙德結構，構造包含信源DOA信息的線性多項式；最后通過加權最小二乘法重復迭代得到線性預測系數成功求得相干信源的波達方向。

通過理論分析和仿真結果，可以得出VLPEA算法的估計精度明顯優于MODE-TOEP和MODE-DIFF算法，且與ML算法和OG-SBI算法接近，但是復雜度遠低于上述兩種算法，因此VLPEA算法具有較好的估計精度、分辨力且實時估計更可靠。