0 引 言
近年来，随着雷达、无线通信、声呐和地震勘探等领域对信号处理技术的不断追求，波达方向(Direction of Arrival,DOA）估计问题逐渐成为学术界和工业界研究的热点。传统的DOA估计算法例如旋转不变性技术(Estimating Signal Parameter via Rotational Invariance Techniques,ESPRIT），或者是多重信号分类(Multiple Signal Classification,M USIC），因其高分辨率和较低的计算复杂度被广泛应用于各种信号探测设备。
然而在实际雷达应用中，面对非高斯噪声和杂波干扰，传统算法因依赖协方差矩阵的精确估计而性能严重衰减。而高阶累积量（Higher-Order Cumulants,HOC）不仅可以极大扩展阵列孔径，还对高斯噪声具有天然的抑制能力，这使得算法能有效提取非高斯信号的统计特征，在低信噪比(Signal-to-Noise Ratios,SNR）下仍能保持鲁棒性［1］。
文献［2］采用互质L型阵列，构造四阶累积量并利用张量分解实现了闭式的二维DOA估计。该算法不仅解决了样本噪声扰动问题，还提高了阵列自由度，但其张量分解运算复杂度极高。文献［3］ 提出了一种基于四阶累积量的双网络多源DOA估计算法，可以在不需要精确阵列结构信息和无需阵列校准的情况下完成多源DOA估计，但因其天线分置并不适合小规模部署。文献［4］先对阵列接收信号特征分解，提取主信息分量，构建高斯窗修正的虚拟波束约束，再结合四阶累积量进行DOA估计，但其存在检测成功率不够稳定的问题。由以上研究可知，基于高阶累积量的DOA算法，虽可成倍扩展阵列孔径且抑制噪声干扰，但大多存在运算复杂度高的缺点。
此外，在实际复杂地形中因信号多径传播等干扰因素，会出现信号协方差秩亏损问题。而交替方向乘子法(Alternating Direction Multiplier Method,ADMM）作为一种分布式优化框架，特别适合求解高维、非光滑或带约束的DOA问题。该算法通过变量分裂将复杂问题拆解为多个子问题，并交替迭代求解，可在较少迭代次数内实现良好的收敛性［5］。进一步灵活地嵌入Toeplitz结构约束、非凸惩罚函数约束，将显著提升模型的适应性，使得在剧烈噪声干扰、信号相干、阵元失效、数据缺失的情况下仍能利用现有数据实现DOA估计，确保了DOA估计的稳定性。
文献［6］提出一种基于虚拟阵列扩展的加权空间平滑算法，通过虚拟阵列扩展补偿子阵划分导致的孔径损失，构建基于子阵互协方差信息的加权矩阵，在低信噪比下有较高精度。文献［7］提出一种基于二维重构的改进Toeplitz矩阵集重构法，通过将Toeplitz矩阵集与共轭转置矩阵相乘并求和，避免额外的去噪处理。该算法虽有较低复杂度，但面对多源相干信号性能较差。文献［8］将接收信号协方差矩阵分解为低秩无噪协方差矩阵和稀疏噪声协方差矩阵，构造凸优化问题求解有限采样下协方差矩阵估计误差问题。但其误差参数对性能影响显著，且其算法中的凸优化问题涉及高维矩阵分解，因此实时性受限。
本文结合降维高阶累积量与Toeplitz矩阵重构，通过理论分析和仿真实验，提出一种新型多源相干信号DOA估计算法。首先构建四阶累积量扩展阵列孔径，削弱高斯噪声干扰；之后提出一种科学高效的降维方法，既有更好的噪声鲁棒性，还有效降低了运算复杂度；然后提出一种自适应信噪比ADMM低秩Toeplitz矩阵重构法，将重构信号矩阵经MUSIC谱峰搜索得到DOA估计结果。不同场景下的仿真实验表明，本文算法在较低的计算复杂度下，在非高斯噪声环境中对多源相干信号亦有精确的估计结果。
1 阵列信号模型
如图1所示，一个配置了线性阵列(Uniform Linear Array,ULA）的雷达系统，由1个发射和M个接收阵元组成，阵列阵元间距为d=λ 2，λ为信号波长。假设t时刻存在P个远场窄带信号，分别位于方向θ1,θ2,…,θP。
图1 阵列信号模型
回波信号经匹配滤波后，接收数据可表示为
式中，A=[a1,a2,…,ap]∈ℂM×P是联合导向矩阵，其中ap是第p个目标的导向向量：
而s(t)=[s1(t),s2(t),…,sP(t)]T∈ℂP×1是匹配滤波后的信号矩阵。n(t)∈ℂM×1是加性噪声向量，其中包含高斯噪声、杂波干扰、脉冲噪声和相位噪声等多种分量。当接收信号包含L个快拍时，此时的观测数据矩阵X∈ℂM×L可表示为
式中，S=[s(t1),s(t2),…,s(tL)]∈ℂP×L，N=[n(t1),n(t2),…,n(tL)]∈ℂ M ×L。
2 算法原理
2.1 MUISC算法
以图1阵列信号模型为例，对于式（3）的阵列接收信号，可得信号的理想协方差为
式中，RS=E[S(t)SH(t)]为信号协方差矩阵，σ2I为噪声协方差矩阵，σ2为噪声方差，I为单位矩阵。显然协方差矩阵由信号S(t)与噪声N(t)两部分组成。此时信噪比估计公式为
式中，‖X‖F为X的F范数，即
实际应用中，由于采样数据有限长，我们通常使用观测数据样本的协方差进行DOA估计：
定义MUSIC算法功率谱为
其前P个峰值即DOA角度θ。同时根据上述推导，MUSIC算法至多可估计M -1个信源。
2.2 降维高阶累积量算法
对于零均值的平稳复随机过程x，其一种四阶累积量(Forth Order Cumulant,FOC）定义为
式中，cum(·)表示变量的累积量，(·)*表示矩阵复共轭，E[x(t1)x(t2)x*(t3)x*(t4)]表示x(t)的四阶矩，E[x(ti)x(tj)]表示x(t)的二阶矩。
对于高斯随机过程y，其高阶矩可描述为
显然，高斯随机变量的高阶累积量(k≥2)恒等于0，且其高阶矩只取决于二阶矩σ2［9］。
由式（4）可知，协方差矩阵RXX可分解为噪声和信号子空间协方差矩阵的组合。同理，四阶累积量可实现类似分解，根据式（9）的四阶累积量形式，此时令第p个目标的累积量矩阵为
再将其扩展至L个快拍的情况可得：
式中，Xp=apsp+np∈ℂM×L为第p个目标的观测数据矩阵，sp=[sp(1),sp(2),…,sp(L)]∈ℂ1×L和np∈ℂ1×L分别为第p个目标的多快拍信号和噪声，CS和CN分别为目标信号的四阶累积量；cP为按式（13）分解后，CS的对角线元素，⊙表示Khatri-Rao积。其中：
显然，对比原阵列，四阶累积实现了孔径扩展。与2.1节MUSIC算法类似，定义空间谱：
对式（17）进行谱峰搜索，即可得到DOA角度。但对于原信号模型中M个接收阵元，可知累积量矩阵C共有M4个元素。因此式（17）的计算过程复杂度极高。为降低复杂度，文献［10］给出了式（18）所示的逐元素快速去冗余办法：
式(18）的解算过程，如式（19）～（20）所示：
式中，R4为降维矩阵。但上述降维过程运算繁琐，且由式（10）可知高斯噪声的四阶累积量为0，因此文献［10］中降维矩阵是在高斯噪声下推导得出，未考虑实际工作情况中矩阵降维前后噪声与信号功率的一致性。因此本文对四阶累积量降维算法提出一种更为高效科学策略。
对于四阶累积量矩阵，定义：
式中，Hm=[0M×(M -m),IM,0M×(m -1)]，IM为M维单位矩阵，0M×(M -m)为M×(M -m)的零矩阵，m=1,2,…,M -1。此时导向矢量可写为
式中，F(θ)=[f(θ1),f(θ2),…,f(θP)]∈ℂ(2M-1)×P为去冗余阵元后虚拟阵列的阵列流型矩阵，其中f(θ)=[e-j(M-1)sin θ,e-j(M -2)sin θ,…,ej(M -1)sin θ]。
考虑到降维变换前后信号及噪声功率的一致性，定义降维变换矩阵：
表1为本文降维法与文献［10］逐元素降维法复杂度对比。逐元素降维法对累积量矩阵的行列两个维度分别进行降维，相比之下本文算法只计算一次降维变换矩阵，显然有更低的运算复杂度。
表1 降维算法运算复杂度对比
2M -1 2M 2M -1 2M算法 复杂度本文恒功率降维法 O[2()]文献［10］逐元素降维法 O[4()]
将式（23）代入式（13）得降维的四阶累积量：
此时累积量矩阵CRD从M2×M2降维至(2M -1)×(2M -1)。对CRD特征分解，取其2M -1-P个小特征值对应的特征向量构造噪声子空间UN，结合降维的导向矢量矩阵F(θ)，定义降维累积量空间谱：
对式（25）进行谱峰搜索，即可得到DOA角度。表2为降维前后运算复杂度的对比，可见降维后运算复杂度得到指数级的降低。
表2 降维前后运算复杂度对比
M2 3+O■■()■■P2M4 2+O[]2M2log M()(P2())2M -1 2 2+降维后 2M -1 3+O■■■■■■■■O[(2M -1)log(2M -1)]降维前后 复杂度降维前 O■■()■■O[]
2.3 ADMM低秩Toeplitz矩阵重构算法
ADMM是一种用于求解带有约束的凸优化问题的迭代方法，适用于式（26）的问题形式：
式中，x,z为优化变量，f,g为待优化凸函数。
它将问题拆分为若干个子问题，分别求解，使得最终解满足所有约束。其优点是可以将复杂问题分解为较简单的子问题，每个子问题都能有较为简单的更新公式［11］。
在阵列信号处理中，我们通常会利用式（7）的多快拍数据计算样本协方差矩阵。对于均匀线阵，理想协方差矩阵应满足半正定性（所有特征值非负）和Toeplitz结构（沿对角线元素相同）。但因噪声干扰和信号相干，实际协方差矩阵不满足上述特性，因此需要进行矩阵重构，辅以Toeplitz约束，以解决信号矩阵的秩亏问题［12］。
假设R为重构矩阵。其同时满足Toeplitz结构且与
能量相近，考虑如下优化问题：
式中，
为拟合项，
以F范数衡量重构矩阵R与样本协方差
的差异；λ‖R‖*为低秩正则项，‖R‖*为以核范数代替的凸松弛约束函数；λ为正则化参数，用于平衡数据拟合项与低秩约束；
表示所有Toeplitz矩阵的集合，
表示所有半正定矩阵的集合。
在动态环境中，SNR变化会影响噪声功率和信号相干性，因此本文算法引入λ自适应调整机制以优化重构性能。根据式（5）得到SNR估计值后，对其进行指数衰减映射：
式中，λmin和λmax为算法给定的上下界，α为衰减速率。即低信噪比时强化低秩约束，抑制噪声引起的矩阵秩升高；高信噪比保真数据拟合，避免过度约束导致信号失真。
采用ADMM求解上述优化问题，首先引入一个辅助变量Z，式（27）问题等价为
即将数据拟合项作用于R，低秩约束和结构约束作用于Z，再利用Z=R将二者相关联。
根据式（29）构造增广拉格朗日函数：
式中，
Y,R -Z
为内积，目的为确保约束Z=R被满足，其中Y是拉格朗日乘子；
是额外的惩罚项，ρ＞0是惩罚参数。
ADMM迭代求解矩阵重构的更新步骤如下：
1）固定当前Zk和Yk，更新R，即求解
式(31）为一个最小二乘问题，其梯度为
其闭式解为
解为Toeplitz集合上的加权平均。
2）固定Rk+1和Yk，更新Z
式中，U=Rk+1+Y
，对U做奇异值阈值(Singular Value Threshold,SVT）操作，即先对U执行奇异值分解：
结合引入的动态调整λ的策略，再对得到的奇异值同步进行软阈值处理，将每个奇异值σi更新为
由此重构得
之后先对Z'取均值进行Toeplitz投影，再对该矩阵特征分解，投影到半正定集合。重复此步骤，直至收敛得到Zk+1。
缩小R与Z差距，逼近Z=R约束。
图2为矩阵重构前后的对比，由图可知重构矩阵与原矩阵的误差很小，既实现了矩阵的半正定性和Toeplitz性，还尽量保留了原矩阵信息。
图2 矩阵重构前后对比
表3展示了ADMM法Toeplitz重构和经典Toeplitz重构复杂度的对比，可见本算法主要复杂度由迭代次数k和信源数P决定，换言之ADMM带来的计算开销很小。
表3 不同矩阵重构法复杂度对比
算法 复杂度ADMM法Toeplitz矩阵重构 O[]kPM2经典Toeplitz矩阵重构 O[]M 2
综上，本文基于降维累积量和ADMM低秩重构的多源相干DOA估计算法如表4所示。
表4 降维累积量ADMM低秩Toeplitz矩阵重构算法
输入：观测数据矩阵X；输出：DOA估计结果θ；①计算观测数据矩阵X的四阶累积量C4，并得到虚拟阵列的导向矢量矩阵B；②构造降维变换矩阵H和T，对导向矢量B和四阶累积量C4进行降维变换，得到降维矩阵CRD和降维导向矢量F；③初始化R0=CRD，Z0=CRD，Y0=0；for k=0,1,2,…RUYTiΣ͂fo kUk‖e=++R=p11Rlk==id+tkiz1+[Ya投R-1ĝk+[Z+影+mYkρρa和+kx((1ρR(Z半‖σ2F；kk+正i＜-1-ε定-Yλ bkZ投rρρeka影,+)0k]1) )迭]；(，1Z代+'得=ρU)Z；Uk Σ+͂1U；VUH；end④对迭代矩阵R特征分解，取前P个大特征值对应特征向量构建噪声子空间UN；⑤使用MUSIC算法，利用降维导向矢量F和噪声子空间UN谱峰搜索得到信源θ。
3 仿真与分析
本文对传统DOA算法在实际应用中的性能衰减做出了诸多改进。本节对算法进行仿真验证，设置不同的仿真场景，与同类算法进行横向对比，验证了本算法在复杂场景中优良的性能表现。
实验1 验证本文算法经过矩阵重构之后，对相干信号的估计表现。仿真参数如下：非高斯噪声条件、SNR = -10、快拍数为100、发射和接收阵元分别为1和8、本文算法迭代1次、非高斯相干信源。仿真结果如图3和图4所示。
图3 矩阵重构前后估计性能对比
图4 矩阵重构前后RMSE对比
由图3可知，未经矩阵重构的算法面对相干信号会丢失目标，而经过重构的算法在低信噪比和少快拍下都保持了良好的性能。最后根据图4的均方根误差(Root Mean Square Error,RMSE）变化曲线可看出非高斯噪声和相干信源条件下经过矩阵重构算法的RMSE水平更低，消融实验结果说明了矩阵重构的必要性。
实验2 验证算法在高斯噪声下，分别在不同信噪比时对相干信号的估计表现。仿真参数如下：高斯噪声条件、快拍数为100、发射和接收阵元分别为1和8、本文算法迭代1次、非高斯相干信源。仿真结果如图5和图6所示。
图5 高斯噪声下相干信号估计性能对比
图6 高斯噪声下相干信源不同信噪比RMSE
由图5（a）可知，在SNR = 0时，文献［14］算法与本文性能相近，而文献［13］算法会丢失目标。图5（b）当SNR = -10时，前两者算法均有较大误差。且文献［13］算法因空间平滑存在孔径损失，因此在图5（c）中，信源数增加至8时，该算法无法运行。最后根据图6的RMSE变化曲线可看出高斯噪声条件下本文算法的RMSE水平更低，实验结果说明了本文算法解相干的有效性。
实验3 验证算法在高斯噪声下，分别在不同信噪比时对非相干信号的估计表现。仿真参数如下：高斯噪声条件、本文算法迭代1次、快拍数为100、发射阵元为1、接收阵元为8、8个非高斯非相干信源。仿真结果如图7和图8所示。
图7 高斯噪声下非相干信号估计性能
图8 高斯噪声非相干信源不同信噪比RMSE
由图7（a）可知，在SNR = 0下文献［15］算法与本文算法性能相近。但在图7（b）中，随着信噪比持续降低，因波束指向性变差，表现逐渐弱于本算法。极端情况下，因其较高的旁瓣水平，可能会出现伪峰，造成严重的估计错误。对比图8中曲线可看出，高斯噪声条件下，在信噪比变化区间，尤其是低信噪比段，相比其他算法，本算法的RMSE水平更低，实验结果说明了本文所提算法的稳定性。
实验4 验证算法在非高斯噪声的环境中在不同信噪比下的估计性能表现。仿真参数如下：非高斯噪声条件（50%高斯噪声，40%杂波干扰，10%脉冲噪声，10%脉冲噪声率、10%相位噪声率）、本文算法迭代次数为1、快拍数为100、发射阵元为1、接收阵元为8、8个非高斯非相干信源。仿真结果如图9至图12所示。
图9 非高斯噪声下非相干信号估计性能
对比图9（a）和图9（b），非高斯噪声下文献［12］算法劣化严重，文献［9］算法受益于传播算子的特性，某些情况下相比之前有更好的性能，但低信噪比下会丢失目标。
由图10可知，非高斯噪声下文献［14］算法会出现极大误差，即便SNR = 0时也会丢失目标。
图10 非高斯噪声，SNR=0，相干信号估计性能
对比图11和图12可知，无论是否为高斯噪声或相干信号，本文算法的RMSE都更低，实验结果说明了本文所提算法具有较好的鲁棒性。
图11 非高斯噪声下相干信源不同信噪比RMSE
图12 非高斯噪声非相干信源不同信噪比RMSE
实验5 验证各算法在高斯噪声的环境中，分别对相干信号在不同快拍数下的估计性能。仿真参数如下：高斯噪声条件、本文算法迭代次数为1、信噪比分别为-10和0、快拍数为100～1 000、发射阵元为1、接收阵元为8、信源为非高斯相干信号。仿真结果如图13所示。
图13 高斯噪声相干信源不同快拍数RMSE
根据图13（a）可知，各算法的误差都会随快拍数增加而减小，其中文献［14］算法在高信噪比和多快拍条件下与本算法接近。但在图13（b）中低信噪比的条件下，本文算法在少快拍数下明显更具优势，随快拍数有更快的收敛趋势。究其原因，本算法得益于引入的ADMM法对多快拍数据的迭代处理，其估计结果的RMSE水平更低，实验结果对比说明了本文算法的高效性。
实验6 验证各算法在非高斯噪声的环境中分别对相干和非相干信号在不同快拍数下的估计性能。仿真参数如下：非高斯噪声条件（同实验4）、本文算法迭代次数为1、信噪比为0、快拍数为100～1 000、发射阵元为1、接收阵元为8。仿真结果如图14所示。
图14 非高斯噪声，不同快拍数RMSE
由图14可知，在非高斯噪声场景下，本文算法的凭借迭代矩阵重构，在1～200快拍区间有更好的性能表现。实验结果表明本算法在更少快拍下有更低的初始RMSE水平，实验结果验证了本文算法的优越性。
实验7 验证本算法在高斯与非噪声的环境中，分别在不同迭代次数下的估计性能表现。仿真参数如下：信噪比分别为-10和0、本文算法迭代次数为1～100、快拍数为100、发射阵元为1、接收阵元为8、非高斯相干信源。仿真结果如图15所示。
图15 高斯与非噪声相干信源不同迭代次数RMSE
由图15（a）可知，在低信噪比下，1～20次的迭代即可显著降低估计结果的RMSE；而在高信噪比的条件下，1～2次迭代即可得到理想效果。由图15（b）可知，相比高斯噪声的环境，非高斯噪声条件下虽然准确度有所降低，但是体现出与高斯噪声环境中相似的趋势。综上实验结果可知本算法大部分情况下只需1～10次的迭代即可得到优良性能，很好平衡了算法复杂度与精确度。
4 结束语
本文研究了非高斯噪声干扰下多源相干信号的DOA估计问题，提出了一种结合降维高阶累积量和自适应信噪比ADMM法低秩Toeplitz矩阵重构的算法。首先计算阵列接收数据的四阶累积量，以实现削减高斯噪声和扩展阵列孔径；之后提出一种高效科学的降维策略，不仅提高了在非高斯噪声环境下的稳定性，还有效降低了运算复杂度；之后通过ADMM法高效迭代，灵活嵌入自适应信噪比约束和低秩Toeplitz约束，得到重构矩阵；最后利用MUSIC进行谱峰搜索得到目标精确的估计值。算法不仅实现了孔径扩展，还可在不损失阵列孔径的前提下实现多源相干信号探测，即便在低信噪比和非高斯噪声干扰的条件下依然保持良好的估计性能，实验验证了本算法在复杂干扰下DOA估计有着更强鲁棒性。本研究提出的算法不仅为雷达的阵列信号处理提供了新的解决思路，还为非理想环境中DOA估计平衡性能与效率提供了创新性技术路径。
［1］ MASSOUD A,IQBAL U,NOURELDIN A.Warped Fourth Order Cumulant Beamforming for Estimating DOA from Non-Uniform Linear Arrays of DIFAR Sonobuoys［C］//2024 6th International Conference on Communications,Signal Processing,and Their Applications,Istanbul,Turkiye:IEEE,2024:10794275.
［2］ ZHENG Hang,ZHOU Chengwei,WANG Yong,et al.Sample Fourth-Order Cumulant Tensor Denoising for DOA Estimation with Coprime L-Shaped Array［C］//2021 55th Asilomar Conference on Signals,Systems,and Computers,Pacific Grove,CA,USA:IEEE,2021:1668-1672.
［3］ ZHANG Chenhao,WANG Wenjie,HONG Xi,et al.A Fourth-Order Cumulant-Based Multi-Sources DOA Estimation in UAV Collaborative Systems［J］.IEEE Signal Processing Letters,2024,31:251-255.
［4］ 徐峰,杨小鹏,赵毅.基于快速张量分解的波束空间MI-MO雷达二维DOA估计算法［J］.信号处理,2022,38（1）:1-8.
［5］ 李荣禄，汤建龙，袁永强.基于TOEPLITZ重构的压缩感知嵌套阵列DOA估计［J］.雷达科学与技术,2024,22（3）:334-340.
［6］ 孙志国,柏乔森,白永珍,等.基于虚拟阵列扩展的改进加权空间平滑算法［J］.系统工程与电子技术,2023,45（1）:250-256.
［7］ JIN Aisuo,ZHANG Fei.2D DOA Estimation of Coherent Sources Based on Reconstruction of Toeplitz Matrix Sets［C］//2020 Chinese Automation Congress,Shanghai，China:IEEE,2020:508-513.
［8］ 王洪雁,于若男.基于稀疏和低秩恢复的稳健DOA估计方法［J］.电子与信息学报,2020,42（3）:589-596.
［9］ 窦慧晶,郭威,谢东旭,等.基于四阶累积量的改进传播算子DOA估计算法［J］.北京工业大学学报,2024,50（9）:1061-1068.
［10］ 刁鸣,吴小强,李晓刚.基于四阶累积量的测向方法研究［J］.系统工程与电子技术,2008,30（2）:226-228.
［11］ KLIOUI Y.Toeplitz-Hermitian ADMM-Net for DOA Estimation［PP/OL］.arXiv,2025.https://doi.org/10.48550/arXiv.2502.13896.
［12］ WEN Kelei,TIAN Ye,QIN Yunbai,et al.Toeplitz Rectification Based DOA Estimation Under the Coexistence of Mutual Coupling and Nonuniform Noise［J］.AEU-International Journal of Electronics and Communications,2023（1）:154453.
［13］ TIAN Gege,ZHAO Chenglu,LI Jingxiao,et al.The DOA Estimation Method Based on Improved SSMUSIC Algorithm［C］//2022 3rd China International SAR Symposium,Shanghai,China:IEEE,2022:9971459.
［14］ 张豪,甄冬,刘英辉,等.采用四阶累积量多重矩阵重构的低信噪比相干声源波达方向估计［J］.声学学报,2023,48（2）:337-346.
［15］ 刘甲磊,马佳智,施龙飞.虚拟波束四阶累积量DOA估计方法［J］.系统工程与电子技术,2022,44（7） :2134-2142.