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

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

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

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

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

1 均勻矩形陣列和稀疏矩形陣列的接收信號模型

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

x=As+v，

(1)

其中，x(l)為陣列接收信號的第l個采樣，稱為第l個快拍，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)]分別為觀測信號、源信號和噪聲的快拍構成的矩陣；A=[a(θ1,φ1),a(θ2,φ2), …,a(θK,φK)]，a(θk,φk)=ax(αk)?ay(βk)為第k個源信號方向的導向矢量，k=1,2,…,K，?為克羅內克積，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為轉置運算，且有

αk=cosθksinφk,

(2a)

βk=sinθksinφk。

(2b)

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

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

(3)

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

圖1 均勻矩形陣列的信號模型

圖2 稀疏矩形陣列與對應的均勻矩形陣列關系圖

2 基于均勻矩形陣列或稀疏矩形陣列的二維無格波達方向估計方法

2.1 二階特普利茨矩陣表達

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

(4)

(5)

(6)

2.2 二階特普利茨矩陣重構

(7)

根據舒爾補理論[23]，優(yōu)化問題式(7)等價為

(8)

(9)

(10)

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

(11)

(12)

2.3 二維波達方向估計

An-Ψn=An-,

(13)

(14)

(15a)

(15b)

2.4 對比分析

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

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

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

3 仿真實驗

源信號使用零均值且方差為1的獨立同分布的復高斯信號，每個陣元添加等功率的復高斯白噪聲。接收信號信噪比定義為rSNR=10 lg(K/σ2)。使用均方根誤差作為算法性能評價標準，方位角與俯仰角的總均方根誤差可表示為

(16)

3.1 均勻矩形陣列情形

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

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

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

圖3 均勻矩形陣列下的均方根誤差仿真實驗結果

(3)相鄰源的影響。設置源數(shù)目為2，兩個源的俯仰角均為50°，第1個源的方位角為-1°到1°隨機產生，第2個源的方位角為第1個源的方位角加一個角度間隔，這里實驗設置該角度間隔分別為2°，4°，…，16°，圖3(c)顯示了快拍數(shù)為200，信噪比為6 dB時的仿真實驗結果?？梢钥闯觯臃稊?shù)最小化表現(xiàn)較差，其他方法的表現(xiàn)差別不大。

3.2 稀疏矩形陣列情形

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

(1)考慮快拍數(shù)的影響。設置源數(shù)目為2，真實方位角分別為-40.17°和80.23°，真實俯仰角分別為50.30°和60.25°。設置信噪比為3 dB，對不同快拍數(shù)進行獨立實驗，結果如圖4(a)所示?？梢钥闯觯c均勻矩形陣列情形不同，隨著快拍數(shù)的增加，快速無格最大似然方法的均方根誤差值不再明顯減小，這是因為該情形下重構的特普利茨矩陣并不具有正定性，導致波達方向估計性能下降；低秩矩陣重構、原子范數(shù)最小化、GL-SPICE與文中方法的均方根誤差值隨著快拍數(shù)增加而減小，原子范數(shù)最小化與文中所述方法的均方根誤差更接近克拉美羅界。

(2)不同信噪比的影響。設置快拍數(shù)為100，對不同信噪比進行獨立實驗，結果如圖4(b)所示。可以看出，與不同快拍數(shù)情形相似，快速無格最大似然方法的均方根誤差值隨著信噪比增加不再明顯減小，文中方法在大部分信噪比下均具有最小的均方根誤差，更接近克拉美羅界。

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

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

(5)各算法的平均運行時間。使用上述5個源信號的設置，令快拍數(shù)為200，信噪比為9 dB，5種方法的200次平均運行時間如表1所示(計算機的CPU為英特爾i5-8500，內存的頻率和容量分別為2 666 MHz和16 GB)?？梢钥闯觯焖贌o格最大似然方法的平均運行時間最短，文中方法運行時間最長，驗證了節(jié)2.4的計算復雜度分析。

圖4 稀疏矩形陣列下的均方根誤差仿真實驗結果

算法FGMLLRMRANMGL-SPICE文中方法平均運行時間/s0.020 20.161 50.212 50.192 10.415 9

4 結束語

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