基于解耦原子范数最小化的二维DOA估计

2021-04-08彭加强郑桂妹

空军工程大学学报 2021年1期

彭加强，郑桂妹

(空军工程大学防空反导学院，西安，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个。