阵元失效下稀疏阵列的二维DOA估计算法

摘 要：本文针对二维稀疏阵列在阵元失效条件下， 因数据缺失导致虚拟阵列连续性被破坏及自由度下降的问题， 提出了一种二维DOA估计算法。 首先基于二维差分共阵构建虚拟阵列， 然后利用解耦原子范数最小化理论， 以矩阵填充的形式恢复协方差矩阵数据， 实现对虚拟阵列中丢失虚拟阵元的内插， 最后采用SS-MUSIC算法进行多信源的二维DOA估计。 所提方法弥补了物理阵元失效所造成的影响， 恢复了原始虚拟阵列的完整孔径特性， 保持了虚拟阵列的自由度， 从而确保了较高精度的二维DOA估计性能。 仿真实验结果表明， 在相同阵元数量及阵元失效情况下， 本文提出的算法相比已有方法能有效地估计更多信源， 并在小快拍数和低信噪比条件下表现出更高的稳健性， 最大限度地保留并利用了稀疏阵列在二维DOA估计中的自由度优势。

关键词：  二维DOA估计； 稀疏阵列； 差分共阵； 阵元失效； 解耦原子范数最小化； 矩阵填充

中图分类号： TJ765； TN911.7

文献标识码： A

文章编号：  1673-5048（2024）02-0114-09

DOI： 10.12132/ISSN.1673-5048.2024.0042

0 引  言

波达方向（Direction of Arrival， DOA）估计是阵列信号处理领域内的一个极其重要的分支， 长期以来受到学者们的广泛关注。 在过去的几十年间， DOA估计理论不断发展， 从雷达探测、 无源定位、 无线通信、 卫星导航， 乃至物联网环境中的智能感知， DOA估计技术在众多领域均展现出广泛而深远的应用价值［1-3］。 以多重信号分类（Multiple Signal Classification， MUSIC）算法［4］为代表的子空间分解类算法凭借高精度和超分辨的DOA估计能力为人所熟知， 并获得了广泛的应用与实践验证。 然而这类算法的信号接收模型中通常依赖于均匀阵列， 其理论上可估计的最大信源数目受限于阵元的个数。

因此， 在日益复杂的空间电磁环境下， 针对信号密集场景的精确测向问题， 基于稀疏阵列的欠定DOA估计成为新的研究热点［5］， 稀疏阵列突破了阵元间距必须为入射信号半波长整数倍的限制， 通过设计阵列构型以构建大孔径虚拟阵列， 从而显著提升了阵列自由度（Degree of Freedom， DOF）， 进而能够实现超过实际物理阵元数量的DOA估计［6］。 Pal等学者开创性地提出嵌套阵列的概念， 并基于此发明出一种欠定DOA估计算法［7］。 该算法首先通过对原始接收数据进行协方差矩阵的向量化处理， 形成与虚拟阵列相对应的单快拍数据， 随后执行空间平滑（Spatial Smoothing， SS）操作获得满秩的协方差矩阵， 使用MUSIC算法完成角度估计。 随后Pal等人进一步拓展， 又提出了适用于二维空间的平面嵌套阵列设计及相应DOA估计算法［8-9］。 尽管SS-MUSIC算法通常可以应用于大多数类型的稀疏阵列DOA估计中， 其有效性却高度依赖于虚拟阵列结构的连续性［10］， 一旦虚拟阵列呈现出不连续性， 即其间存在孔洞， 那么阵列的连续自由度扩展将受到显著限制。

通过使用协方差矩阵重构技术， 能够对虚拟孔洞位置的空缺数据信息进行近似化填充， 从而扩大虚拟阵列孔径。 一系列基于互质线阵内插的DOA估计算法被提出， 文献［11］通过核范数最小化构建凸优化模型， 填补虚拟阵列孔洞并提升DOA估计精度， 但核范数方法可能引入一定近似误差； 文献［12］则进一步将核范数最小化问题转化为迹最小化以求解， 实现离网格DOA估计； 文献［13］对虚拟阵列内插算法进行了完整分析， 优化模型限制了此过程中的误差积累。 近年来， 原子范数理论因其对信号空域稀疏特性的有效利用， 在DOA估计领域取得了显著进展， 原子范数最小化（Atomic Norm Minimization， ANM）作为一种无网格方法， 避免了一般压缩感知方法中的基失配问题。 文献［14］将ANM算法应用于互质线阵的虚拟阵列内插， 实现了高精度的无网格DOA估计。 文献［15-16］则针对二维DOA估计提出解耦原子范数最小化（Decoupled ANM， DANM）算法， 基于二维阵列接收数据模型， 将二维DOA估计问题转化为两个独立的一维DOA估计； 文献［17］中对DANM算法做出改进， 使之适用于多快拍情形下的二维DOA估计。 但这两种方法都局限于均匀矩形阵列， 且最大可估计信源数受限于较短一边的阵元数。 文献［18］基于互质平面阵列， 采用DANM算法进行二维虚拟阵列内插， 可实现欠定的二维DOA估计。

航空兵器 2024年第31卷第2期

司伟建， 等： 阵元失效下稀疏阵列的二维DOA估计算法

在实际应用环境中， 由于硬件老化、 环境干扰等多种难以预测的因素导致性能退化乃至完全失效。 一旦阵元随机发生故障， 虚拟阵列的连续性很可能会因此受到显著破坏， 形成位置和大小不定的孔洞， 进而大幅降低连续自由度， 并对DOA估计精度产生严重影响［19］。 面对阵元失效情况下的DOA估计难题， 类似地， 可通过内插方法填补阵元失效后虚拟阵列中的孔洞， 近似恢复丢失虚拟阵元的数据信息， 以维持虚拟阵列的连续性， 从而最大限度地利用阵列原始的完整孔径优势。 文献［20］采用Toeplitz矩阵重构技术恢复协方差矩阵； 而文献［21］则基于重加权l2， 1范数最小化算法恢复了完整的阵列数据矩阵， 实现了一维DOA的有效估计。 然而这两种方法在将算法扩展到二维场景时仍面临较大困难。

鉴于当前研究现状， 阵元失效问题的研究重心主要集中在一维阵列上， 而在二维稀疏阵列DOA估计领域， 关于阵元失效下的虚拟阵列孔洞内插问题， 现有的研究尚不充分。 因此， 本文提出一种阵元失效下的二维稀疏阵列DOA估计算法， 该算法利用二维差分共阵形成虚拟阵列， 通过建立基于DANM理论的矩阵填充问题， 恢复协方差矩阵数据， 对阵元失效导致的孔洞内插， 以恢复原虚拟阵列的完整自由度， 并结合SS-MUSIC算法进行二维DOA估计。 对比传统算法， 本文方法能够避免舍弃部分虚拟阵元而造成的自由度损失， 有效维持了稀疏阵列二维DOA估计的优势性能。

1 信号模型

考虑空间中K个窄带远场独立信号入射到XOY平面中一个由M个标量天线构成的稀疏平面阵列上， 如图 1所示。 在平面直角坐标系中， 阵列中任一阵元的位置都可表示为（xid， yid）的形式， 其中： xi， yi∈瘙綄， i=1， 2， …， M； d=λ/2， 即半波长。 令ni=（xi， yi）∈瘙綄2， 代表单个阵元位置， 构成该稀疏平面阵列的所有物理阵元的二维位置整数集合记为P。

假设第k个信号的入射角度为（θk， φk）， 其中： θk为仰角， 即入射方向与Z轴的夹角， θk∈［0， π/2］； φk为方位角， 即入射方向在XOY平面上的投影与X轴的夹角，  φk=［0， 2π）。 第k个信号的导向矢量可表示为

aP（θk， φk）=［ax1， y1（θk， φk）， …， axM， yM（θk， φk）］T（1）

式中： axi， yi（θk， φk） = ejπ（xicosφksinθk+yisinφksinθk）， i = 1， 2， …， M， j是虚数单位， k=1， 2， …， K。

假设噪声为加性高斯白噪声， 且噪声与信号相互统计独立， 则在t时刻时， 稀疏平面阵列的信号接收数据模型为

XP（t）=APS（t）+N（t）=∑Kk=1aP（θk，φk）sk（t）+N（t） （2）

式中：  S（t）=［s1（t）， s2（t）， …， sK（t）］T为信号向量；N（t）=［n1（t）， n2（t）， …， nM（t）］T为噪声向量； 阵列流形矩阵AP=［aP（θ1， φ1）， …， aP（θk， φk）］∈瘙綇M×K。

当考虑全部T个快拍数据时， 式（2）转化为

XP=APS+N （3）

式中： XP∈瘙綇M×T， S∈瘙綇K×T， N∈瘙綇M×T。

计算阵列接收数据的协方差矩阵RP=E［XPXHP］， 其中E［·］表示数学期望， 将式（3）代入可得：

RP=APRSAHP+σ2NI=∑Kk=1σ2kaP（θk，φk）aHP（θk，φk）+σ2NI（4）

式中： RS为信号的自协方差矩阵， 其对角线元素表示信号功率； σ2N为输入噪声N的方差， I为单位阵。

由于采样数量有限， 通常采用最大似然估计方法对式（4）进行近似计算， 可表示为

R^P=1T∑Tt=1XPXHP（5）

对R^P进行向量化操作， 可得

z=vec（R^P）=（A*P⊙AP）p+σ2Nι（6）

式中： p =［p1， …， pK］T，  ι =vec（ I ）； ⊙表示Khatri-Rao积。 可将 z 视作一个单快拍的虚拟阵列接收信号矢量，  A*P⊙ AP则是虚拟阵列的流形矩阵， 大小为M2×K， 其第k列为 a*P（θk， φk） aP（θk， φk）， 表示Kronecker积， 将式（1）代入后展开可得

a*P（θk， φk）aP（θk， φk）=a*P（～x， k， ～y， k）

aP（～x， k， ～y， k）=［ejπ［（x1－x1）～x， k+（y1－y1）～y， k］， …，

ejπ［（xM-x1）～x， k+（yM-y1）～y， k］， ejπ［（x1-x2）～x， k+（y1-y2）～y， k］， …，

ejπ［（xM-x2）～x， k+（yM-y2）～y， k］， …，  ejπ［（x1-xM）～x， k+（y1-yM）～y， k］， …，

ejπ［（xM－xM）～x， k+（yM－yM）～y， k］］T（7）

式中： 为便于表示， 令～x， k=cosφksinθk， ～y， k=sinφksinθk。 通过观察式（7）中波程差的形式， 可以将二维虚拟阵列的阵元位置坐标表示为

D={（xi－xj， yi－yj）|（xi， yi）， （xj， yj）∈P，

i， j=1， 2， …， M}（8）

即D为P中任意两个坐标之间差的集合（包括i=j的情况）， 但集合中元素有唯一性， 同一个虚拟阵元位置可能是由多对物理阵元的位置差而得， 所以D中元素个数一定小于M2， 但虚拟阵列的阵元数量较原物理阵列的大幅提升将使二维欠定DOA估计得以实现。

对应地， 式（7）向量中的元素也会有重复， 于是根据与D中各阵元对应关系， 对式（6）中向量z进行去冗余操作， 即将其中对应同一虚拟阵元的元素取平均值再按对应坐标大小重新排序［10］， 则对应D的单快拍虚拟接收数据为

zD=ADp+σ2NιD（9）

式中： AD为虚拟阵列D接收数据的流形矩阵， 其大小为|D|×K（|·|表示集合的势， 即集合中元素的个数）， ιD表示式（6）中向量化单位阵vec（I）对应D的重排序。

2 阵元失效下二维DOA估计

2.1 二维差分共阵

根据式（8）对二维虚拟阵列的表示， 可以给出以下概念。 对于由集合P={ni=（xi， yi）|i=1， 2， …， M}指定的二维阵列， 其差分共阵D定义为P阵元位置坐标两两之间差的集合［22-23］：

D={m|m=ni－nj， ni， nj∈P}（10）

式中： m=（mx， my）∈瘙綄2表示虚拟阵元的平面位置坐标。

如此由M个物理阵元形成的差分共阵的阵元数最大可达|D|=M（M－1）+1［7］。 由式（10）易知， 差分共阵的阵元位置是关于原点（0， 0）中心对称的， 所以， 当分别取得虚拟阵元坐标中mx和my的最大和最小值时， 有mxmax=－mxmin， mymax=－mymin。 于是可以定义一个包含D中所有阵元的最小均匀矩形阵列V：