一種用于矩形陣列的二維波達(dá)方向估計(jì)方法

王劍書，樊養(yǎng)余，杜 瑞，呂國(guó)云

(西北工業(yè)大學(xué) 電子信息學(xué)院，陜西西安，710129)

波達(dá)方向(Direction Of Arrival, DOA)估計(jì)是陣列信號(hào)處理領(lǐng)域的一個(gè)重要研究方向，并廣泛應(yīng)用于大量實(shí)際場(chǎng)景，比如雷達(dá)、聲吶、麥克風(fēng)陣列和通信系統(tǒng)等。通過(guò)傳感器陣列接收單個(gè)或多個(gè)目標(biāo)源的信號(hào)，形成陣列信號(hào)，而波達(dá)方向估計(jì)則是通過(guò)對(duì)陣列信號(hào)進(jìn)行處理求得目標(biāo)源方向的技術(shù)方法。早期的波達(dá)方向估計(jì)一般通過(guò)波束形成實(shí)現(xiàn)，這類方法實(shí)現(xiàn)簡(jiǎn)單，應(yīng)用廣泛，但對(duì)角度相鄰的源的分辨能力欠佳?，F(xiàn)代高分辨波達(dá)方向估計(jì)方法主要包括子空間類方法與稀疏重構(gòu)方法。子空間類方法主要包括多重信號(hào)分類(MUltiple SIgnal Classification, MUSIC)[1]和旋轉(zhuǎn)不變參數(shù)估計(jì)技術(shù)(Estimation of Signal Parameters via Rotational Invariance Techniques, ESPRIT)[2]等，稀疏重構(gòu)方法主要包括l1范數(shù)最小化[3-5]和稀疏貝葉斯學(xué)習(xí)[6]等。

上述方法主要基于一維波達(dá)方向估計(jì)研究，盡管部分方法可以使用二維(Two-Dimensional, 2-D)陣列，對(duì)于二維波達(dá)方向估計(jì)(同時(shí)估計(jì)方位角和俯仰角)，大部分方法需要對(duì)空間角度進(jìn)行網(wǎng)格化，對(duì)于角度的搜索也擴(kuò)展到二維，這大大增加了計(jì)算負(fù)擔(dān)。利用特殊結(jié)構(gòu)的二維陣列，不少高效的二維波達(dá)方向估計(jì)方法被提出。對(duì)于均勻矩形陣列(Uniform Rectangular Array, URA)，文獻(xiàn)[7]提出了基于酉ESPRIT的二維波達(dá)方向估計(jì)方法。對(duì)于均勻圓陣，文獻(xiàn)[8-9]通過(guò)對(duì)其陣列流型進(jìn)行貝塞爾函數(shù)重新表達(dá)，提出有效的類ESPRIT算法。對(duì)于L形陣列，文獻(xiàn)[10-12]利用兩個(gè)線陣的互相關(guān)矩陣，構(gòu)造了有效的二維波達(dá)方向估計(jì)方法。文獻(xiàn)[13-14]使用均勻矩形陣列或稀疏矩形陣列(Sparse Rectangular Array, SRA)，利用均勻矩形陣列協(xié)方差矩陣的二階特普利茨結(jié)構(gòu)，重構(gòu)了其協(xié)方差矩陣，并通過(guò)二維ESPRIT方法進(jìn)行波達(dá)方向估計(jì)。上述使用特殊二維陣列的波達(dá)方向估計(jì)方法無(wú)需對(duì)空間角度網(wǎng)格化，可以直接估計(jì)方位角和俯仰角。這種無(wú)網(wǎng)格的波達(dá)方向估計(jì)，稱為無(wú)格波達(dá)方向估計(jì)。上述使用均勻矩形陣列或稀疏矩形陣列的方法有以下缺點(diǎn)：文獻(xiàn)[7]的方法不能用于稀疏矩形陣列情形；文獻(xiàn)[13]的快速無(wú)格最大似然方法(Fast Grid-less Maximum Likelihood, FGML)本質(zhì)為加權(quán)的最小二乘法，由于未對(duì)協(xié)方差矩陣進(jìn)行正定約束，在稀疏矩形陣列下估計(jì)性能一般；文獻(xiàn)[14]使用的低秩矩陣重構(gòu)(Low-Rank Matrix Reconstruction, LRMR)法，以核范數(shù)進(jìn)行凸松弛，求解結(jié)果并不是最優(yōu)。近年來(lái)發(fā)展起來(lái)無(wú)格一維波達(dá)方向估計(jì)方法(包括原子范數(shù)最小化(Atomic Norm Minimization, ANM)法[15-17]、加權(quán)原子范數(shù)法[18]和無(wú)格稀疏迭代協(xié)方差估計(jì)(Grid-Less SParse Iterative Covariance-base Estimation, GL-SPICE)方法[19-21]等)，表現(xiàn)出了非常良好的波達(dá)方向估計(jì)性能，尤其是文獻(xiàn)[18]提出的使用log-det函數(shù)的稀疏測(cè)度，若將其應(yīng)用于均勻矩形陣列或稀疏矩形陣列情形，可以進(jìn)一步提高二維波達(dá)方向估計(jì)的性能。

筆者提出一種基于均勻矩形陣列或稀疏矩形陣列的二維無(wú)格波達(dá)方向估計(jì)方法。首先對(duì)協(xié)方差矩陣的二階特普利茨結(jié)構(gòu)進(jìn)行更簡(jiǎn)潔的表達(dá)；然后根據(jù)文獻(xiàn)[18]提出的稀疏測(cè)度，提出該二階特普利茨矩陣的重構(gòu)方法；最后通過(guò)二維ESPRIT方法[22]進(jìn)行二維波達(dá)方向估計(jì)。筆者提出的方法需要多次求解半定規(guī)劃(Semi-Definite Programming, SDP)問(wèn)題，計(jì)算復(fù)雜度相對(duì)較高，但由于二階特普利茨矩陣的重構(gòu)結(jié)果比其他方法有更好的稀疏性，故能獲得更好的波達(dá)方向估計(jì)結(jié)果。在仿真實(shí)驗(yàn)中，分別在均勻矩形陣列和稀疏矩形陣列下，對(duì)不同快拍數(shù)、信噪比(Signal-to-Noise Ratio, SNR)和相鄰源角度間隔的均方根誤差(Root-Mean-Square Error, RMSE)結(jié)果進(jìn)行分析，驗(yàn)證了筆者提出的方法具有良好的波達(dá)方向估計(jì)性能。

1 均勻矩形陣列和稀疏矩形陣列的接收信號(hào)模型

使用均勻矩形陣列或稀疏矩形陣列，假設(shè)各陣元各向同性。為方便描述，先考慮均勻矩形陣列。如圖1所示，將均勻矩形陣列置于xOy平面，x軸和y軸上陣元數(shù)目分別為Mx和My，記M0=MxMy；該陣列接收K個(gè)獨(dú)立窄帶源信號(hào)，最小陣元間距設(shè)為窄帶信號(hào)中心頻率的半波長(zhǎng)λ/2，第k個(gè)入射信號(hào)的方位角和俯仰角分別為θk∈[-π,π]和φk∈[0,π/2]；令m=(mx-1)My+my，其中mx=1,2,…,Mx，my=1,2,…,My，則第m個(gè)陣元的坐標(biāo)為((mx-1)λ/2,(my-1)λ/2)。該陣列信號(hào)快拍數(shù)據(jù)模型可以表示為

x=As+v，

(1)

其中，x(l)為陣列接收信號(hào)的第l個(gè)采樣，稱為第l個(gè)快拍，l=1,2,…,L，L為快拍數(shù)；x=[x(1),x(2), … ,x(L)]、s=[s(1),s(2), …,s(L)]和v=[v(1),v(2), …,v(L)]分別為觀測(cè)信號(hào)、源信號(hào)和噪聲的快拍構(gòu)成的矩陣；A=[a(θ1,φ1),a(θ2,φ2), …,a(θK,φK)]，a(θk,φk)=ax(αk)?ay(βk)為第k個(gè)源信號(hào)方向的導(dǎo)向矢量，k=1,2,…,K，?為克羅內(nèi)克積，ax(αk)=[1,exp(jπαk), …, exp(jπ(Nx-1)αk)]T，ay(βk)=[1, exp(jπβk), …, exp(jπ(Ny-1)βk)]T，j為虛數(shù)單位，(·)T為轉(zhuǎn)置運(yùn)算，且有

αk=cosθksinφk,

(2a)

βk=sinθksinφk。

(2b)

假設(shè)各陣元的噪聲為獨(dú)立同分布的高斯白噪聲，其功率為σ2。若使用均勻矩形陣列的部分M(M>K)個(gè)陣元構(gòu)造的稀疏矩形陣列，則陣列信號(hào)可以表示為

xΩ=ΓΩx=ΓΩAs+vΩ,

(3)

其中，ΓΩ∈{0,1}M×M0為選擇矩陣，xΩ和vΩ分別為稀疏矩形陣列接收信號(hào)矩陣和噪聲矩陣。均勻矩形陣列可以視為稀疏矩形陣列的特殊情形，此時(shí)對(duì)應(yīng)的ΓΩ為M0維單位矩陣。圖2給出了一個(gè)稀疏矩形陣列的示例，此時(shí)Mx=6，My=5，M=10，M0=30，根據(jù)此稀疏矩形陣列與均勻矩形陣列陣元位置的關(guān)系可知，ΓΩ由M0維單位矩陣IM0的第1、2、5、6、8、15、19、21、27和30行構(gòu)成。

圖1 均勻矩形陣列的信號(hào)模型

圖2 稀疏矩形陣列與對(duì)應(yīng)的均勻矩形陣列關(guān)系圖

2 基于均勻矩形陣列或稀疏矩形陣列的二維無(wú)格波達(dá)方向估計(jì)方法

2.1 二階特普利茨矩陣表達(dá)

對(duì)于模型式(3)，協(xié)方差矩陣為

(4)

(5)

(6)

2.2 二階特普利茨矩陣重構(gòu)

(7)

根據(jù)舒爾補(bǔ)理論[23]，優(yōu)化問(wèn)題式(7)等價(jià)為

(8)

(9)

(10)

優(yōu)化問(wèn)題式(10)的目標(biāo)函數(shù)是非凸的，但可以使用優(yōu)化最小化(Majorization-Minimization, MM)算法[24]，將其分為Q次半定規(guī)劃問(wèn)題求解，其中第q次求解的優(yōu)化問(wèn)題為

(11)

(12)

2.3 二維波達(dá)方向估計(jì)

An-Ψn=An-,

(13)

(14)

(15a)

(15b)

2.4 對(duì)比分析

將筆者提出方法的計(jì)算復(fù)雜度與估計(jì)性能進(jìn)行對(duì)比分析。對(duì)比的方法為最近提出的FGML[13]、LRMR[14]、ANM[17]和GL-SPICE[19]4種二維無(wú)格波達(dá)方向估計(jì)方法。

先考慮算法的計(jì)算復(fù)雜度，由復(fù)數(shù)乘法次數(shù)的階數(shù)進(jìn)行表達(dá)。二維ESPRIT方法可以使用截?cái)嗟钠娈愔捣纸猓?jì)算復(fù)雜度為O(K2M0)[22]，5種方法均使用該方法進(jìn)行二維波達(dá)方向估計(jì)?？焖贌o(wú)格最大似然方法有閉式解，計(jì)算復(fù)雜度為O{(2Mx-1)3(2My-1)3}，在5種方法中最低；其他3種方法需要迭代求解，低秩矩陣重構(gòu)的計(jì)算復(fù)雜度為O{(2MxMy-2Mx-2My+1)3T}，原子范數(shù)最小化與GL-SPICE的計(jì)算復(fù)雜度為O{(K2+2MxMy-2Mx-2My+1)3T},筆者提出方法的計(jì)算復(fù)雜度為O{(K2+2MxMy-2Mx-2My+1)3TQ}，其中T為各優(yōu)化問(wèn)題的迭代次數(shù)，并假設(shè)其在各方法中相等。可以看出，筆者提出的方法由于需要求解半定規(guī)劃問(wèn)題Q次，在5種方法中具有最高的計(jì)算復(fù)雜度。

再考慮波達(dá)方向估計(jì)性能。原子范數(shù)最小化實(shí)際上是l1范數(shù)最小化的連續(xù)域?qū)崿F(xiàn)[17]；GL-SPICE是協(xié)方差矩陣的擬合方法[21]；低秩矩陣重構(gòu)是利用無(wú)噪聲條件下協(xié)方差矩陣低秩信息構(gòu)造的方法[14]；快速無(wú)格最大似然方法并未約束協(xié)方差矩陣的正定性[13]；筆者提出的方法使用的稀疏測(cè)度是由原子l0范數(shù)通過(guò)log-det函數(shù)平滑替代而來(lái)的，在參數(shù)τ趨近于0時(shí)，其趨近原子l0范數(shù)[18]。因此，相比其他方法，筆者提出的方法能夠提供更稀疏的解，從而可以獲得更好的波達(dá)方向估計(jì)性能。節(jié)3的仿真實(shí)驗(yàn)可以驗(yàn)證該結(jié)論。

3 仿真實(shí)驗(yàn)

源信號(hào)使用零均值且方差為1的獨(dú)立同分布的復(fù)高斯信號(hào)，每個(gè)陣元添加等功率的復(fù)高斯白噪聲。接收信號(hào)信噪比定義為rSNR=10 lg(K/σ2)。使用均方根誤差作為算法性能評(píng)價(jià)標(biāo)準(zhǔn)，方位角與俯仰角的總均方根誤差可表示為

(16)

3.1 均勻矩形陣列情形

使用4×4的16元均勻矩形陣列。除了筆者提出的方法，同時(shí)仿真了二維ESPRIT方法(記為2-D ESPRIT)、快速無(wú)格最大似然方法、低秩矩陣重構(gòu)、原子范數(shù)最小化、GL-SPICE和對(duì)應(yīng)的克拉美羅界(Cramér-Rao Bound, CRB)[8]進(jìn)行對(duì)比。

(1)快拍數(shù)的影響。設(shè)置源數(shù)目為2，真實(shí)方位角分別為-40.17°和80.23°，真實(shí)俯仰角分別為50.30°和60.25°。設(shè)置信噪比為3 dB，對(duì)不同的快拍數(shù)進(jìn)行獨(dú)立實(shí)驗(yàn)，結(jié)果如圖3(a)所示?？梢钥闯?，2-D ESPRIT、原子范數(shù)最小化、GL-SPICE與筆者提出的方法在所有快拍數(shù)條件下的均方根誤差均接近克拉美羅界，而快速無(wú)格最大似然方法與低秩矩陣重構(gòu)只在快拍數(shù)不小于70時(shí)接近克拉美羅界，在低快拍數(shù)時(shí)均方根誤差結(jié)果一般。

(2)信噪比的影響。同上設(shè)置兩個(gè)源信號(hào)，設(shè)置快拍數(shù)為200，對(duì)不同信噪比進(jìn)行獨(dú)立實(shí)驗(yàn)，結(jié)果如圖3(b)所示?？梢钥闯?，6種方法在-9 dB到12 dB的信噪比條件下的均方根誤差均接近克拉美羅界，性能比較接近。

圖3 均勻矩形陣列下的均方根誤差仿真實(shí)驗(yàn)結(jié)果

(3)相鄰源的影響。設(shè)置源數(shù)目為2，兩個(gè)源的俯仰角均為50°，第1個(gè)源的方位角為-1°到1°隨機(jī)產(chǎn)生，第2個(gè)源的方位角為第1個(gè)源的方位角加一個(gè)角度間隔，這里實(shí)驗(yàn)設(shè)置該角度間隔分別為2°，4°，…，16°，圖3(c)顯示了快拍數(shù)為200，信噪比為6 dB時(shí)的仿真實(shí)驗(yàn)結(jié)果。可以看出，原子范數(shù)最小化表現(xiàn)較差，其他方法的表現(xiàn)差別不大。

3.2 稀疏矩形陣列情形

使用圖2所示10元稀疏矩形陣列。由于2-D ESPRIT無(wú)法直接使用，該情形下仿真快速無(wú)格最大似然方法、低秩矩陣重構(gòu)、原子范數(shù)最小化、GL-SPICE、克拉美羅界與文中所述方法進(jìn)行對(duì)比。

(1)考慮快拍數(shù)的影響。設(shè)置源數(shù)目為2，真實(shí)方位角分別為-40.17°和80.23°，真實(shí)俯仰角分別為50.30°和60.25°。設(shè)置信噪比為3 dB，對(duì)不同快拍數(shù)進(jìn)行獨(dú)立實(shí)驗(yàn)，結(jié)果如圖4(a)所示。可以看出，與均勻矩形陣列情形不同，隨著快拍數(shù)的增加，快速無(wú)格最大似然方法的均方根誤差值不再明顯減小，這是因?yàn)樵撉樾蜗轮貥?gòu)的特普利茨矩陣并不具有正定性，導(dǎo)致波達(dá)方向估計(jì)性能下降；低秩矩陣重構(gòu)、原子范數(shù)最小化、GL-SPICE與文中方法的均方根誤差值隨著快拍數(shù)增加而減小，原子范數(shù)最小化與文中所述方法的均方根誤差更接近克拉美羅界。

(2)不同信噪比的影響。設(shè)置快拍數(shù)為100，對(duì)不同信噪比進(jìn)行獨(dú)立實(shí)驗(yàn)，結(jié)果如圖4(b)所示?？梢钥闯?，與不同快拍數(shù)情形相似，快速無(wú)格最大似然方法的均方根誤差值隨著信噪比增加不再明顯減小，文中方法在大部分信噪比下均具有最小的均方根誤差，更接近克拉美羅界。

(3)5個(gè)源信號(hào)的情形。其真實(shí)方位角分別為-100°、-50°、0°、80°和120°，真實(shí)俯仰角分別為50.34°、65.48°、55.95°、63.66°和60.25°。分別設(shè)置信噪比為9 dB和快拍數(shù)為200，對(duì)不同快拍數(shù)和信噪比進(jìn)行獨(dú)立實(shí)驗(yàn)，結(jié)果分別如圖4(c)和(d)所示?？梢钥闯觯闹兴龇椒ㄔ诖罂炫臄?shù)或高信噪比下的性能明顯更好，能接近克拉美羅界。

(4)相鄰源的影響。源信號(hào)與第3.1節(jié)相鄰源實(shí)驗(yàn)一致。圖4(e)顯示了快拍數(shù)為200，信噪比為6 dB時(shí)，5種方法在不同方位角間隔下的仿真實(shí)驗(yàn)結(jié)果?？梢钥闯觯S著角度間隔的增大，5種方法的均方根誤差均有所減?。豢焖贌o(wú)格最大似然方法需要方位角間隔不小于14°時(shí)，才有可以接受的均方根誤差值，其他方法表現(xiàn)遠(yuǎn)優(yōu)于快速無(wú)格最大似然方法；文中所述方法在角度間隔為4°時(shí)仍有很低的均方根誤差，具有最好的相鄰源波達(dá)方向估計(jì)性能。

(5)各算法的平均運(yùn)行時(shí)間。使用上述5個(gè)源信號(hào)的設(shè)置，令快拍數(shù)為200，信噪比為9 dB，5種方法的200次平均運(yùn)行時(shí)間如表1所示(計(jì)算機(jī)的CPU為英特爾i5-8500，內(nèi)存的頻率和容量分別為2 666 MHz和16 GB)。可以看出，快速無(wú)格最大似然方法的平均運(yùn)行時(shí)間最短，文中方法運(yùn)行時(shí)間最長(zhǎng)，驗(yàn)證了節(jié)2.4的計(jì)算復(fù)雜度分析。

圖4 稀疏矩形陣列下的均方根誤差仿真實(shí)驗(yàn)結(jié)果

算法FGMLLRMRANMGL-SPICE文中方法平均運(yùn)行時(shí)間/s0.020 20.161 50.212 50.192 10.415 9

4 結(jié)束語(yǔ)

筆者提出一種基于均勻矩形陣列或稀疏矩形陣列的二維無(wú)格波達(dá)方向估計(jì)方法。該方法使用樣本協(xié)方差矩陣和log-det稀疏測(cè)度對(duì)相關(guān)的二階特普利茨結(jié)構(gòu)矩陣進(jìn)行重構(gòu)；然后通過(guò)二維ESPRIT方法進(jìn)行二維波達(dá)方向估計(jì)。相比快速無(wú)格最大似然方法、均方根誤差、原子范數(shù)最小化與GL-SPICE，文中方法能提供更稀疏的解，具有更好的波達(dá)方向估計(jì)性能，但計(jì)算復(fù)雜度相對(duì)較高。仿真實(shí)驗(yàn)中，筆者提出的方法具有在均勻矩形陣列和稀疏矩形陣列下的各情形中均具有很低的均方根誤差，證明筆者提出方法良好的波達(dá)方向估計(jì)性能。在未來(lái)的研究中，可以考慮將該方法中使用的優(yōu)化最小化算法替換為其他復(fù)雜度更低的優(yōu)化算法，從而節(jié)省算法運(yùn)行時(shí)間。