基于解耦原子范數最小化的二維DOA估計

彭加強，鄭桂妹

(空軍工程大學防空反導學院，西安，710051)

2D DOA估計技術應用十分廣泛，比如MIMO雷達的波離角和波達角聯合估計、二維陣列的方位角和俯仰角聯合估計等都可以歸為2D DOA估計問題。得益于1D DOA估計問題的成功應用，子空間類超分辨算法被順利地推廣到2D DOA估計問題，比如2D酉信號參數旋轉不變估計技術(unitary estimation of signal parameters via rotational invariance techniques,U-ESPRIT)[1]，2D多重信號分類算法(multiple signal classification,MUSIC)[2]等。該類算法在已知信源數目、獲得多快拍采樣數據及信號源非相干的情況下可以取得較好的估計性能，但其中任一條件不滿足就會導致其估計性能下降，甚至失效。

近年來，基于壓縮感知技術框架的有效替代算法被引入到2D DOA估計領域，取得了豐碩的研究成果。傳統的壓縮感知算法將信號源可能出現的空域按維度劃分為有限個網格點，在真實目標方位準確地落在既設網格上時展現出了優良的估計性能，能夠適用于單快拍、信源相干、數據缺失等復雜場景。比如正交匹配追蹤(orthogonal matching persuit,OMP)[3],稀疏貝葉斯學習(sparse bayesian learning,SBL)[4]等算法。反之，如果真實信號源沒有落在既設網格上則會造成網格失配的問題，估計性能也會隨之下降，甚至失效。另外，傳統壓縮感知算法還必須滿足成對等距特性(pairwise isometry property,PIP)[5]及高密度網格劃分。為克服以上問題，一種新的基于原子范數理論和范德蒙德分解定理的無網格連續壓縮感知技術被提出，稱為原子范數最小化(atomic norm minimization,ANM)[6]。ANM通過拓普利茲(Toeplitz)矩陣將觀測數據中陣列流行矩陣的范德蒙德結構投射到對應的半定規劃(semi-definite programming,SDP)模型，并通過對SDP優化模型的求解獲得恢復的信號實現超分辨。相比于傳統壓縮感知算法，其無需對空域進行網格劃分，有效避免了網格失配的問題以及PIP限制。但根據Caratheodory的理論，Toeplitz矩陣的范德蒙德分解在高維空間中不成立，因此1D ANM不能直接擴展到2D DOA估計。幸運的是，通過接收數據的矢量化操作及一種包含2個維度范德蒙德結構的雙重Toeplitz矩陣的構造，Chi等人成功解決了該問題，并率先將矢量化ANM(VANM)應用于2D DOA估計[7]，但對接收測量數據的矢量化操作及雙重Toeplitz矩陣的構造使得該算法計算代價巨大，不能應用于實際中的場景。基于對偶的2D ANM[8]從VANM的對偶問題出發進行求解，但并沒有減輕VANM的高計算量。

為減輕繁重的計算負擔，Tian等人提出一種新的解耦原子范數最小化算法(decoupled atomic norm minimization,DANM)[9-11]。DANM將VANM中的矢量原子集替換為矩陣原子集，并推導出相應的SDP模型，該模型天然地將VANM中的雙重Toeplitz矩陣解耦為2個分別包含一維范德蒙德結構的Toeplitz矩陣，從而將2D DOA估計問題轉換為2個1D DOA估計問題，該算法在保持ANM類算法優良估計性能的同時顯著降低了計算復雜度，相較于VANM降低了幾個數量級。但原始的基于均勻矩形陣列(uniform rectangle array,URA)的DANM算法目前只能工作于單快拍，對于多快拍DANM的研究只見于特殊的陣列結構，如文獻[12～13]基于L型陣列2個子陣的互協方差矩陣對DANM進行了相應地改進，使其能夠適用于多快拍的場景，文獻[14]針對互質陣列對DANM進行了相應地改進。同時，其他對DANM的研究也局限于單快拍，如文獻[15]利用交替方向乘子法(alternating direction method of multipliers,ADMM)研究了單快拍DANM的快速算法，文獻[16]利用DANM進行MIMO雷達角度和距離的單快拍聯合估計。

為使DANM適用于多快拍，本文提出一種改進的解耦原子范數最小化方法。

1 信號模型

考慮空間遠場K個窄帶信號作用于一個N×M的均勻矩形陣列，陣元間距為半波長，如圖1所示，其L次快拍采樣數據表示為sk,l∈CK×L。鑒于第k個入射信號的俯仰角φk、方位角θk同其與x、y軸之間的夾角αk、βk具有如下關系：

(1)

(2)

即求出αk、βk便可根據式(1)～(2)得到θk、φk，因此本文使用αk、βk進行信號建模分析。則x、y維的陣列導向矢量、流行矩陣分別為:

(3)

(4)

Ax=[ax(α1),ax(α2),…,ax(αK)]

(5)

Ay=[ay(β1),ax(β2),…,ay(βK)]

(6)

第l次快拍數據及所有快拍數據可以表示為:

(7)

X=[X(1),X(2),…,X(L)]
Y=[Y(1),Y(2),…,Y(L)]

(8)

2D DOA估計就是要將所有的αk和βk從觀測數據X或Y中恢復出來。本文中我們主要以無噪數據的形式進行建模，但也會對有噪的情況進行討論。

圖1 URA信號模型

2 2D DOA估計中的原子范數最小化算法

2.1 基于矢量的原子范數最小化算法

根據文獻[7]，接收信號數據X(l)可以被矢量化為以下表達式:

(9)

相應地原子集AV可以表示為:

AV={ay(βk)?ax(αk),αk,βk∈[-90,90]}=

{α(γ),γ∈[-90,90]×[-90,90]}

(10)

利用矩陣變量P=[ul1,l2]∈C(2N-1)×(2M-1)構造雙重Toeplitz矩陣T2D(P)，其中-N

(11)

(12)

‖X(l)‖AV=

(13)

‖X(l)‖AV=

(14)

從式(11)～(14)中可以得出VANM優化模型中半定約束矩陣的維數為(NM+1)×(NM+1)，這就直接導致了VANM的高計算量，當維數N和M較大時甚至到了不可接受的地步。

2.2 解耦原子范數最小化算法

根據文獻[9]，由(7)可以得到另一種矩陣形式的原子集合:

{Aγ,γ∈[0,2π]×[0,2π]}

(15)

其中的每一個原子為秩1矩陣，在單快拍情況下對應的原子范數為:

(16)

為從式(16)中求出各維度的DOA，引入以下定理。

定理對于一個N×M的數據矩陣：

(17)

定義最小角度間隔為Δmin,x=mini≠j|sinαx,i-sinαx,j|，Δmin,y=mini≠j|sinβy,i-sinβy,j|，如果它們滿足：

(18)

則式(17)為式(16)的最優解。進而，式(16)可通過式(19)進行有效求解。

‖X(l)‖AM=

(19)

式中：T(ux)和T(uy)表示一重Toeplitz矩陣，分別使用ux和uy作為其第1行進行構造。在有噪情況下，式(19)將變為以下SDP求解模型：

(20)

在得到T(ux)和T(uy)后，便可通過以下分解得到x維和y維的DOA，

(21)

式中：Dx、Dy為對角矩陣，在得到各維度DOA后再通過配對程序得到最終的2D DOA。

3 本文算法

3.1 算法步驟

為使DANM適用于多快拍，本文對其進行如下改進。

步驟1將式(19)中的約束條件由1個改為2個對等條件，得到如下SDP求解模型：

(22)

式中：Z1∈CN×N、Z2∈CM×M表示埃爾米特Toeplitz矩陣。而X(l)、Y(l)在如下表述下可分別視為x維和y維的1D多快拍接收數據[11]。

(23)

雖然文獻[11]中提到一種次優的完全分解的解耦方法，即對兩個維度分別使用1D ANM進行計算，但本文算法并非該次優方法的復制。原因如下：

1)本文算法只需1個優化求解模型，而文獻[11]中方法需要2個。

2)文獻[11]分析到該完全分離的次優解耦方法忽視了2個維度的聯合信息，而本文算法將2個維度的信息約束于1個目標函數，保留了2個維度的聯合信息。

步驟2為使多快拍數據能夠代入式(19)中運算求解而不增加模型的維度，本文參考文獻[17]中V.C大快拍降維部分內容，分別使用X、Y的協方差矩陣Rx、Ry替換式(19)中的單快拍數據X(l)、Y(l),得到最終的多快拍DANM SDP求解模型如下。

(24)

(25)

在本文中加權因子λ的取值與文獻[6]的3.2節選擇正則化參數中陣元個數大于3時的取值一致。

步驟4通過如下配對程序進行配對得到2D DOA。

2)通過Pi,j中K個較大值的下標索引得到2D DOA。

3.2 復雜度計算

給出本文算法、DANM和基于對偶的2D ANM算法的計算復雜度，以便為后續數值仿真環節提供理論支撐。

又根據文獻[8]，基于對偶的2D ANM的計算復雜度為O((NM+L)3.5log(1/ε)+UVL(NM+1))。其中L表示快拍數，U、V分別表示2個維度DOA搜索的次數。

4 數值仿真

本文數值仿真部分在一臺Intel(R) Core(TM) i5-5200U@2.2 GHz處理器上完成，SDP模型使用CVX工具箱進行求解，統計誤差分析采用均方根誤差(root mean squared error,RMSE)及其平方(均方誤差)。其中，N表示蒙特卡羅仿真次數。

RMSE=

(26)

4.1 角度估計性能分析

假設N=M=10,K=3,x、y維入射角度分別為αx=[-35°,2°,31°],βy=[-29°,5°,37°]。DANM、OMP及本文算法均采用單快拍，SNR取0 dB；U-ESPRIT采用20次快拍，SNR取10 dB，OMP網格間隔取2°。從圖2中可以看出DANM和本文算法可以準確地估計出3個信號源的2D DOA，而OMP算法遭受了嚴重的網格失配問題，這一結果與信號源角度和網格間隔的設定相對應。可見，在[-40°,40°]×[-40°,40°]的空域范圍內，只有x維中的2°落在網格上，但其y維對應的角度5°又不在網格上，因此導致OMP算法整體估計失效。而U-ESPRIT算法作為子空間類經典超分辨算法在既設條件下的估計精度相比于本文算法和DANM明顯需要更大的快拍和更高的SNR。

圖3中，對U-ESPRIT、DANM及本文算法在SNR取0∶5∶20 dB時對3個信源的2D DOA估計RMSE進行了50次蒙特卡羅仿真。其中U-ESPRIT采用200次快拍，本文算法采用5次快拍，DANM采用單快拍。從仿真結果來看，圖3進一步體現了ANM類算法相比子空間類算法在估計精度上的優勢，同時，基于多快拍的DANM相比原DANM在估計精度上取得了較大的提升。

圖2 不同算法的2D DOA估計

圖3 RMSE性能比較

4.2 運行時間比較

假設K=2,在圖4中N=M=8∶22，DANM采用單快拍，本文算法采用200次快拍，2個獨立ANM求解的次優DANM[11]采用200次快拍；在圖5中N=M=8∶20，基于對偶的2D ANM和本文算法均采用5次快拍。為了更公平地進行比較，運行時間均以cvx_cputime進行統計，該時間也是各算法所需運行時間的主要部分。從圖4來看，雖然本文算法采用200次快拍，但在N=M=22時所需運行時間相比單快拍的DANM僅增加了0.812 s，而次優DANM的運行時間遠大于本文算法，近似為DANM的2倍。在圖5中可以清晰地看到，本文算法相比基于對偶的2D ANM在運行時間上顯著降低，在2個維度的維數都增加到20時，基于對偶的2D ANM需要1 749.343 75 s，而本文算法僅需11.406 3 s。

圖4 解耦ANM算法運行時間比較

圖5 本文算法與基于對偶的2D ANM在多快拍下的運行時間比較

4.3 無噪情況下的稀疏恢復能力比較

假設N=M=10,K=4，2個維度的DOA分別為αx=[-5°,8°,17°,31°]，βy=[0°,10°,23°,37°]。DANM采用單快拍，本文算法采用10個快拍，針對每個稀疏信號分別采用50次蒙特卡羅仿真。仿真中采用數據壓縮的形式構造稀疏信號，數據壓縮比定義為：

(27)

(28)

(29)

從圖6中可以看到，DANM在壓縮比為36%時才能精確恢復所有數據，而本文算法在壓縮比為16%時便可以實現，體現出更強的稀疏恢復能力。

圖6 稀疏恢復性能比較

4.4 快拍數對本文算法的影響

鑒于ANM類算法天然地適用于小快拍場景，同時考慮到實際應用場景中計算效率的問題，在本次實驗中，假設快拍數L=5∶5∶100,SNR取20 dB，2個維度的DOA分別為αx=[3°,17°,31°]，βy=[10°,23°,37°]，陣元數目與4.3節中相同，采用50次蒙特卡羅仿真。

圖7 快拍數對本文算法估計性能的影響

從圖7中可以看出，隨著快拍數的增長，2個維度的估計誤差總體呈現下降趨勢，快拍數10和40為曲線的2個拐點。在快拍數達到10次以后，估計誤差基本處于0.01～0.001之間，該結果表明本文算法同樣適用于小快拍，保持了ANM類算法的優勢；在快拍數達到40次以后，估計誤差又出現了明顯的區域性下降。因此，在實際應用中可根據陣元數、信源數及實時性要求等條件，選擇大于等于10次的快拍數進行數據處理。

4.5 信源數估計能力分析

本次實驗采用均方誤差作為評判依據，并設定RMSE<0.01視為能夠正確估計，2個維度的DOA按照αx=[0°,10°,…]，βy=[10°,20°,…]進行設定，SNR取20 dB，陣元數目與4.3節中相同，采用50次蒙特卡羅仿真。根據文獻[11]，理論上本文算法能夠正確估計出9個信源，但限于壓縮感知類算法對稀疏性的要求，從圖8可見，本文算法只能準確估計出5個。

圖8 本文算法在不同信源數下的估計性能

5 結語

本文詳細闡述了將DANM方法由單快拍推廣至多快拍的一種改進方法和步驟，并對改進方法的角度估計性能、稀疏恢復能力、計算復雜度、對快拍數的依賴性等方面進行了對比分析。數值仿真結果表明，本文算法在保留DANM高效運算能力的同時，提高了有噪信號的估計精度和無噪稀疏信號的恢復能力。