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

彭加強，鄭桂妹

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

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

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

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

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

1 信號模型

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

(1)

(2)

即求出αk、βk便可根據(jù)式(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次快拍數(shù)據(jù)及所有快拍數(shù)據(jù)可以表示為:

(7)

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

(8)

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

圖1 URA信號模型

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

2.1 基于矢量的原子范數(shù)最小化算法

根據(jù)文獻[7]，接收信號數(shù)據(jù)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)構(gòu)造雙重Toeplitz矩陣T2D(P)，其中-N

(11)

(12)

‖X(l)‖AV=

(13)

‖X(l)‖AV=

(14)

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

2.2 解耦原子范數(shù)最小化算法

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

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

(15)

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

(16)

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

定理對于一個N×M的數(shù)據(jù)矩陣：

(17)

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

(18)

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

‖X(l)‖AM=

(19)

式中：T(ux)和T(uy)表示一重Toeplitz矩陣，分別使用ux和uy作為其第1行進行構(gòu)造。在有噪情況下，式(19)將變?yōu)橐韵耂DP求解模型：

(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多快拍接收數(shù)據(jù)[11]。

(23)

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

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

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

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

(24)

(25)

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

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

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

3.2 復雜度計算

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

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

4 數(shù)值仿真

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

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網(wǎng)格間隔取2°。從圖2中可以看出DANM和本文算法可以準確地估計出3個信號源的2D DOA，而OMP算法遭受了嚴重的網(wǎng)格失配問題，這一結(jié)果與信號源角度和網(wǎng)格間隔的設定相對應?？梢?，在[-40°,40°]×[-40°,40°]的空域范圍內(nèi)，只有x維中的2°落在網(wǎng)格上，但其y維對應的角度5°又不在網(wǎng)格上，因此導致OMP算法整體估計失效。而U-ESPRIT算法作為子空間類經(jīng)典超分辨算法在既設條件下的估計精度相比于本文算法和DANM明顯需要更大的快拍和更高的SNR。

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

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

圖3 RMSE性能比較

4.2 運行時間比較

假設K=2,在圖4中N=M=8∶22，DANM采用單快拍，本文算法采用200次快拍，2個獨立ANM求解的次優(yōu)DANM[11]采用200次快拍；在圖5中N=M=8∶20，基于對偶的2D ANM和本文算法均采用5次快拍。為了更公平地進行比較，運行時間均以cvx_cputime進行統(tǒng)計，該時間也是各算法所需運行時間的主要部分。從圖4來看，雖然本文算法采用200次快拍，但在N=M=22時所需運行時間相比單快拍的DANM僅增加了0.812 s，而次優(yōu)DANM的運行時間遠大于本文算法，近似為DANM的2倍。在圖5中可以清晰地看到，本文算法相比基于對偶的2D ANM在運行時間上顯著降低，在2個維度的維數(shù)都增加到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次蒙特卡羅仿真。仿真中采用數(shù)據(jù)壓縮的形式構(gòu)造稀疏信號，數(shù)據(jù)壓縮比定義為：

(27)

(28)

(29)

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

圖6 稀疏恢復性能比較

4.4 快拍數(shù)對本文算法的影響

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

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

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

4.5 信源數(shù)估計能力分析

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

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

5 結(jié)語

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