0 引 言
方向到达角（Direction of Arrival, DOA）估计是阵列信号处理领域中研究的重要问题之一，旨在通过传感器阵列接收的空域信号实现空间辐射源方位的精准定位，在雷达、声呐、无线通信等领域具有重要应用［1］。其理论发展历经数十年革新，逐步形成了从经典谱估计到高分辨率子空间方法，再到现代稀疏重构技术的完整体系。
DOA估计目的是通过传感器阵列接收信号，精准定位空间辐射体的方位角。早期DOA估计主要基于波束形成技术，通过空间扫描计算不同方向的输出功率，但受限于瑞利准则，分辨率难以突破物理孔径的限制。为优化波束权重，Capon提出的最小方差无失真响应（MVDR）算法［2］通过对协方差矩阵逆矩阵加权，在抑制干扰的同时提升分辨能力，是经典波束形成的重要改进。然而，传统方法在相干信号或者低信噪比场景下的性能急剧下降。为了突破瑞利限，20世纪80年代，以多重信号分类（MUSIC）［3］和旋转不变子空间（ESPRIT）［4］为代表的子空间方法应运而生；MUSIC算法通过分解协方差矩阵，利用噪声子空间与噪声子空间的正交性构造空间谱函数，借助峰值搜索实现超分辨估计；ESPRIT算法则基于传感器阵列的平移不变结构，直接通过信号子空间的特征值分解求解DOA，避免了谱搜索的高计算量，提高角度分辨率。为进一步提升估计精度与鲁棒性，基于MVDR方法下提出了一种最大似然方法（ML）的子空间拟合算法［5］，目的是找到使接收信号的概率密度函数最大的方位角，随之先后提出了基于确定性假设的最大似然方法（DML）和考虑统计随机性的最大似然方法（SML）。然而，此类方法涉及高维非凸优化，计算复杂度限制了其工程适用性。为了减少计算量，提高计算效率，在角度估计中引进了快速傅里叶变换（FFT）算法，文献［6］介绍了一种基于FFT算法下的稀疏阵列的空间频率的DOA估计方法，可以同时对不同的入射角度进行角度估计，极大地降低了测向的计算量，易于工程实现。
21世纪初，压缩感知（CS）理论的突破为DOA估计开辟了新路径。文献［7］基于ℓ1范数的DOA估计方法忽略弱信号的问题，提出一种结合ℓ1范数与ℓ2范数平方的组合惩罚项模型，通过最大化-最小化算法求解，提高了角度估计精度。然而上述方法需要通过离散化角度网格构建稀疏模型，将DOA估计转化为凸优化或贪婪迭代，存在网格失配的问题，会导致性能显著退化，制约其在连续角度空间的应用。针对此瓶颈，原子范数最小化（ANM）［8］通过构建连续字典的稀疏表示框架，将无限维优化问题转化成为有限维半定规划（SDP），实现了无网格的超分辨估计，成为近年来的研究热点。文献［9］在原子范数最小化的ℓ0原子范数问题中，提出了基于Circle混沌映射无网格DOA估计算法，降低均方根误差。文献［10］在原子范数最小化的ℓ1原子范数问题中，提出了利用交替方向乘子法（ADMM）的快速迭代求解范数问题的方法，提高了估计精度，降低了计算量。但是用ADMM算法解决原子范数最小化问题需要先将原子范数转化为半定规划问题，涉及到大规模矩阵运算和求逆计算，导致计算量和内存消耗过大，抗噪声性能较差。
为此，本文在原子范数最小化的交替方向乘子法（ANM-ADMM）基础上，引入了最大最小化（MM）优化算法，提出了一种ANM-ADMM-MM算法来估计稀疏阵列下的DOA，相比于传统的ANM-ADMM算法，本文算法提高了计算效率和抗噪声性能，降低估计误差。
1 稀疏阵列信号模型
目前，均匀阵列下的DOA估计算法已趋于成熟，但该算法因依赖大量阵元，并且阵列单元需要布局严谨，导致硬件成本增加。因此少阵元的稀疏阵列的DOA估计技术一直是学术界与工程领域的研究热点［11］。若在均匀阵列基础上移除部分阵元单元，将对应通道数据置为0，便可构建出稀疏阵列，具体形式如图1所示。
图1 稀疏阵列示意图
设X∈M×L由L个相同的角度分量的复信号组成的无噪声的数据矩阵。X的第(m，l)个元素，记为xml，其形式为
式中，(m，l)∈[M]×[L]，M表示阵元个数，L表示快拍数，K表示角度分量的个数，θ∈[-90°，90°]表示角度分量，Skl表示第k个正弦分量在第l个快拍数的复振幅。为了简化式（1）的表示，假设以采样周期为单位，X的每一列是一个由M个均匀样本组成的测量向量。那么，信号模型式（1）可以用矩阵-向量的形式重写为
式中，A(θ)=[a(θ1),…,a(θK)]，其中，a(θ)=[1,ej2πdsinθ/λ,…,e-j2π(M-1)dsinθ/λ]T，因阵元1，2，6，7，9，13和14在a(θ)设置为0，则a(θ) = [0,0,…,e-j2π(M-1)dsinθ/λ]T，skl=[s1,s2,…,sk]T为矩阵S={skl}的第k行。
为了减少硬件的复杂度，采用1位ADC来分别对接收信号的实部和虚部进行采样和量化。首先，定义复值参考阙值矩阵H是从一组离散的8个值随机生成的，这8个值均匀分布在[-hmax,hmax]值域上，hmax的定义:
然后通过量化过程将有噪声的接收信号的样本与定义的复值参考阙值矩阵H进行差值比较计算，得到有符号的观测矩阵Y∈{yml}，Y的每个元素yml的表达式为
式中，eml和hml分别为未知的加性噪声和已知的阙值矩阵H的元素。signc(·)定义为
式中，signc(·)定义为
所以，Y={yml}∈{1+j,1-j,-1+j,-1-j}作为本文算法的输入观测矩阵。在许多雷达侦测和通信应用中，Y的角度分量K的数量远远小于阵元数M，基于这一特性，通过合理减少阵元数量不仅能够降低硬件成本，还能进一步优化DOA估计的性能，实现资源的有效利用和精度提升。本文的目标是从观测矩阵Y中通过改进的原子范数恢复信号，从而估计角度值θ。
2 算法描述
2.1 原子范数最小化
原子范数最小化的核心是通过构建合适的目标函数与约束条件，从而实现对原始信号的有效重构和分析。信号矩阵X的原子集合XΑ为
式中，ϕ∈C1×L，并且‖ϕ‖2=1。
信号矩阵X可以表示在原子集合上的最稀疏分解代价［12-13］为
通过凸松弛，可以转换为SDP问题:
其中，u=[u1,…, uM]T；T(u)为M×M的Toeplitz矩阵，定义为
根据Vandermonde原理分解［14］，任何半正定规划（SDP）公式的Toeplitz矩阵都可以重写为T(u)=
，因此由 Toeplitz矩阵可估计出信号的角度值θ。
2.2 基于ANM-ADMM-MM的DOA估计
MM是一种迭代优化算法，用于解决难以直接优化的目标函数问题。其核心思想是通过构造一个易于优化的替代函数来逐步逼近目标函数的局部最优解［15］。在每次迭代中，MM优化算法不是直接对原始目标函数进行求解，而是通过优化一个更简便的辅助函数来间接地求得最优解。
为了降低计算量，提高计算效率，在ADMM的基础上引入MM优化算法，在复杂问题中找到一个更简单的替代函数，称这个算法为ANM-ADMM-MM算法。加入MM优化算法思想如下，首先，观测矩阵Y的负对数似然函数［16］为
式中，η为噪声因素，Φ(x)为标准正态分布的累积函数（cdf），
。
由
关系可知，负对数似然函数
可以被替换为
，然后添加原子范数正则式（8）构建复合优化目标:
其中，前面一项是数据拟合的负对数似然函数，后面一项是矩阵
的原子范数项，μ是正则化参数，用于平衡两个分量之间的偏差和方差权衡。μ过大，原子范数项主导优化过程，模型会过度正则化，恢复出的信号
原子数量过少，幅度被压缩，无法很好地拟合数据中的真实结构，导致估计出的角度数量少于真实值；μ过小，负对数似然函数主导优化过程，模型会欠正则化，容易拟合数据中的噪声和随机波动，出现虚假的、强度弱的角度估计（假峰）。选择理想的μ，可以在原子范数项和负对数似然函数取得最佳的平衡，使角度估计的均方误差最小化。正则化参数μ的理论经验公式［17］:
其中，α的取值范围是［1,10］。当SNR=0 dB，α=1时，μ=4.816 4；当SNR=15 dB，α=10时，μ=8.565 0，所以μ的取值选择［4.816 4，8.565 0］可以平衡原子范数项和负对数似然函数。故后续的仿真实验，正则化参数选取μ=5。
角度估计值θ可以由式（12）估计出来。将式（9）代入式（12），替换原子范数函数得到目标函数问题为
由式（14）可以看出，这是多个非线性复数项的SDP问题，难以用凸优化工具箱CVX解决。此时，需要用ADMM方法迭代来求解问题（14），恢复信号
，从中得到角度估计值θ。首先引入辅助变量矩阵Z∈C(M+L)×(M+L)和辅助参数τ，将式（14）改为
式中tr(·)表示矩阵的迹。式（15）中的
加入MM优化算法，首先定义c(x)=lnΦ(x)，则式（11）可以表示为
对所有的x,u∈R，x和u的关系［18］:
并且，有且仅当x=u等号成立。令:
将式（18）带入式（17）可以得到
同理对虚部项可以得到类似式（19），对目标函数
求和:
式中，
。对式（20）的二次部分进行配方得到
将式（21）代入式（20）得到
由{yml}∈{1+j,1-j,-1+j,-1-j}，可以得到
，所以式（22）的二次项可以表示为
式中，
。令Qi={ymlqml}，则式（22）可以写为
根据式（15）和式（24）得到目标函数优化:
式（25）的增广拉格朗日函数［19］为
其中，Λ和λ是拉格朗日函数的两个拉格朗日乘数；ρ>0是ADMM的惩罚参数，是ADMM算法中平衡原始可行性和对偶可行性的关键参数。ρ过大，强调原始可行性，原始残差会快速收敛，但对偶残差收敛缓慢，导致解的质量不高；ρ过小，对偶残差收敛相对较快，但原始残差收敛缓慢，导致算法运行时间过长。选择理想的ρ，可以在原始残差和对偶残差的收敛速度之间找到一个平衡点，使得总收敛速度最快。在本文稀疏优化问题中，正则化参数选择μ=5较大，此时惩罚参数选择较小可使原始可行性和对偶可行性平衡。故后续的仿真实验，惩罚参数选取ρ=0.3。
表1为求解式（15）的闭式解，其中tr+n(·)表示矩阵第几个对角线之和，+是主对角线上方的子对角线，PS+(·)表示在正半定矩阵下的投影。Z和Λ辅助变量的定义为
表1 ANM-ADMM-MM算法闭式解
输入:观测矩阵Y，复值参考阙值矩阵H，正则化参数μ和ADMM惩罚参数ρ。参数初始化:x~1=0M×L，η1=τ1=0，λ1=1，i=1，Z1=Λ1=0(M+L)×(M+L)，I=EL×L。迭代更新:1）更新Q:由式（16）到式（24）的推导过程更新；2）更新η:ηi+1=[trρ(tHr(HH)(H-HρQ)i++ρρ(Zρxi++0Λ.5xi))]R+λi(ρ+0.5)+τiρ(ρ+0.5)ρtr(HHH)+ρ(ρ+0.5)；3）更新~X:~Xi+1=(0.5ηi+1H+0.5Qi+ρZxi+Λxi)/(0.5+ρ)；4）更新W:Wi+1=(ρZ1i+Λi1-2μMI)/ρ；5）更新u:ui1+1=(tr(ρZ2i+Λρ2iM)-0.5μM)；uni+1=tρr+(nM(ρZ-2in++Λ12i)), n≥2；6）更新Z:Zi+1=PS+(Gi+1-ρ1Λi)，Ti+1=★★★★★★★W~Xii++11(T~X(iu+i+1)1H)★★★★★★★；7）更新τ:τi+1=max{ηi+1-ρ1λi,0}；8）更新Λ和λ:Λi+1=Λi+ρZi+1-ρGi+1；λi+1=λi+ρ(τi+1-ηi+1)；i=i+1；输出:u_hat=ui+1；X_hat=~Xi+1/ηi+1。
2.3 计算复杂度分析
MUSIC算法的计算复杂度由协方差矩阵的计算O(M2L)、特征分解O(M3)和谱搜索O(IM2)组成，复杂度为O(M2L+M3+IM2)，主导项为O(IM2)，I是角度网格的点数；ESPRIT算法的计算复杂度无需谱搜索，其他和MUSIC算法一样，加上SVD分解O(K3)，复杂度为O(M2L+M3+K3)；ML算法的计算复杂度要迭代计算，单次梯度计算O(M2)，复杂度为O(IM2)；ADMM-MM算法的计算复杂度每次迭代计算半正定矩阵的投影和矩阵求解O(M3+M2L)，复杂度为O(T(M3+M2L))，T为ADMM-MM算法的迭代次数；而本文算法引入MM算法降低了ADMM的迭代次数，复杂度为O(T1(M3+M2L))，其中T>T1。
3 仿真分析与验证
仿真实验1:DOA估计精度分析
均方误差作为评估不同算法DOA估计精度性能的一个指标，定义为
式中，K表示角度估计的目标个数，θk表示第k个目标的真实角度值，
表示第k个目标的角度估计值。下面实验围绕MSE来评估角度估计的精度性能。
本文算法与ANM-ADMM算法、MUSIC算法、ESPRIT算法、ML算法的DOA估计精度进行对比分析。初始条件阵元数M=16，关闭阵元1，2，6，7，9，13和14，快拍数100，惩罚参数ρ=0.3，正则化参数μ=5。分别在信噪比为0 dB和15 dB对角度[-30°，5°，45°]的目标进行单次的DOA估计，结果分别如表2、表3和图2、图3所示。由表2和表3可以看出，本文算法与其他算法在相同条件对多目标的DOA估计中MSE最低，估计的精度最高。相较于ANM-ADMM算法，本文算法在信噪比为0 dB下，DOA估计精度提高约55%；在信噪比为15 dB下，DOA估计精度提高约35%。由图2和图3可以看出，本文算法和ANM-ADMM算法与其他3种传统算法相比，识别率更高，更加突出功率谱主峰的位置。
表2 信噪比为0 dB下的DOA参数估计 （°）
真实值-30545---ESPRIT-36.5562.22355.02650.406方法信号1信号2信号3MSE本文算法-29.0864.61444.0160.650 9 ANM-ADMM-28.5234.75443.5501.448 2 MUSIC-27.524.5141.496.236 9 ML-28.914.7243.840.870 7
表3 信噪比为15 dB下的DOA参数估计 （°）
真实值-30545---MUSIC-29.324.4144.080.552 3方法信号1信号2信号3MSE本文算法-29.2404.91244.6680.286 8 ANM-ADMM-29.1214.59544.3780.441 2 ESPRIT-33.7393.27448.0824.689 3 ML-29.214.8544.450.316 4
图2 信噪比为0 dB下的DOA估计
图3 信噪比为15 dB下的DOA估计
仿真实验2:抗噪声性能分析
均方根误差作为评估不同算法DOA估计抗噪声性能的一个指标，定义为
式中，J为算法的蒙特卡洛试验次数。下面实验围绕RMSE来评估角度估计的抗噪声性能。
为了验证本文算法的抗噪声性能，在仿真实验1的条件下将本文算法和MUSIC算法、ANM-ADMM算法、ESPRIT算法、ML算法的RMSE进行500次蒙特卡洛实验的对比验证，结果如图4所示。由图4可以看出，在-5~15 dB，本文算法较其他算法均方根误差最小，抗噪声性能最好。相对比在信噪比为15 dB条件下，在信噪比为0 dB条件下本文算法对多目标的DOA估计抗噪声性能较差。本文算法在信噪比为0 dB下，均方根误差约为0.93°左右；在信噪比为15 dB下，均方根误差约为0.32°左右，均满足对多目标DOA估计的强抗噪声性能在均方根误差1°以内，即该算法在信噪比为0~15 dB下均可实现对多目标有效的DOA估计，并且抗噪声性能好。
图4 500次蒙特卡洛实验的RMSE分析图
仿真实验3:运行时间对比分析
为了验证本文算法的计算量比ANM-ADMM算法低，在仿真实验1的条件下，统计次数设置为20次，将本文算法和MUSIC算法、ANM-ADMM算法、ESPRIT算法、ML算法的运行时间进行实验的对比验证，结果如表4和图5所示，表4为不同算法平均运行时间，图5为不同算法单次运行时间结果。通过表4和图5可知，虽然ESPRIT算法平均运行时间最短，但其对多目标的DOA估计精度最差，进而可能会导致多目标检测错误；MUSIC算法与ML算法运行时间虽然较本文算法短，但差距不是很大，而本文算法对多目标的DOA估计精度和抗噪声能力要强于它们，因此本文算法兼顾高精度、抗噪声强与计算量较低的优点；本文算法在ANM-ADMM算法的基础上进行改进，降低了计算量，与ANM-AD-MM算法相比较，运行时间下降约74%。
表4 不同算法平均运行时间
算法平均运行时间/s本文算法0.169 6 ANM-ADMM0.610 2 MUSIC0.127 3 ESPRIT0.000 4 ML0.075 9
图5 不同算法单次运行时间结果图
综上所述，在DOA估计精度性能上，本文算法的估计精度相比其他4种算法更高。对比ANM-ADMM算法，本文算法的DOA估计精度在信噪比为0 dB下，提高约55%；在信噪比为15 dB下，提高约35%。在DOA抗噪声性能上，本文算法在整体上的RMSE比其他4种算法低，具有更好的抗噪声性能。在DOA估计运行时间分析上，本文算法的运行时间虽然比MUSIC算法、ESPRIT算法、ML算法要长，并且ESPRIT算法运行时间最短，但估计的精度和抗噪声性能最差，不能有效地进行多目标DOA估计，MUSIC算法和ML算法的估计精度和抗噪声性能均没有本文算法好；相对比ANM-AD-MM算法，运行时间下降约74%。
4 结束语
本文算法在稀疏阵列下实现DOA估计，结合了ADMM算法和MM优化算法的优点。仿真实验表明，本文算法具有DOA估计精度高、抗噪声强、计算量较低等优点，相比于MUSIC算法、ESPRIT算法和ML算法，提高了DOA估计的精度和增强了抗噪声的能力；相比于ANM-ADMM算法，DOA估计精度提高约35%，运行时间下降约74%。
［1］ 杨洁艺.基于非均匀阵列的DOA估计算法研究［D］.西安:西安电子科技大学，2024.
［2］ ZUO Weiliang, XIN Jingmin, LIU Changnong, et al. Improved Capon Estimator for High-Resolution DOA Estimation and Its Statistical Analysis［J］. IEEE Journal of Automatica Sinica, 2023, 10（8）:1716-1729.
［3］ 陈柱学，吴建新.非均匀线阵的快速求根MUSIC方法［J］.雷达科学与技术，2014，12（2）:166-170.
［4］ 郑志东，张剑云.基于ESPRIT的MIMO雷达测向方法［J］.雷达科学与技术，2009，7（3）:205-209.
［5］ 陈婷，燕慧超，王鹏，等.基于SCAGOA优化BP神经网络和极大似然算法的DOA估计研究［J］.中北大学学报（自然科学版），2021，42（2）:177-186.
［6］ 杨青丽，成志强，熊瑛.基于空间频率估计的稀疏阵列测向［J］.电子信息对抗技术，2024，39（4）:19-24.
［7］ 李宝山，徐海文，陈晨,等.一种包含组合范数惩罚项的波达方向稀疏估计方法［J］.四川大学学报（自然科学版），2024，61（3）:93-100.
［8］ 张胜楠，刘峥，谢荣，等.基于ANM的非均匀圆阵二维DOA估计［J］.雷达科学与技术，2023，21（6）:630-636.
［9］ 龙伟军，徐艺卓，张玉禄，等.基于原子范数的非均匀阵列到达角估计［J］.现代电子技术，2024，47（17）:10-18.
［10］ 程智勇，陈胜垚，吴文，等.基于原子范数最小化的单比特稀疏双极子阵列的波达角估计［J］.电波科学学报，2023，38（2）:227-236.
［11］ 张燕，程先宬，陈伯孝，等.大孔径超稀疏阵列综合算法研究与应用［J］.雷达科学与技术，2020，18（4）:452-456.
［12］ 徐明，唐倩.基于原子范数的波达方向与时延联合估计方法［J］.通信学报，2024，45（10）:153-162.
［13］ 陈根华，罗晓萱.基于原子范数的互质阵列协方差矩阵重构算法［J］.太赫兹科学与电子信息学报，2023，21（3）:332-339.
［14］ 江首德，鄢社锋，毛琳琳，等.基于迭代原子范数最小化的均匀圆阵方位估计［J］.系统工程与电子技术，2025，47（6）:1746-1756.
［15］ 徐乃清，张劲东，李晨，等.基于MM算法的脉冲串模糊函数设计方法［J］.雷达科学与技术，2020，18（6）:599-604.
［16］ GIANELLI C,XU Luzhou,LI Jian,et al. One-Bit Compressive Sampling with Time-Varying Thresholds: Maximum Likelihood and the Cramér-Rao Bound ［C］//2016 50th Asilomar Conference on Signals, Systems and Computers,Pacific Grove, CA, USA:IEEE, 2016:399-403.
［17］ BHASKAR B N, TANG Gongguo, RECHT B.Atomic Norm Denoising with Applications to Line Spectral Estimation［J］.IEEE Transactions on Signal Processing,2013, 61（23）:5987-5999.
［18］ REN Jiaying, ZHANG Tianyi, LI Jian, et al. Sinusoidal Parameter Estimation from Signed Measurements via Majorization-Minimization Based Relax［J］. IEEE Transactions on Signal Processing, 2019,67（8）:2173-2186.
［19］ 王立府，王鹏.基于深度学习的原子范数最小化无网格DOA估计方法［J］.测试技术学报，2025，39（2）:218-229.