一種寬帶相干信號的二維DOA估計新方法

朱立東，鄧超升

(電子科技大學通信抗干擾技術國家級重點實驗室，四川成都611731)

0 引言

在無線移動通信系統中，同頻干擾(Co-Channel Interference，CCI)、多址干擾(Multiple Address Interference，MAI)是影響系統容量和性能的因素之一。智能天線技術為解決頻率資源不足、提高移動通信系統容量和系統服務質量提供了一個有效的解決途徑。在智能天線技術的研究中，入射信號的DOA(Direction Of Arrival)估計作為一個重要的研究內容受到廣泛關注。根據入射信號帶寬的不同，DOA估計方法可分為兩大類:窄帶信號DOA估計和寬帶信號DOA估計。常用的寬帶信號DOA估計方法大致分為信號子空間方法和最大似然估計方法。在信號子空間方法的幾個分類中，非相干信號子空間方法(Incoherent Signal-Subspace Method，ISSM)[1]是最簡單的寬帶信號處理方法，該方法在高信噪比時簡單有效，但在低信噪比下對某些頻段信號的DOA估計效果非常不理想，降低了整體的估計性能，且該方法不能處理相干源。基于角度預估計的相干信號子空間方法(Coherent Signal-Subspace Method，CSSM)[2]通過構造聚焦矩陣，將預估角度處每個頻點的陣列流形對齊到聚焦頻率處，通過陣列流形替換實現DOA估計，該方法在克服ISSM方法缺陷的同時容易引入角度預估計誤差，從而影響DOA估計性能。為進一步消除角度預估計誤差帶來的影響，近年來出現了很多無需角度預估計的寬帶波達方向估計算法[3]。文獻[4-6]提出基于波場模型的陣列流形內插方法(Array Manifold Interpolation，AMI)，將陣列流形分為僅與陣列結構和頻率有關的采樣矩陣和僅與波達方向有關的向量，并且對二維及三維陣列的波場模型分別進行建模;文獻[7]針對稀疏均勻圓陣，給出基于流形分離技術(MST)的DOA估計方法;文獻[8]將基于流形分離技術的波達方向估計方法應用到非均勻圓陣上;文獻[9]在直線陣的基礎上，給出一種采用傅里葉-勒讓德級數展開構造變換矩陣的寬帶聚焦估計方法;文獻[10]研究了不同陣列結構對空間角度估計的影響，并用仿真分析給出幾種常見陣列在低俯仰角下空間到達角的估計性能，驗證均勻圓陣在信號空間到達角估計中的優越性。

在以上文獻和方法的啟發下，提出一種均勻圓陣下寬帶相干信號的二維來波方向估計新方法。先使用多項式展開的方法實現陣列流形矢量的頻率與角度的分離，再選取簡易的與頻率有關的采樣矩陣構造聚焦矩陣，在保證聚焦前后無損失的基礎上減小計算復雜度。此外，針對均勻圓陣在俯仰角較小或較大時對空間角估計精度較差的問題提出有效的解決方法。仿真分析表明，所提方法可以獲得較優的2D-DOA估計效果，且該方法可推廣至任意的陣列結構，并明顯改善2D-DOA估計性能。

1 信號模型和算法描述

1.1 信號模型

對于任意一個含有M個陣元的幾何平面陣列，設第m個陣元的位置為(xm，ym，0)，接收N(M≥N)個位于遠場的寬帶平穩隨機信號，噪聲為加性帶通高斯白噪聲，且噪聲與源信號相互獨立。則第m個陣元的接收信號在時域中的表達式為

其中，sn(t)為第n個信號源;nm(t)為第m個陣元上的加性高斯白噪聲;τmn= (xmcosφn+ymsinφn)sinθn/c為第m個陣元接收源信號相對于參考點陣元接收信號的傳播時間延遲，其中，c為光速，(θn，φn)為第n個源信號的俯仰角和方位角，m=1，2，…，M ，n=1，2，…，N 。

對陣列的接收信號進行離散傅里葉變換，并將整個信號帶寬劃分為J個子頻帶，各子頻帶的中心頻率分別為 f1，f2…，fj，…，fJ，則接收陣列在頻域內對應于頻率fj的輸出信號可以表示為

其中，j=1，2，…，J;X(fj)=[X1(fj)，X2(fj)，…，Xm(fj)，…，XM(fj)]T，Xm(fj)為第 m 個陣元對應于頻率 fj的輸出;A(fj)=[a1(fj)，a2(fj)，…，am(fj)，…，aM(fj)]T為對應于頻率fj的陣列流形矩陣，am(fj)為導向矢量，其表達式為am(fj)= [ei2πfjτm1，ei2πfjτm2，…，ei2πfjτmn，…，ei2πfjτmN]，其中;S(fj)和 N(fj)分別為信號和噪聲在第fj頻段上的輸出矩陣，S(fj)=[S1(fj)，S2(fj)，…，SM(fj)]T，N(fj)=[N1(fj)，N2(fj)，…，NM(fj)]T。

1.2 常規算法

文獻[11]對傅里葉變換在平面、柱、球陣波束形成和DOA估計中的應用進行全面分析。對于平面陣列，其陣元坐標可通過相應的映射關系將其轉換為球坐標，轉換關系為

在球坐標下，設陣元坐標向量為 p=Rm[sinθmcosφm，sinθmsinφm，cosθm]T，同時設頻率 fj對應的第n個信號Sn(fj)的波束向量為k=[kjsinθncosφn，kjsinθnsinφn，kjcosθn]T，其 中 kj=－2π/λj=－2πfj，則信號Sn(fj)在第m個陣元處的接收信號可以表示為

因此，設 N個信號的來波方向為 Ζ(N)={(θ1，φ1)，…，(θn，φn)，…，(θN，φN)}，第m個陣元對應于頻率fj的導向矢量αm(fj)可以寫為αm(fj)=[αm(fj，θ1，φ1)，…，αm(fj，θn，φn)，…，αm(fj，θN，φN)]。其中，αm(fj，θn，φn)= exp{ikjRm[sinθnsinθmcos(φn－φm)+cosθncosθm]}，αm(fj)中的元素 αm(fj，θn，φn)可用球形貝塞爾函數及球函數的級數展開表示。

其中jξ(·)為 ξ階球形貝塞爾函數，Ylξ(·)為 ξ階l度的球函數。Ylξ(θ，φ)可以表示為

其中，Plξ(·)為ξ階l度的連帶勒讓德函數，表達式為

歸一化的球函數可以表示為

針對球貝塞爾函數，其與貝塞爾函數(第一類)之間的轉化關系如下

將式(8)和式(9)帶入式(5)，則 αm(fj，θn，φn)可以表示為

其中，C=(2ξ+1)(ξ－l)!/((ξ+l)!)。

從式(11)可知，αm(fj，θn，φn)可以分為兩項，一項僅與陣列的結構和頻率有關，另外一項與信號的來波方向有關。因此，陣列流形矩陣A(fj)可以分為兩個部分，一部分與陣列結構和頻率有關，另一部分只與波達方向有關。設僅與陣列結構和頻率有關的采樣矩陣 為 G(fj)= [G1，…，Gm，…，GM]T，Gm=[g00，g1－1，g01，g11，…，]T，故導向矢量 αm(fj)又可以表示為αm(fj)=Gm·Ψ，Ψ為關于ψln的向量。Gm的表達式中對于矩陣G(fj)，根據貝塞爾函數的性質，可以選取適當的截斷個數Nc(通常取Nc=2×(2πfjRm/c))，得到截斷后的G(fj)為M×(Nc+1)2維的矩陣，之后可選取相應的方法構造聚焦矩陣。幾種常見的方法有:基于波場模型的陣列流形內插方法[5]和旋轉信號子空間的陣列流形內插方法[12]。使用聚焦矩陣可獲得性質類似于窄帶相關矩陣的寬帶相關矩陣，進而獲得寬帶信號的DOA信息。

2 改進算法

考慮一個以圓心為參考點的均勻圓陣，設陣列的半徑為寬帶信號起始頻率對應波長的p倍，則Rm為定值，即 Rm=pλ0=pc/f0。又 kj=－2π/λj=－2πfj/c，則 kjRm=2πfjp/f0。此時，式(10)中球形貝塞爾級數項的系數可化簡得到

將(12)式帶入(10)式，αm(fj，θn，φn)可以化簡得到

通過帶入歸一化球函數的表達式，可以將αm(fj，θn，φn)進一步寫為

其中，C1為常數項，表達式為C1=(2ξ+1)(ξl)!/(4π(ξ+l)!)。通過進一步化簡與常數項合并，αm(θn，φn，fj)可以表示為

其中，C1為常數項，表達式為1)(ξ－l)!/((ξ+l)!)。

在此基礎上，考慮信號帶寬。對于寬帶信號，其相對帶寬定義為:ffoc=2(fH－fL)/(fH+fL)，其中fH和fL分別為寬帶信號的最高頻率和最低頻率，fL=f0。若將整個頻率帶劃分成J個子頻率帶，則第j(j=0，1，…，J－1)個子頻率可以表示為

令ρ=fH/fL，顯然ρ＞1，由此可以得到

在均勻圓陣圓心的軸上增加一個全向陣元，設該陣元的編號為 M+1 ，為了保證 αM+1(fj，θn，φn)的可分離性，RM+1不能為0。設RM+1=qλ0，選取q為一個較小的數，可使包含陣元M+1的立體陣趨近于平面陣，此陣列布局易于工程實現。對于該陣元，其在球坐標系下的坐標可以表示為:RM+1=qλ0，θM+1=0°，φM+1可以為任意值。當陣元處于Z軸上時，φM+1取值已無意義，但考慮到構造(fj)的需要，因此對其進行賦值。在φM+1=0°的情況下，對于陣元M+1，αM+1(fj，θn，φn)=exp(ikjRM+1cosθn)，顯然，該項可以用球貝塞爾函數及球函數的級數形式表示，并劃分成兩個部分:一部分只和該陣元的位置和頻率有關，另一部分只和來波方向和對應該頻率的衰減系數有關。

(19)式的擬合解為

其中，Uj、Vj分別為 G(f0)G(fj)H的左、右奇異矢量組成的矩陣。求得聚焦矩陣Tj后，進而可求得對應于每個頻率的聚焦變換后的相關矩陣 ^R(fj)，其表達式為

對聚焦后的數據相關矩陣進行相加求平均，可得到相關矩陣R，其表達式為

由式(15)可知，R可視為一個窄帶相關矩陣，且包含原寬帶信號的所有信息。因此，可按照窄帶信號DOA估計方法的步驟對該相關矩陣進行處理，從而得到寬帶信號的DOA估計信息。此外，由于R的構成過程中對所劃分的各個頻率分量進行平滑處理，因而R具有解相干性能，能處理寬帶相干信號。

根據以上理論分析，采用文中方法構造聚焦矩陣的寬帶相干二維DOA估計方法的步驟可總結如下:

(1)將陣列接收信號的快拍數據分為K段，每段含有L個快拍，同時將頻帶劃分為J個頻段;

(2)分別對每段數據進行DFT變換，得到J個頻段的陣列接收數據;

(3)分別對J個頻段的接收數據構造頻域采樣空間相關矩陣:

(4)對G(f0)G(fj)H進行奇異值分解，求得Uj、Vj，再利用式(20)求出聚焦變換矩陣Tj;

(5)根據式(22)得到組合加權后的相關矩陣，并進行特征值分解，得到噪聲子空間和信號子空間;

(6)利用二維MUSIC算法進行來波方向估計。

所述算法在計算聚焦變換矩陣Tj的過程中，盡可能剔除與頻率無關的乘積項來化簡采樣矩陣(fj)，并且巧妙地結合均勻圓陣的特點，進一步化簡(fj)，在提高聚焦效果的同時減小了計算量。此外，根據空間角估計精度隨俯仰角的變化情況，可適當調整陣列結構，能夠改善估計性能。

3 數值仿真分析

假設陣列的接收信號為遠場寬帶信號，中心頻率為f0=100 Hz，信號相對帶寬為40%，采樣頻率fs=256 Hz，噪聲為加性零均值帶通高斯白噪聲，噪聲與源信號相互獨立，陣列接收數據的快拍數為4096個，把數據分為32段，每段含有128個快拍。信號帶寬為40 Hz，可把整個帶寬劃分為41個子頻帶，選擇中心頻率f0為聚焦頻率。根據貝塞爾函數的性質，選取截斷點數Nc為9。為便于觀察和分辨，將常規多項式展開方法(general polynomial expansion method)簡記為GPEM;文中方法(improved polynomial expansion method)簡記為IPEM;改進后的均勻圓陣(improved uniform circular array)簡記為IUCA。

仿真1:驗證文中算法與常規算法的均方根誤差。假設進行n次獨立重復實驗，N個信號源的空間角的均方根誤差定義為式中θd和φd為信號的真實俯仰角和方位角，^θdi和^φdi為估計值。設接收陣列的陣元個數為11，陣列構型分別為UCA和IUCA，仿真使用兩個完全相干的遠場時域平穩寬帶信號，其均值為0，信號到達角分別為(40°，40°)、(20°，200°)，進行 100 次蒙特卡洛仿真，統計文中算法與常規算法在不同信噪比條件下的空間角的均方根誤差，對比曲線如圖1所示。

仿真2:統計文中算法與常規算法的成功分辨率。使用兩個完全相干的信號進行仿真，其來波方向分別為(40°，50°)和(60°，75°)，其他仿真條件同仿真 1。當兩個信號的俯仰角和方位角的估計值均小于0.2°時視為成功估計，進行200次蒙特卡洛仿真，成功分辨性能曲線如圖2所示。

仿真3:驗證文中方法在其他任意陣列上的空間角估計性能。仿真中使用的任意陣列的布局如圖3所示，對任意陣的改進操作為:將任意陣中編號為M的陣元置于參考點上方合適位置，總體陣元個數不改變。兩個相干信號的來波方向分別為(20°，50°)和(35°，55°)，其他仿真條件同仿真1。得到的均方根誤差性能曲線如圖4所示。

仿真4:驗證文中算法對兩個來波方向相近的寬帶相干信號的分辨性能。設兩個相干信號來波方向分別為(10°，200°)和(15°，200°)，文中算法在11 個陣元組成的IUCA的陣列構型下對兩個信號的俯仰角分辨性能如圖5所示。

圖1 不同方法在不同信噪比下的空間角估計性能

圖2 不同方法的成功分辨概率

圖1分別對比了改進的算法和常規算法分別在11陣元均勻圓陣、10陣元均勻圓陣與圓心軸上的陣元組成的11陣元陣列下的空間角估計性能。無論是使用均勻圓陣(UCA)或是改進的均勻圓陣(IUCA)，IPEM算法在空間角的均方根誤差上的性能都優于GPEM。在UCA下，兩種算法的均方根誤差在13～20 dB內大致趨于穩定，上下波動在0.01°范圍以內。其中，IPEM算法的均方根誤差在20 dB時最小，大小為0.1227°;GPEM算法的均方根誤差在20 dB時最小，大小為0.1468°。兩種算法的均方根誤差的差值在0.0220°和0.0707°之間波動。在IUCA下，使用兩種算法得到的DOA估計的均方根誤差逐漸減小，隨著信噪比從0 dB到20 dB的變化，IPEM算法得到的均方根誤差從0.1441°降低到0.054°，GPEM算法得到的均方根誤差從0.2031°降低到0.079°。兩種算法的均方根誤差的差值在0.0250°和0.0589°之間波動。此外，GPEM算法在ICUA下的均方根誤差相對于在UCA下降低的幅度大于0.04°。同時，IPEM算法在IUCA下的性能提升幅度同樣大于 0.04°，并且最大達到0.0704°。

由圖2可以看出，在UCA下，IPEM算法在信噪比為－1 dB時的成功分辨概率為81.5%，GPEM算法的成功分辨率為73%，IPEM算法在GPEM算法的基礎上成功分辨概率提升了8.5%。在信噪比為－5～14 dB，改進后的 IPEM與常規的多項式展開方法GPEM相比成功分辨概率都有所提升，提升幅度范圍為1%～8.5%。在陣列構型為IUCA的情況下，IPEM算法和GPEM算法在信噪比高于6 dB條件下的成功分辨概率都為100%。在信噪比為－5～5 dB，IPEM算法的成功分辨概率明顯高于GPEM算法的成功分辨概率，提升幅度范圍為0.25%～15%。其中，在信噪比為－5 dB時提升幅度為15%。總體來說，IPEM算法在相同陣列構型和相同信噪比的條件下，成功分辨概率高于GPEM算法。此外，同一算法在不同陣列構型下的成功分辨概率也有所不同。在信噪比為－5～15 dB，兩種算法在ICUA陣列下的成功分辨率比在UCA陣列下有所提升，IPEM算法最高提升了12.75%，GPEM算法最高提升了11.25%。

分析出現圖1和圖2中統計結果的原因可知，兩種算法的區別僅在于用于構造聚焦矩陣的采樣矩陣G(fj)不同，在構造G(fj)的過程中，重點在于將不同頻率上的陣列流形統一到聚焦頻率上的陣列流形，IPEM算法盡可能剔除與頻率無關的項，并巧妙結合均勻圓陣的特性，進一步將 ^G(fj)化簡，提高了聚焦效果。在改進的陣列結構下，兩種算法的性能在均勻圓陣的基礎上提升明顯?？臻g角的均方根誤差不僅與算法有關，還與寬帶信號的入射角度有關。仿真(1)中使用的俯仰角20°屬于較低的俯仰角，如文中分析，該角度處sin(θ)的值較小，但sin(θ)的變化率較大，此時對方位角的估計精度不高，相當于有很小的陣列孔徑，改進陣列在原陣列的基礎上加上了cos(θ)項，該項的加入改善了低俯仰角下的方位角的精度，進而提升了空間角估計性能。在仿真(2)中使用的相干信號中，其中有一個信號的俯仰角為60°，屬于較高俯仰角，該俯仰角的正弦值sin(θ)較大，但是sin(θ)的變化率很低，此時對俯仰角的估計精度不高，影響了算法的成功分辨性能。改進后的陣列有效地解決了均勻圓陣的不足，明顯提高了陣列的成功分辨性能。

圖3 某任意幾何陣列示意圖

圖3中將原點作為參考點，各點的坐標表示各個陣元相對于參考點的位置，x軸和y軸的單位長度為寬帶相干信號起始頻率所對應的波長λ0。在圖3所示的陣列構型下統計IPEM算法得到的兩個寬帶相干信號的俯仰角與方位角的均方根誤差曲線如圖4。由圖4可以看出，在兩個不同的陣列構型下，IPEM算法在俯仰角上的均方根誤差變化不明顯，尤其是當信噪比高于9 dB以后俯仰角的變化幅度都在0.004°范圍內。但在方位角上的均方根誤差變化明顯，最大差值達到0.0516°，在改進陣列下，均方根誤差平均減小了0.0417°。以上數據表明，方法可以推廣至任意陣列，在改進陣列的基礎上結合文中算法IPEM能夠提高寬帶相干信號的DOA估計性能。根據文中分析，對于來波方向分別為(20°，50°)和(35°，55°)的兩個相干信號，其方位角的估計精度較低，如圖4所示。若存在參考點上方的陣元，對具有較低俯仰角的接收信號而言，該陣列結構對俯仰角的估計精度影響不明顯，但是可以有效地降低方位角的估計誤差，從而提高總體的DOA估計性能。

圖4 文中算法在任意陣列下的性能

圖5 文中算法的分辨性能

由圖5可以看出，在信噪比為0 dB、5 dB和10 dB的條件下，文中算法對仿真所設的兩個相干寬帶信號的俯仰角估計的平均誤差分別為0.205°、0.135°和0.06°。通過仿真結果可以看出，采用文中方法進行寬帶相干信號的二維DOA估計時，在信噪比為0 dB以上能夠成功分辨出空間中相隔較近的相干源。以方位角相同、俯仰角相隔5°的兩個寬帶相干信號為例，隨信噪比的增加，信號的空間譜峰的幅度不斷增大。在信噪比從0 dB增長至10 dB的過程中，譜峰高度從70.52和102.8增大至1244和1754，且譜峰變得更加尖銳，譜峰所對應的俯仰角與真實俯仰角的平均誤差從0.205°逐漸減小至0.06°，說明文中算法隨著信噪比的不斷增加分辨性能也逐漸提高。

4 結束語

均勻圓陣具有幾乎無方向性、二維性和無模糊范圍廣等特點，應用均勻圓陣進行寬帶測向能夠得到良好的效果。以均勻圓陣為接收模型，首先使用多項式展開的方法對陣列流形矢量進行頻率與角度的分離，通過分析陣列流形矢量的級數展開式，結合均勻圓陣的特點，構造較簡易的與頻率相關的采樣矩陣構造聚焦矩陣，可以達到減少計算量的目的。在此基礎上，對陣列構型進行一定調整可明顯改善均勻圓陣或其他陣列在俯仰角較大時或較小時對空間角角估計精度較差的問題。數值仿真分析顯示，文中方法在寬帶相干信號的2D-DOA估計中，能夠降低信號空間到達角度估計的均方根誤差，提高相干信號的成功分辨概率。同時，該算法對空間中相隔較近的相干源具有良好的分辨性能，且該方法也可推廣到任意的陣列結構。

[1] Wax M，Shan T J，Kailath T.Spatio-temporal spectral analysis by eigenstructure methods[J].IEEE Transactions on Acoustics Speech and Signal Processing，1984，32:817－827.

[2] M A Doron，A J WEISS.On focusing matrices for wide-band array processing[J].IEEE Trans on Signal Process-ing，1992，40(6):1295－130.

[3] Bo L，Xiong J，Luo L.A novel wideband DOA estimation method using direction-free focusing matrix[C]//Computer Science and Network Technology(ICCSNT)，2013 3rd International Conference on.IEEE，2013:1065－1069.

[4] M A Doron，E.Doron.Wavefield modeling and array processing.Ⅰ:Spatial sampling[J].IEEE Trans Signal Process，1994，42(10):2549－2559.

[5] M A Doron，E Doron Wavefield modeling and array processing.Ⅱ:Algorithms[J].IEEE Transaction Signal Process，1994，42(10):2560－2570.

[6] M A Doron，E Doron.Wavefield modeling and array processing.Ⅲ:Resolution capacity[J].IEEE Trans Signal Process，1994，42(10):2571－2580.

[7] 潘捷，周建江，汪飛.基于流形分離技術的稀疏均勻圓陣快速DOA估計方法[J].電子與信息學報，2010，32(4):963－966.

[8] 盧海杰，章新華，熊鑫.流形分離在非均勻圓陣上的應用[J].兵工學報，2011，32(9):1113－1117.

[9] 林晉美，侯衛民，涂英，等.傅立葉-勒讓德級數展開的寬帶聚焦[J].聲學技術，2007，12，26(6):1235－1239.

[10] 顧建峰.空間角度估計的陣列結構及其算法的研究[D].成都:電子科技大學，2004.

[11] 陳豐.Fourier變換在平面、柱、球陣波束形成和DOA估計中的應用[D].杭州:浙江大學，2006.

[12] 李剛，趙春暉，寧海春.用于寬帶測向的廣義陣列流行內插算法[J].彈箭與制導學報，2008，28(1):277－279.