基于L型和差嵌套阵的二维波达方向估计

朱大君 何培宇 喻伟闯

（四川大学电子信息学院，四川成都 610065）

1 引言

在阵列信号处理领域，波达方向估计（Direction of Arrival，DOA）估计是至关重要的部分，在例如雷达，射电天文学和无线通信［1-4］等领域有着广泛的研究与应用。传统的DOA 估计研究大多集中在均匀线性阵列（Uniform Linear Array，ULA）上，由于ULA的自由度受阵元数限制，导致具有N个传感器的ULA，最多可检测N-1个信源。

为了提高检测能力，有学者提出稀疏阵列。2010 年Vaidyanathan 和Pal 提出了两种稀疏阵列，嵌套阵［5］（Nested Array，NA）和互质阵［6］（Coprime Array，CA），通过向量化接收到的信号协方差矩阵，使用N1+N2个阵元的NA 就可以达到2N2(N1+1) -1 的自由度（Degrees of Freedom，DOF）；2N1+N2-1 个阵元的CA 可以达到2N1(N2+1) -1 的自由度，远大于ULA 的N1+N2的自由度。但是阵元位置间的差分阵（difference co-array，DCA）存在孔洞，导致阵列的DCA 中央连续部分，才能应用空间平滑MUSIC［7-9］（Spatial Smoothing Multiple Signal Classification，SS-MUSIC）算法来估计入射信号的DOA。近期，有学者提出了一些改进版本的稀疏阵列，如超级嵌套阵［10］（Super Nested Array，SNA）、增强嵌套阵［11］（Augmented Nested Array，ANA）和广义嵌套阵［12］（Generalized Nested Array，GNA）等，来提高阵列自由度。

但这些稀疏阵列仅考虑差分阵列DCA，未考虑求和阵列［13-14］（sum co-array，SCA）。适当地利用阵列的SCA，可以填补DCA 中的一些孔洞，得到比仅考虑DCA 的稀疏阵列获得更高的DOF。2017 年Xinghua Wang 等人提出VCAM［15］（Vectorized Conjugate Augmented MUSIC，VCAM）算法，对物理阵元同时进行DCA 和SCA 操作，称为DSCA（difference and sum co-array，DSCA）。但是这些稀疏阵列，其DSCA的虚拟阵元间仍然存在很多冗余阵元和孔洞，导致阵列自由度扩展有限。为了更好的解决这些问题，本文提出了一种基于求和求差的改进嵌套阵，记为INA（Improved Nested Array，INA），使用N个阵元，就可以实现O(N2)的自由度。

另外，在DOA 估计应用时，不仅需要知道入射信号的方位角，还需要知道俯仰角，即需要进行二维的DOA 估计。在［16］中已经证明，L 型阵列不仅有简单阵列布置结构，由两个互相垂直的子阵组成，而且具有更好的估计性能。但是传统的L 型阵列的每个子阵都是ULA，具有低分辨率和物理阵元限制的问题。此外，二维（Two-Dimensional，2D）DOA 估计时需要使用2D-MUSIC算法的进行二维的角度搜索，其计算复杂度高。

因此本文将提出的INA 扩展到L 阵，进行二维的DOA 估计。此L 型阵列的两个子阵都是INA，通过VCAM 算法对阵元位置之间的差分、求和操作达到虚拟扩展阵元数目的效果，虚拟阵元数目具有确定的表达式。此L 阵能够提升阵列天线的自由度，更加准确的估计出更多的信源波达方向，突破了物理阵元的限制。为了避免二维角度搜索［17-18］，采用PSCM［19-20］（Pair-matching Signal Covariance Matrices，PSCM）算法，对L 阵两个子阵的一维DOA 估计角度，进行角度的匹配。文中第3.2 节给出自由度的数学表达式和第4 节实验仿真来说明此阵列的优异性。

本文符号说明，(.)T、(.)*和(.)H分别表示矩阵的转置、共轭和共轭转置。vec(.)代表向量化操作，⊙和⊗分别表示Khatri-Rao积和Kronecker积。

2 阵列信号模型

2.1 INA模型

常规的均匀线阵（ULA）进行DOA 估计时需要满足空间采样定理，即阵元间间隔不能大于接收信号的半波长，这限制了阵列的分辨率，对于N个阵元的ULA最多估计N-1个信号。为了解决上述问题，获得更高的自由度和分辨率，本文提出了一种基于求和求差的改进嵌套阵，记为INA。阵列结构如图1 所示，阵列由子阵1 和子阵2 两部分组成，共有N个阵元，N=N1+N2。子阵1 是一个密集的均匀线阵，共有N1个阵元，阵元间距为d=λ/2，λ是入射信号波长；子阵2 是一个稀疏线阵，共有N2个阵元，子阵2的0位置阵元和子阵1的N1-1位置阵元的间距也为d，子阵2前N2-1个阵元之间的间距为(2N1+2)d，最后两个阵元之间的间距为(2N1+1)d。

根据图1 所示的INA 阵元位置结构，可以定义SINAd为阵元位置，其中SINA是阵元位置集

阵列孔径LINA为

2.2 L型INA信号模型

如图2 所示的L 型INA，该阵列由位于x轴和z轴上的两个子阵组成，每一个子阵都是由2.1 提出的INA 组成，原点被两个子阵共用，总的阵元个数为2N-1，N=N1+N2。位于x轴和z轴上阵元的位置可以表示为(xid，0)和(0，zid)，其中xi，zi∈SINA，i=1，2，…，N。

假定有K个远场不相关的窄带信号s(t)以角度(θk，φk)，k=1，2，…，K入射到图2 所示的阵列，其中θk∈(0，π)是信号s(t)与x轴的夹角，称为方位角；φk∈(0，π)是信号s(t)与z轴的夹角，称为俯仰角。则对于x轴和z轴上第t个快拍的阵列输出可建模为：

其中s(t)=[s1(t)，s2(t)，…，sK(t)]T∈CK×1是信号源矢量，nx(t)=[nx1(t)，nx2(t)，…，nxN(t)]T∈CN×1和nz(t)=[nz1(t)，nz2(t)，…，nzN(t)]T∈CN×1分别是x轴和z轴上的噪声矢量，彼此是相互独立的均值为零方差为的高斯白噪声，并且与信号源不相关。Ax(θ)和Az(φ)是x轴和z轴上的阵列流形，x(t)和z(t)是x轴和z轴上第t个快拍的接收信号，可以表示为

ax(θk)和az(θk)，k=1，2，…，K分别是第k个信号在x轴和z轴上的导向矢量，表示为

3 阵列孔径虚拟扩展

3.1 VCAM算法

利用VCAM 算法可以生成2.1 节所提出的INA的虚拟扩展阵列，扩展后的虚拟阵列为均匀线阵，虚拟阵元之间的间距为d=λ/2。由N1+N2个物理阵元可虚拟出4N1N2+4N2-3个连续的虚拟阵元，即阵列的自由度由N1+N2提升到了4N1N2+4N2-3。下面介绍VCAM算法。

对于2.2 节中所提出的x轴上的阵列，在位置(xmd，0)和(xnd，0)，m，n=1，2，…，N上阵元的输出可以表示为xm(t)和xn(t)。则xm(t)和xn(t)的时间自相关函数可以表示为

当xn=0 时，即取(0，0)位置阵元作为参考阵元，向量化可以得

联合式子（12）和式子（13）可以得

用将Ts和Ns表示为伪采样周期和伪快照的数量，r(τ)的伪数据矩阵可以表示为

R的协方差矩阵表示为

同理对于z轴上的和差嵌套阵也可以通过VCAM算法生成其虚拟扩展阵列。

3.2 INA的虚拟扩展

利用VCAM 算法可以生成INA 的DSCA，阵元间 距d=λ/2，由N1+N2个物理阵元可虚拟出4N1N2+4N2-3 个连续的虚拟阵元，虚拟阵元的位置集SVINA在内连续

当N1=N2=4 时，即一共有8 个阵元的INA 的物理阵元结构，和其经过DSCA 虚拟后阵元结构，记为INA（DSCA）。如图3 所示，连续阵元的区间为［-38，38］，即阵列的自由度由8提升到了77。

为了说明本阵列的优越性，与同样是8 个阵元的嵌套阵NA、互质阵CA 和增强嵌套阵ANA 来进行比较。对NA 和CA 同时进行DCA 和DSCA 操作，对ANA 仅进行DCA 操作，为了方便，将其相对应的虚拟阵列分别记为NA（DCA）、NA（DSCA）、CA（DCA）、CA（DSCA）、ANA（DCA）。上述阵列的结果如图4 所示，可以发现同样是8 个阵元的情况下，下面三种阵列的最大连续范围（图（d））是［-28，28］，自由度为57，比本文提出的INA 阵列的自由度要少20，说明了本阵列的具有更高的自由度。

除此之外，同时对比图（a）和图（b）、或者图（c）和图（d），可以发现同一阵列经过DSCA 后的虚拟阵元数是要大于只进行DCA 的虚拟阵元数，说明了阵列进行DSCA后能够获得更高的自由度。

当NA（DSCA）、CA（DSCA）、ANA（DCA）以及INA（DSCA）的物理总阵元数N不变时，合理分配N1、N2的值可得到虚拟阵元数的最优表达式如表1所示。注意下面只列举了当N为偶数时的表达式，N为奇数的表达式与N为偶数的数量级是一致的。

观察表1可知，对于本文提出来的INA 的DSCA的自由度是O(N2)级别的，而NA 和CA 的DSCA 自由度仅为O(N2/2)级别，虽然ANA 的DCA 自由度为O(2N2/3)，但是还是小于O(N2)。表1 在数学表达式上说明了本阵INA具有更高的自由度。

表1 N阵元不同阵列的最大自由度Tab.1 The maximum DOF of different arrays of N elements

3.3 基于L型INA的二维DOA估计

通过VCAM 可以生成INA 的和差嵌套阵，其计算复杂度约为O(T(4N1N2+4N2-3)4)，其中T代表快拍数，N1、N2分别代表子阵1、2 的阵元数，式子（17）是向量化接收信号的协方差矩阵，矢量z等价于经过DSCA 虚拟化后的和差嵌套阵的接收信号。但是由于式子中矩阵p的秩为1，无法直接用MUSIC算法进行正确的DOA 估计，需要先进行秩的恢复，可以采用文献［9］提出的Toeplitz 矩阵重构的方法来恢复矩阵的秩。

通过对物理阵元位置的DSCA 操作，会得到一些冗余阵元，矢量z中的每个元素都对应一个虚拟阵元，因此z中也会存在冗余元素，对z中的元素进行去冗余，并按阵元位置从小到大排序，取中间连续的部分为得下式。

采用Toeplitz 矩阵重构的方法会损失一半的自由度，因此基于和差嵌套阵的Toeplitz 矩阵重构法最多可以估计-1个信号的入射角度。

基于L型INA的二维DOA估计步骤总结如下：

4 实验仿真

实验1 将对比INA（DSCA）与NA（DSCA）、CA（DSCA）、ANA（DCA）分别扩展至L 型阵上的二维DOA估计效果，L阵x轴和z轴的两个子阵相同，其物理阵元排列方式同3.2节，每个子阵的阵元设为N=8，则L阵总阵元数为2N-1=15，四种阵列采用Toeplitz矩阵重构法时的理论自由度分别为38、23、28、22。先对L阵的两个子阵分别采用MUSIC算法进行一维的DOA 估计，然后再采用PSCM 算法进行一维角度配对，进而完成L阵的二维DOA估计。

进行蒙特卡洛实验，引入均方根误差（RMSE）来评估上述四种阵的DOA 估计性能。定义RMSE为

其中W代表进行独立重复实验的次数，K代表总的入射信号个数，而分别代表第w次独立实验中对于第k信号的方位角θk和俯仰角φk的估计值。

设有K=3 个远场不相关的窄带信号其入射方向为(45°，70°)、(80°，40°)和(130°，90°)。

a）研究均方根误差与信噪比的关系，信噪比SNR 在［-5 dB，15 dB］内均匀变化，快拍数为500，谱峰搜索步长为0.01°，进行300 次蒙特卡洛实验。实验结果如图5 所示，随着SNR 增大，这四种阵列的DOA 估计性能有明显的提高，且INA（DSCA）相比其他三个虚拟阵列一直有着更低RMSE 这是因为INA（DSCA）具有更多的连续阵元，更高的自由度。

b）研究均方根误差与快拍数的关系，快拍在[100，1000]内均匀变化，信噪比SNR=5 dB，谱峰搜索步长为0.01°，进行300次蒙特卡洛实验。实验结果如图6 所示，随着快拍数增加，这四种阵列的DOA估计性能有明显的提高，且INA（DSCA）相比其他三个虚拟阵列一直有着更低RMSE，这是因为INA（DSCA）具有更高的自由度。

实验2 研究INA（DSCA）的L阵可估计信源数，并与NA（DSCA）、CA（DSCA）和ANA（DCA）做对比，其他L 阵的结构同实验1。设有11个远场不相关的窄带信号其入射方向为（30°，75°）、（45°，80°）（50°，140°）、（60°，55°）、（65°，85°）、（70°，120°）、（80°，60°）、（90°，90°）、（120°，95°）、（135°，35°）和（140°，50°），信噪比SNR=10 dB，快拍数为500，谱峰搜索步长为0.01°。进行DOA 估计的结果如图7 所示，结果表明每个子阵是8阵元INA（DSCA）组成的L阵可以准确估计出11个入射信号的二维角度，突破了物理阵元个数的限制。且对比NA（DSCA）、CA（DSCA）和ANA（DCA）有更高的估计的估计精度，这是因为INA（DSCA）具有更高的自由度。8阵元本阵列理论自由度38，这里只估计11 个入射信号的原因是，PSCM 这类的匹配算法，尽管是针对稀疏阵列的角度配对，可以减少计算复杂度，在信源数不多时可以实现角度的准确配对，但是当信源数远远大于阵元数会出现匹配失败的现象。

5 结论

本文针对传统L型均匀阵列二维波达方向估计中可估计信源数目受限于阵元数、分辨率低等问题，提出了一种L 型INA 的阵列结构。INA 由一个密集线阵和一个稀疏线阵两部分组成，通过VCAM算法对INA 的物理阵元位置进行求和求差（DSCA）的虚拟扩展操作，由N1+N2个物理阵元可得到4N1N2+4N2-3 的自由度，相比NA（DSCA）、CA（DSCA）和ANA（DCA）具有更高的自由度。并与这三种阵列进行二维DOA 估计的实验仿真对比，结果显示此阵具有更高的估计精度和分辨率。以及利用此阵去估计超过子阵元数的信源数，实验表明此阵具有高估计精度和能突破物理阵元限制的结论。此外，文中采用的VSCM算法，是将二阶统计量转换为四阶统计量来进行运算，即阵列自由度和估计精度提高是以增加一定的计算复杂度为代价的。