南京航空航天大学学报  2017, Vol. 49 Issue (5): 635-644   PDF    
互质阵中空间谱估计研究进展
张小飞, 林新平, 郑旺, 翟会     
南京航空航天大学电子信息工程学院,南京,210016
摘要: 阵列信号处理技术是信号处理领域一个重要分支,广泛应用于雷达、声呐、水电天文等领域。空间谱估计是阵列信号处理中一个主要研究热点,为了避免角度模糊,要求阵元间距要小于等于载波波长的一半。互质阵突破半波长的限制,其阵元间距一般大于半波长。因此,互质阵能够在阵元数固定的情况下,获得更大的阵列孔径。互质阵空间谱估计算法可以利用互质特性辨识目标源,且能够获得更高的测向精度与分辨率。互质阵空间谱估计研究逐渐成为当今阵列信号处理领域的热点。本文总结了互质线阵与互质面阵下的空间谱估计的研究进展, 分析了互质阵列相对于普通均匀阵列在阵列孔径扩展以及提高空间自由方面的优势。仿真结果展现了互质阵下高效的波达方向估计性能。
关键词: 互质线阵     互质面阵     空间谱估计     波达方向估计    
Research Progress on Spatial Spectrum Estimation Based on Coprime Array
ZHANG Xiaofei, LIN Xinping, ZHENG Wang, ZHAI Hui     
College of Electronic and Information Engineering, Nanjing University of Aeronautics & Astronautics, Nanjing, 210016, China
Abstract: Array signal processing technology is an important part in the field of signal processing, which is widely used in radar, sonar and hydropower astronomy. Spatial spectrum estimation is a major research hotspot in array signal processing and is confined to the arrays with inter-element spacing smaller than half wavelength of the carrier wave to avoid ambiguous angles. Generally, the inter-element spacing of coprime arrays is larger than half wavelength. As a result, with fixed number of sensors, the coprime arrays can obtain extended array aperture and hence the corresponding spectrum estimation methods can achieve better direction of arrival estimation performance and resolution, where the targets can be detected uniquely with coprime property. Therefore, the investigations of spatial spectrum estimation with coprime arrays have aroused considerable attentions in array signal processing. In this paper, we summarize the existing researches based on coprime linear array and coprime planar array and analysis the advantages of the coprime array relative to the uniform array in expanding array aperture and increasing spatial degree of freedom. Simulation results validate the effectiveness and superiority of the methods with coprime arrays.
Key words: coprime linear array     coprime planar array     spatial spectrum estimation     direction of arrival estima tion    

阵列信号处理的基本原理是利用收集到的有用的空间信号,对其进行信号处理,同时对噪声信号进行抑制,从而得到源信号的特征及其所包含的有用的信息[1]。空间谱估计[2-4],也称为波达方向(Direction of arriva l, DOA)估计,主要是为了估计空域中接收到的信号的来波方向, 在阵列信号处理领域中是一个主要研究方向。

经典的超分辨率DOA估计方法主要有多重信号分类(Multiple signal classification, MUSIC)算法[5]、借助旋转不变性进行信号参数估计(Estimation of signal parameters v ia rotational invariance technique, ESPRIT)算法[6],这些都是基于特征结构的子空间类方法。子空间类方法都是最初均是针对传统阵列提出的,并且要求阵元间距小于载波的半波长,以此避免出现角度模糊。由于将线性空间概念引入到DOA估计中,子空间方法实现了超分辨率。子空间类方法的DOA估计分辨率与传感器阵列有效孔径成正比,即在相同阵元数目的情况下,随着阵元间距的增大,DOA估计算法可以获得更高的测向精度和分辨率[7-9]。实际中,无论是军用雷达还是民用通信[10-11],都对测向精度与分辨率有很高的要求。传统阵列需要通过增加阵元数目来得到大的阵列孔径,但是这会增加系统成本和算法的复杂度。此外,传统阵列的DOA估计算法中还存在着欠定问题,即当待估计信源数目大于物理阵元数目时,上述子空间类的DOA估计算法将会失效。因此,传统阵列由于阵元间距的限制,在许多航空航天及军事民用领域已经无法满足需求。

文献[12]中提出了互质线阵的概念。互质线阵由阵元数分别为MN的两个均匀线阵组成,其中,MN为互质数,子阵的阵元间距分别为/2和/2,λ为载波波长。特别地,文献[12]已经证明了通过阵元总数为M+N-1的互质阵可以获得O(MN)的自由度。与最小冗余阵相比,互质阵可以得到阵元位置与自由度的闭式解,阵列结构简单,无需复杂的搜索即可获得优化的阵列结构。近年来,许多学者在互质阵列上作了大量科学研究。文献[13]中提出了一种基于互质特性的联合MUSIC估计算法。文章利用子阵阵元数的互质特性,通过对比两个子阵的DOA估计结果,证明了DOA估计值的存在性与唯一性,从而消除角度模糊的问题。该方法虽然在实现上较为简单,但是将两个子阵单独运作,可辨识的目标数减少至少一半,而且阵列间的互信息损失也会一定程度上削弱DOA估计的性能。文献[14]中,作者提出了一种基于协方差矩阵矢量化的DOA估计算法。该算法通过对接收信号协方差矩阵进行矢量化操作进行数据重构,从而得到一个通过虚拟阵列接收到的单快拍信号。特别地,该虚拟阵列的孔径与自由度都大于实际物理阵列,从而实现孔径与自由度的拓展。为了解决重构信号的相干性,作者引入了空间平滑技术。为了解决空间平滑技术带来的自由度损失的问题,文献[15~17]引入压缩感知技术,充分利用互质阵列提供的最大孔径与自由度。

事实上,空间目标源定位是一个二维DOA(仰角与方位角)估计的问题,国内外学者基于传统的紧凑阵列提出了二维(Two-dimensional, 2D) MUSIC算法、2D ESPRIT算法、二维传播算子方法(2D propagator method, 2D PM)及这些算法的众多改进算法。在互质阵的概念提出后,国内外学者提出了互质L型阵列、互质十字阵和互质面阵列,并给出了一系列无模糊且高效的2D DOA估计算法。文献[18]研究了互质L型阵列下的2D DOA估计问题,利用两个L型子阵的互质特性,仰角和方位角的估计值可以唯一确定。基于传统L型阵列,文献[19]提出了一种互质十字阵,并提出了一种改进的高效子空间(Improved computationally efficient subspace algorithm, ICESA)算法。利用接收数据的四阶累计量得到数据的协方差矩阵,该算法比2D ESPRIT算法复杂度低、角度估计性能更好。文献[20]中,作者提出了一种互质面阵,并利用2D MUSIC算法进行角度估计。由于2D MUS IC算法需要复杂度很高的二维谱峰搜索,作者利用模糊DOA值与真实DOA值之间的正弦线性关系,提出了一种局部谱峰搜索算法,大大地降低了算法的复杂度。在文献[20]的基础上,文献[21]通过降维转换,将二维局部谱峰搜索降低至一维局部谱峰算法,进一步降低了算法的复杂度,且不影响DOA估计性能。由于上述算法都是将互质阵的两个子阵独立运算,自由度被减少了至少一半,文献[22]利用和阵列与差阵列,得到一个大于实际阵列的虚拟阵列,提高了自由度与阵列孔径。作者引入稀疏恢复的方法,充分地利用互质面阵所能提供的最大自由度与阵列孔径,角度估计精度得到了极大的提升。

论文第一节主要介绍了互质阵列提出的背景以及发展互质阵列DOA估计算法的重要意义。第二节介绍了互质线阵模型及基于互质线阵的DOA估计算法。第三节介绍了互质面阵模型及基于互质面阵的2D DO A估计算。第四节为总结与期望,对整篇论文的概括性总结以及对未来互质阵列发展的展望。

符号表示:本文中用(·)T表示矩阵转置,(·)H表示矩阵共轭转置,大写字母X表示矩阵,小写字母x(·)表示矢量,Xi, j表示矩阵的第(i, j)个元素,Ip表示为P×P的单位矩阵,⊗表示Kronecker积,⊙表示Khatri-Rao积,diag(x)表示以x中元素构成的对角矩阵。

1 互质线阵空间谱估计

目前关于互质阵列空间谱估计方法主要有两种,一是利用联合互质阵列子阵列的思想进行空间谱估计,二是对互质阵列得到的协方差矩阵进行矢量化后构造虚拟阵列的空间谱估计方法。

1.1 互质线阵拓扑结构与数据模型

互质线阵的基础拓扑结构如图 1所示[14],互质线阵主要由两个阵元间距大于半波长的均匀线阵所构成,记为子阵一和子阵二,阵元数分别为M1M2,其阵元间距为d1=M2λ/2,d2=M1λ/2,其中M1M2互为质数。

图 1 互质线阵基础结构拓扑图 Figure 1 Topological structure of coprime linear array

互质阵的两个子阵均以原点为参考点,阵元位置可记为如下集合

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{L}}s = \left\{ {\left( {{m_1}{d_1}\left| {{m_1} = 0,1,2, \cdots ,{M_1} - 1} \right.} \right.} \right\} \cup }\\ {\left\{ {{m_2}{d_2}\left| {{m_2} = 0,1,2, \cdots ,{M_2} - 1} \right.} \right\}} \end{array} $ (1)

假设空间远场有K个不相关窄带信号入射到上述互质线阵上,其入射角度分别为[θ1, θ2, …, θK]时,阵列接收信号可以表示为

$ \mathit{\boldsymbol{x}}\left( t \right) = \mathit{\boldsymbol{As}}\left( t \right) + \mathit{\boldsymbol{n}}\left( t \right) $ (2)

式中:s(t)为信源矢量,n(t)为噪声矢量,A为方向矩阵

$ \mathit{\boldsymbol{A}} = \left[ {\mathit{\boldsymbol{a}}\left( {{\theta _1}} \right),\mathit{\boldsymbol{a}}\left( {{\theta _2}} \right), \cdots ,\mathit{\boldsymbol{a}}\left( {{\theta _K}} \right)} \right] $ (3)

式中:a(θk)= $ \left[ {1,{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda } d \cdot \sin \left( {{\theta _k}} \right)}}, \cdots ,} {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }\left( {M - 1} \right) d \cdot \sin \left( {{\theta _k}} \right)}}} \right]^{\rm{T}}$为方向矢量,dλ分别为阵元间距和载波波长。子阵一和子阵二的方向矢量分别为

$ \begin{array}{l} {\mathit{\boldsymbol{a}}_1}\left( {{\theta _k}} \right) = \left[ {1,{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }{M_2} \cdot d \cdot \sin \left( {{\theta _k}} \right)}}, \cdots ,} \right.\\ \;\;\;\;\;\;\;\;\;\;\;{\left. {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }\left( {{M_1} - 1} \right) \cdot {M_2} \cdot d \cdot \sin \left( {{\theta _k}} \right)}}} \right]^{\rm{T}}} \end{array} $ (4)
$ \begin{array}{l} {\mathit{\boldsymbol{a}}_2}\left( {{\theta _k}} \right) = \left[ {1,{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }{M_1} \cdot d \cdot \sin \left( {{\theta _k}} \right)}}, \cdots ,} \right.\\ \;\;\;\;\;\;\;\;\;\;\;{\left. {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }\left( {{M_2} - 1} \right) \cdot {M_1} \cdot d \cdot \sin \left( {{\theta _k}} \right)}}} \right]^{\rm{T}}} \end{array} $ (5)

可以从式(4,5)中得出,用整个阵列进行DOA估计时,方向矢量为

$ \mathit{\boldsymbol{a}}\left( {{\theta _k}} \right) = {\left[ {\mathit{\boldsymbol{a}}_1^{\rm{T}}\left( {{\theta _k}} \right),\mathit{\boldsymbol{a}}_{21}^{\rm{T}}\left( {{\theta _k}} \right)} \right]^{\rm{T}}} $ (6)

式中:a21T (θk)为a2T(θk)去除第一行元素后构成的新方向矢量。

1.2 DOA估计算法 1.2.1 解模糊方法

以子空间类算法中的ESPRIT算法为例,解模糊方法通过联合互质线阵中的两个均匀子阵DOA估计结果获得测向角,假定子阵i (i=1, 2)得到接收信号xi(t),协方差矩阵Rxi = E[xi(txiH(t)],对协方差矩阵特征分解后得到信号子空间Es, i,再利用信号子空间的前Mi-1行和后Mi-1行得到两个新矩阵Es1, iEs2, i。由旋转不变性得

$ \mathit{\boldsymbol{E}}_{s1,i}^ + {\mathit{\boldsymbol{E}}_{s2,i}} = {\mathit{\boldsymbol{T}}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} T}} $ (7)

式中:T为满秩矩阵,Φ=diag(${{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda } d_i \sin \left( {{\theta _1}} \right)}}, \cdots , {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda } d_i \sin \left( {{\theta _K}} \right)}}} $),di为子阵i阵元间距,Φ中包含着DOA信息。通过对式(7)左边特征分解后就能得到DOA估计值。然而,由于两个子阵的阵元间距大于半波长,两个子阵得到的DOA估计值存在角度模糊问题。

文献[13]指出,DOA模糊值与真实值之间存在如下关系

$ \mathit{\boldsymbol{a}}\left( \theta \right) = \mathit{\boldsymbol{a}}\left( {{\theta ^a}} \right) $ (8)

式中:θ为真实值,θa为对应的模糊值。由上可知

$ 2{\rm{ \mathit{ π} }}{d_i}\sin \theta /\lambda - 2{\rm{ \mathit{ π} }}{d_i}\sin {\theta ^a}/\lambda = 2k{\rm{ \mathit{ π} }}\;\;\;\;\;i = {\rm{1,2}} $ (9)

由式(9)可知,当子阵中阵元间距为di=Miλ/2时,因为一维测向角度范围是(-90 °, 90°),根据|sinθ-sinθa|<2,可以从式(9)中推导出k的取值范围是k∈[-(Mi-1), …, (Mi-1)]。又因为正弦函数的范围,所以k实际只有Mi个取值。当Mi等于1时,k只能取值为0,此时估计结果只有一个,也就是传统阵列形式下的DOA估计结果。当Mi的值大于1时,互质线阵的子阵会产生Mi个模糊值,其中包含对应的真实角度估计值。

另外,文献[13]中给出利用两个子阵得到的相同的角度估计值消除模糊的证明,如下:

假设除了θ为两子阵中相同值外,还有一个值θ′也为两子阵得到的角度估计的相同值。根据式(9)可知阵元数为M1的子阵1阵元间距为d1=M2λ/2有关系式

$ \sin \left( {{\theta _k}} \right) - \sin \left( {{{\theta '}_k}} \right) = 2{k_1}/{M_2} $ (10)

式中:k1=-(M2-1), …, -1, 1, …, (M2-1)。

同理,阵元数为M2的子阵2阵元间距为d2=M1λ/2有关系式

$ \sin \left( {{\theta _k}} \right) - \sin \left( {{{\theta '}_k}} \right) = 2{k_2}/{M_1} $ (11)

式中:k2=-(M1-1), …, -1, 1, …, (M1-1)

根据式(10,11)可得

$ \frac{{2{k_1}}}{{{M_2}}} = \frac{{2{k_2}}}{{{M_1}}} $ (12)

由于M1, M2互为质数,所以不存在k1, k2使式(12)成立。

所以,通过搜索互质线阵中两个子阵得到的相同的角度估计值就能得到唯一的DOA估计结果。

1.2.2 虚拟化阵元方法

上述提出的解模糊方法虽然在DOA估计的子空间类算法中的算法复杂度和算法性能上均要优于均匀线阵,但是以互质子阵的方式估计角度牺牲了空间自由度,即可探测的信源数极大地减少了。文献[14]中提出的虚拟化阵元方法则能明显地提高阵列算法的空间自由度,且不存在角度模糊问题。

假设互质线阵的子阵一和子阵二阵元数分别为M1M2,其中M1M2互为质数。互质线阵得到的信号可以表示为以下形式

$ \mathit{\boldsymbol{x}}\left( t \right) = {\left[ {\mathit{\boldsymbol{x}}_1^{\rm{T}}\left( t \right),\mathit{\boldsymbol{x}}_2^{\rm{T}}\left( t \right)} \right]^{\rm{T}}} $ (13)

式中:x1(t)和x2(t)分别为子阵一和子阵二的接收信号矩阵,互质线阵方向矢量表示为a(θk) = [aM1T(θk), aM2T(θk)]T,接收信号的协方差矩阵为Rxx=E[x(t)xH(t)],对协方差矩阵矢量化后可得

$ \mathit{\boldsymbol{z}} = {\rm{vec}}\left( {{\mathit{\boldsymbol{R}}_{xx}}} \right) = \mathit{\boldsymbol{B}}\left( {{\theta _1},{\theta _2}, \cdots ,{\theta _K}} \right)\mathit{\boldsymbol{P}} + \sigma _n^2\mathit{\boldsymbol{H}} $ (14)

式中:P = [σ12, σ22, …, σK2]Tσk2表示为第k个信号源的信号功率,P可以看作是单快拍的新的信号矩阵,H=vec (I),新的方向矩阵可以表示为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{B}}\left( {{\theta _1},{\theta _2}, \cdots ,{\theta _K}} \right) = \left[ {{\mathit{\boldsymbol{a}}^ * }\left( {{\theta _1}} \right) \otimes \mathit{\boldsymbol{a}}\left( {{\theta _1}} \right),{\mathit{\boldsymbol{a}}^ * }\left( {{\theta _2}} \right) \otimes } \right.}\\ {\left. {\mathit{\boldsymbol{a}}\left( {{\theta _2}} \right), \cdots ,{\mathit{\boldsymbol{a}}^ * }\left( {{\theta _K}} \right) \otimes \mathit{\boldsymbol{a}}\left( {{\theta _K}} \right)} \right]} \end{array} $ (15)

从式(15)中可以得出新方向矢量包含的相位信息有以下3种

$ \begin{array}{l} {{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }d\sin \left( {{\theta _k}} \right)\left( {{M_1}{m_2} - {M_2}{m_1}} \right)}}\\ \;\;\;0 \le {m_2} \le {M_2} - 1,\;\;\;\;0 \le {m_1} \le {M_1} - 1 \end{array} $ (16)
$ \begin{array}{l} {{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }d\sin \left( {{\theta _k}} \right)\left( {{M_1}{m_{21}} - {M_1}{m_{22}}} \right)}}\\ \;\;\;0 \le {m_{21}},\;\;\;\;{m_{22}} \le {M_2} - 1 \end{array} $ (17)
$ \begin{array}{l} {{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }d\sin \left( {{\theta _k}} \right)\left( {{M_2}{m_{11}} - {M_2}{m_{12}}} \right)}}\\ \;\;\;0 \le {m_{11}},{m_{12}} \le {M_1} - 1 \end{array} $ (18)

文献[14]指出,上述公式中的相位包含的其实是一种差分集合形式

$ \begin{array}{l} {S_{{M_1}{M_2}}} = \left\{ {\left( {{M_1}{m_2} - {M_2}{m_1}} \right),0 \le {m_2} \le {M_2} - 1,} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {0 \le {m_1} \le {M_1} - 1} \right\} \end{array} $ (19)

集合中包含的距离原点范围包含于[-M2(M1-1), …, M1(M2-1)]中,其假设单位阵元间距为半波长。需要指出的是该虚拟阵列的部分为连续阵列,且连续部分并没有确定的闭式解。文献[14]在此基础之上提出增广互质线阵结构,其子阵一和子阵二阵元数分别为M2和2M1,如图 2所示,其中M1M2互为质数,且假设M1M2

图 2 增广互质线阵拓扑结构 Figure 2 Topological structure of a ugmented coprime linear array

文献[14, 23]指出在增广互质线阵的差分集合中,连续的阵元部分为[-M1M2, …, M1M2]。这意味着阵元数为2M1+M2-1的互质阵列可以获得2M1M2+1的自由度。

文献[14]利用空间平滑技术与MUSIC算法进行DOA估计。假设z1z消除重复行并且按距离位置排列好后的数据矩阵,z1i是包含有M1M2+1个阵元的第i个子阵列的接收数据矩阵,构造协方差矩阵为

$ {\mathit{\boldsymbol{R}}_i} = {\mathit{\boldsymbol{z}}_{1i}}\mathit{\boldsymbol{z}}_{1i}^{\rm{H}} $ (20)

对所有子阵列协方差矩阵取平均可得

$ {\mathit{\boldsymbol{R}}_{ss}} = \frac{1}{{{M_1}{M_2} + 1}}\sum\limits_{i = 1}^{{M_1}{M_2} + 1} {{\mathit{\boldsymbol{R}}_i}} $ (21)

可以从式(21)中得出Rss= $\boldsymbol{\hat R} ^2$, 其中

$ \mathit{\boldsymbol{\hat R}} = \frac{1}{{\sqrt {{M_1}{M_2} + 1} }}\left( {{\mathit{\boldsymbol{A}}_{11}}\mathit{\boldsymbol{ \boldsymbol{\varLambda} A}}_{11}^{\rm{H}} + \sigma _n^2{\mathit{\boldsymbol{I}}_{{M_1}{M_2}}}} \right) $ (22)

式中:Λ =diag(σ12, σ22, …, σK2),A11表示空间平滑算法中第一个参考子阵列的方向矩阵。上述表达式符合基本子空间类DOA估计算法的协方差分解形式,故能利用它得到正确的DOA估计。

以ESPRIT算法为例,通过对式(22)中得到的$\boldsymbol{\hat R} $特征分解后得到信号子空间,再利用信号子空间的前M1M2行和后M1M2行得到两个新矩阵Es1Es2,构造类似于式(7)中的方程式Φ=(Es1)+Es2,通过对Φ特征分解后得到的特征值就能得到正确的DOA估计值。

虚拟化阵元方法相对解模糊方法有效的提高了空间自由度,但是因为引进了空间平滑算法,同时造成了空间自由度的损失。

下面给出子空间类算法ESPRIT算法下几种方法的DOA估计性能比较。本文中所有的仿真均采用蒙特卡洛仿真来评估算法性能。假设空间远场信号源入射角为θ1=10°, θ2=20°。图 3的仿真条件为普通均匀线阵阵元数M= 11, 互质线阵解模糊方法阵元数M1=6, M2=5, 虚拟化阵元方法互质线阵阵元数M1=3, M2=5。保证总阵元数相同,快拍数L=300条件下的几种方法随着信噪比条件下的DOA估计性能图。

图 3 互质线阵中不同方法下的DOA估计性能 Figure 3 DOA estimation perfor mance under different methods in coprime linear array

可以从图 3中看出,在一般情况下,解模糊方法的DOA估计性能优于虚拟化阵元方法,优于均匀线阵方法,但是随着信噪比的增高,虚拟化阵元方法的DOA估计性能变差。

图 4给出了在图 3相同的阵元数情况下,在信噪比SNR=0 dB条件下,几种方法DOA估计性能随着快拍数的变化情况。随着快拍数的增加,几种方法下的DOA估计性能提升了。其中,相对于传统阵列,互质阵列可以获得更好的角度估计性能。

图 4 互质线阵中不同快拍数下算法性能 Figure 4 Algorithm performance under different snapshots in coprime line ar array

综上所述,相对于均匀线阵,互质线阵不仅在阵列天线孔径与空间自由度上有优势,而且在算法性能上,互质线阵也要优于均匀线阵。

2 互质面阵空间谱估计 2.1 互质面阵拓扑结构与数据模型

文献[20]中描述的互质面阵结构拓扑图如图 5所示,互质面阵的阵列排布较为复杂,和传统均匀面阵相比,互质面阵可以分成两个子阵,记为子阵1和子阵2。其中,子阵1是由阵元数为M1×M1的均匀面阵构成,子阵2由阵元数为M2×M2的均匀面阵构成。其中子阵1在X轴和Y轴方向上的阵元间距都为d1=M2λ/2,子阵2在X轴和Y轴方向上的阵元间距都为d2=M1λ/2。

图 5 互质面阵结构拓扑图 Figure 5 Topological structure of coprime planar array

考虑上述的互质面阵,子阵阵元只在原点位置重合,故互质面阵的阵元总数为T=M12+M22-1。阵元在坐标轴的位置可以表示为如下集合

$ \begin{array}{*{20}{c}} {{L_s} = \left\{ {\left( {m{d_1},n{d_1}} \right)\left| {0 \le m,n \le {M_1} - 1} \right.} \right\} \cup }\\ {\;\;\;\;\left\{ {\left( {p{d_2},q{d_2}} \right)\left| {0 \le p,q \le {M_2} - 1} \right.} \right\}} \end{array} $ (23)

假定空间有K个不相关的窄带信源入射到上述面阵中,并且θk, ϕk表示第k个信源的仰角和方位角。

考虑互质面阵的子阵i(i=1, 2),其在x轴和y轴上阵元的方向矢量可以分别表示为

$ \begin{array}{l} {\mathit{\boldsymbol{a}}_{ix}}\left( {{\theta _k},{\phi _k}} \right) = \left[ {\begin{array}{*{20}{c}} 1\\ {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }{d_i}\sin {\theta _k}\cos {\phi _k}}}}\\ \vdots \\ {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }\left( {{M_i} - 1} \right){d_i}\sin {\theta _k}\cos {\phi _k}}}} \end{array}} \right]\\ {\mathit{\boldsymbol{a}}_{iy}}\left( {{\theta _k},{\phi _k}} \right) = \left[ {\begin{array}{*{20}{c}} 1\\ {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }{d_i}\sin {\theta _k}\sin {\phi _k}}}}\\ \vdots \\ {{{\rm{e}}^{ - {\rm{j}}\frac{{2{\rm{ \mathit{ π} }}}}{\lambda }\left( {{M_i} - 1} \right){d_i}\sin {\theta _k}\sin {\phi _k}}}} \end{array}} \right] \end{array} $ (24)

那么从式(24)中可以得到子阵ix轴和y轴上的方向矩阵分别表示为Aix=[aix(θ1, ϕ1), aix(θ2, ϕ2), …, aix(θk, ϕk)],Aiy=[aiy(θ1, ϕ1), aiy(θ2, ϕ2), …, aiy(θK, ϕK)]。

则在x轴上的阵元接收信号可以表示为

$ {\mathit{\boldsymbol{x}}_{i1}}\left( t \right) = {\mathit{\boldsymbol{A}}_{ix}}\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{n}}_{i1}}\left( t \right) $ (25)

式中:ni1(t)为子阵ix轴上阵元的噪声。s(t)∈CK×1为信源矢量。平行于x轴的第m个分阵列接收信号可以表示为

$ {\mathit{\boldsymbol{x}}_{im}}\left( t \right) = {\mathit{\boldsymbol{A}}_{ix}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{iy}^{m - 1}\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{n}}_{im}}\left( t \right) $ (26)

式中:Φiy = diag(${{\rm{e}}^{ - {\rm{j}}{2{\rm{ \mathit{ π} }}} d_i \sin {{\theta _1}} \sin \phi_1/\lambda}}, \cdots , {{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathit{ π} }} d_i \sin {{\theta _K}}\sin\phi_K/\lambda}}} $),nim(t)为第m个分阵列的噪声。则可得整个子阵i的接收信号xi(t)为

$ {\mathit{\boldsymbol{x}}_i}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_{i1}}\left( t \right)}\\ {{\mathit{\boldsymbol{x}}_{i2}}\left( t \right)}\\ \vdots \\ {{\mathit{\boldsymbol{x}}_{i{M_i}}}\left( t \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_{ix}}}\\ {{\mathit{\boldsymbol{A}}_{ix}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{iy}}}\\ \vdots \\ {{\mathit{\boldsymbol{A}}_{ix}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{iy}^{{M_i} - 1}} \end{array}} \right]\mathit{\boldsymbol{s}}\left( t \right) + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{n}}_{i1}}\left( t \right)}\\ {{\mathit{\boldsymbol{n}}_{i2}}\left( t \right)}\\ \vdots \\ {{\mathit{\boldsymbol{n}}_{i{M_i}}}\left( t \right)} \end{array}} \right] $ (27)

式(27)中的信号也可以表示为

$ \begin{array}{l} {\mathit{\boldsymbol{x}}_i}\left( t \right) = \left[ {{\mathit{\boldsymbol{A}}_{iy}} \odot {\mathit{\boldsymbol{A}}_{ix}}} \right]\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{n}}_i}\left( t \right) = \left[ {{\mathit{\boldsymbol{a}}_{iy}}\left( {{\theta _1},{\phi _1}} \right) \otimes } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {{\mathit{\boldsymbol{a}}_{ix}}\left( {{\theta _1},{\phi _1}} \right), \cdots ,{\mathit{\boldsymbol{a}}_{iy}}\left( {{\theta _K},{\phi _K}} \right) \otimes {\mathit{\boldsymbol{a}}_{ix}}\left( {{\theta _K},{\phi _K}} \right)} \right] \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{n}}_i}\left( t \right) = {\mathit{\boldsymbol{A}}_i}\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{n}}_i}\left( t \right) \end{array} $ (28)

式中:Ai=AiyAixni(t)为整个子面阵i上的噪声。

2.2 DOA估计方法

与互质线阵一样的是,互质面阵的DOA估计方法也包括两种:一种是解模糊方法,另一种是虚拟化阵元方法。

2.2.1 解模糊方法

以子空间类算法中的ESPRIT算法为例,互质面阵中解模糊方法也是通过联合两个子阵估计结果获得DOA估计。通过子阵i得到的式(28)中的接收信号构造协方差矩阵Rxi = E[xi (t)xiH (t)],特征分解后得到信号子空间Esp=(AiyAix)TT为满秩矩阵,通过信号子空间Esp的前Mi(Mi-1)行和后Mi(Mi-1)行得到Esp1Esp2,构造类似于式(7)中的方程Φp = Esp1+Esp2,特征分解后通过特征值就得到y轴方向上的仰角与方位角的估计值。

同样的,信号子空间也可以表示为Esp′=(AixAiy)T,利用Φp特征分解后的特征向量得到,同样的方法处理Esp′就能得到x轴方向上的仰角与方位角的估计值。此时应当注意的是,Esp′是通过Esp右乘$ \boldsymbol{\hat T}$后变换得到的,这样能实现仰角与方位角的自动匹配。两个子阵有两组估计结果,此时估计值中因为阵元间距稀疏性而出现角度模糊问题。

文献[20]中详细推导了互质面阵的解模糊方法的DOA估计理论。与互质线阵一致,模糊产生是由于阵元间距大于半波长引起的。假设有空间角度为(θk, ϕk)(k=1, 2, …, K)的K个非相干窄带信源入射到互质面阵中,θk, ϕk分别代表第k个信源的仰角和方位角。考虑面阵中的子阵列i,其阵元间距假设为di=Miλ/2,那么子阵ix轴方向和y轴方向上的相邻两个阵元所接收信号的相位差分别可表示为

$ \Delta {\varphi _x} + 2{k_x}{\rm{ \mathit{ π} }} = 2{\rm{ \mathit{ π} }}/\lambda {d_i}\sin {\phi _k}\cos {\theta _k} $ (29)
$ \Delta {\varphi _y} + 2{k_y}{\rm{ \mathit{ π} }} = 2{\rm{ \mathit{ π} }}/\lambda {d_i}\sin {\phi _k}\sin {\theta _k} $ (30)

式中:Δφxx轴方向上相邻阵元的相位差,Δφyy轴方向上相邻阵元间的相位差。kxky为整数,又因为θk∈[0, π],ϕk∈[0, π/2],所以有-1≤sinϕkcosθk≤1,0≤sinϕksinθk≤1。则kxky的取值范围分别为

$ {k_x} \in \left[ {\frac{{ - {d_i}}}{\lambda } - \frac{{\Delta {\varphi _x}}}{{2{\rm{ \mathit{ π} }}}},\frac{{{d_i}}}{\lambda } - \frac{{\Delta {\varphi _x}}}{{2{\rm{ \mathit{ π} }}}}} \right] $ (31)
$ {k_y} \in \left[ {\frac{{ - d}}{\lambda } - \frac{{\Delta {\varphi _y}}}{{2{\rm{ \mathit{ π} }}}},\frac{d}{\lambda } - \frac{{\Delta {\varphi _y}}}{{2{\rm{ \mathit{ π} }}}}} \right] $ (32)

同时,kxky还要受0≤(sinϕkcosθk)2+(sinϕksinθk)2≤1的约束。当dλ/2时,kxky只能取零。此时角度估计只有一个解,不存在角度模糊问题。当di=Miλ/2,Mi为大于1的整数时,kxky的取值个数分别为MiMi/2。即使考虑到(kx, ky)配对时,仍然有多个DOA角度满足式(31,32),但是其中只有一个为真实角度,其他为模糊角。

文献[20]中指出,类似于互质线性阵列角度模糊消除原理,互质面阵的角度模糊消除也是通过联合两个子阵间的角度关系,找到相同的真实估计值来达到消除角度模糊的目的。

根据式(29,30),当阵元间距为di=Mjλ/2时,真实DOA估计角度(θp, ϕp)和子阵i的模糊角度(θi, a, ϕi, a)之间满足以下关系

$ \sin {\phi _p}\cos {\theta _p} - \sin {\phi _{i,a}}\cos {\theta _{i,a}} = \frac{{2{k_{i,x}}}}{{{M_j}}} $ (33)
$ \sin {\phi _p}\cos {\theta _p} - \sin {\phi _{i,a}}\sin {\theta _{i,a}} = \frac{{2{k_{i,y}}}}{{{M_j}}} $ (34)

式中:i, j=1, 2;ijki, xki, y都是整数,且取值范围为(-Mj, Mj)和(-Mj/2, Mj/2)。

与线阵一致,考虑只有单个信源入射到互质面阵,可知两个子阵中都能得到正确DOA估计值,此时存在模糊,因此需证明两个子阵中有且只有一个模糊值是相同的。文献[20]给出证明如下:

假设子阵1和子阵2存在两个相同的估计结果($\hat \theta_1\;,\hat \phi_1 $)和($\hat \theta_2\;,\hat \phi_2 $),对子阵1而言,根据式(33,34)有

$ \sin {{\hat \phi }_1}\cos {{\hat \theta }_1} - \sin {{\hat \phi }_2}\cos {{\hat \theta }_2} = \frac{{2{k_{1,x}}}}{{{M_2}}} $ (35)
$ \sin {{\hat \phi }_1}\sin {{\hat \theta }_1} - \sin {{\hat \phi }_2}\sin {{\hat \theta }_2} = \frac{{2{k_{1,y}}}}{{{M_2}}} $ (36)

式中:k1, xk1, y为整数,且其取值范围分别为(-M2, M2)和(-M2/2, M2/2)。同理对于子阵2有

$ \sin {{\hat \phi }_1}\cos {{\hat \theta }_1} - \sin {{\hat \phi }_2}\cos {{\hat \theta }_2} = \frac{{2{k_{2,x}}}}{{{M_1}}} $ (37)
$ \sin {{\hat \phi }_1}\sin {{\hat \theta }_1} - \sin {{\hat \phi }_2}\sin {{\hat \theta }_2} = \frac{{2{k_{2,y}}}}{{{M_1}}} $ (38)

式中:k2, xk2, y为整数,且其取值范围分别为(-M1, M1)和(-M1/2, M1/2)。

因此,我们可以得到

$ \frac{{{k_{1,x}}}}{{{M_2}}} = \frac{{{k_{2,x}}}}{{{M_1}}}\;\;\;\;\;\frac{{{k_{1,y}}}}{{{M_2}}} = \frac{{{k_{2,y}}}}{{{M_1}}} $ (39)

由于M1M2的互质关系,要满足上式则k1, x=k2, x=0,k1, y=k2, y =0,即有θ=$\hat \theta_1=\hat \theta_2 $ϕ=$\hat \phi_1=\hat \phi_2 $,即子阵1和子阵2中只存在一个相同的角度估计值,并且这个值就是唯一正确的DOA估计值,从而消除了阵列测向的角度模糊。

2.2.2 虚拟化阵元方法

互质面阵的虚拟化阵元方法与互质线阵虚拟化阵元方法一致,当使用常规互质面阵时,虚拟阵列会出现不连续部分,这样使DOA估计变得困难,所以,这里直接介绍增广的互质面阵方法,其原理与常规互质面阵一致,只是在阵元数上分布有所不同。

类比于互质线阵,文献[22]提出的虚拟化方法是:假设面阵x轴和y轴方向上分布着2M1x×2M1yM2x×M2y的两个均匀面阵,记为子阵1和子阵2。其阵元间距分别为M2xdM2ydM1xdM1ydd=λ/2,M1xM2x互为质数,M1yM2y互为质数,且M1xM2xM1yM2y。根据文献[22]中的互质面阵,假设M1x = M1y = M1, M2x = M2y = M2

根据文献[22],假定空间有K个不相关的窄带信源入射到上述面阵中,并且θk, ϕk代表第k个信源的仰角和方位角,定义uk =sinθksinϕkvk=sinθkcosϕk。那么可以得到子阵1的接收信号为

$ {\mathit{\boldsymbol{X}}_1}\left( t \right) = \sum\limits_{k = 1}^K {{\mathit{\boldsymbol{a}}_x}\left( {{\theta _k},{\phi _k}} \right)\mathit{\boldsymbol{a}}_y^{\rm{T}}\left( {{\theta _k},{\phi _k}} \right){\mathit{\boldsymbol{s}}_k}\left( t \right) + {\mathit{\boldsymbol{Z}}_1}\left( t \right)} $ (40)

式中:方向矢量可以分别表示为:ax(θk, ϕk) =ax(vk) = [1, …, e-jπ(2M1x-1)M2xvk]Tay(θk, ϕk) = ay (uk) = [1, …, e-jπ(2M1y-1)M2yuk]T,在对接收信号矢量化后可以得到信号

$ {\mathit{\boldsymbol{x}}_1}\left( t \right) = {\rm{vec}}\left( {{\mathit{\boldsymbol{X}}_1}\left( t \right)} \right) = \left( {{\mathit{\boldsymbol{A}}_x} \odot {\mathit{\boldsymbol{A}}_y}} \right)\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{z}}_1}\left( t \right) $ (41)

式中:Ax=[ax(v1), ax(v2), …, ax(vk)],Ay=[ax(u1), ax(u2), …, ax(uk)]。同样的对子阵2作相同处理,可以得到子阵2接收信号在矢量化后的信号为

$ {\mathit{\boldsymbol{x}}_2}\left( t \right) = {\rm{vec}}\left( {{\mathit{\boldsymbol{X}}_2}\left( t \right)} \right) = \left( {{\mathit{\boldsymbol{B}}_x} \odot {\mathit{\boldsymbol{B}}_y}} \right)\mathit{\boldsymbol{s}}\left( t \right) + {\mathit{\boldsymbol{z}}_2}\left( t \right) $ (42)

式中:Bx=[bx(v1), bx(v2), …, bx(vk)],By=[bx(u1), bx(u2), …, bx(uk)]分别为子阵2在平面方向的方向矩阵,同样的可以将子阵2方向矩阵中的方向矢量表示为bx (θk, ϕk) = bx(vk) = [1, …, e-jπ(M2x-1)M1x vk]Tby (θk, ϕk) = by (uk) = [1, …, e-jπ(M2y-1)M1y uk]T。对式(41,42)求互协方差矩阵可得[22]

$ \begin{array}{l} {\mathit{\boldsymbol{R}}_{12}} = E\left[ {{\mathit{\boldsymbol{x}}_1}\left( t \right)\mathit{\boldsymbol{x}}_2^{\rm{H}}\left( t \right)} \right] = \\ \;\;\;\;\;\;\;\;\left( {{\mathit{\boldsymbol{A}}_x} \odot {\mathit{\boldsymbol{A}}_y}} \right){\mathit{\boldsymbol{R}}_s}{\left( {{\mathit{\boldsymbol{B}}_x} \odot {\mathit{\boldsymbol{B}}_y}} \right)^{\rm{H}}} + {\mathit{\boldsymbol{z}}_{12}} \end{array} $ (43)
$ \begin{array}{l} {\mathit{\boldsymbol{R}}_{21}} = E\left[ {{\mathit{\boldsymbol{x}}_2}\left( t \right)\mathit{\boldsymbol{x}}_1^{\rm{H}}\left( t \right)} \right] = \\ \;\;\;\;\;\;\;\;\left( {{\mathit{\boldsymbol{B}}_x} \odot {\mathit{\boldsymbol{B}}_y}} \right){\mathit{\boldsymbol{R}}_s}{\left( {{\mathit{\boldsymbol{A}}_x} \odot {\mathit{\boldsymbol{A}}_y}} \right)^{\rm{H}}} + {\mathit{\boldsymbol{z}}_{21}} \end{array} $ (44)

式中:Rs = E[s(t)sH(t)] = diag(σ12, σ22, …, σk2)为信源矩阵,再对式(43,44)矢量化后得

$ {\mathit{\boldsymbol{r}}_{12}} = {\rm{vec}}\left( {{\mathit{\boldsymbol{R}}_{12}}} \right) = {\mathit{\boldsymbol{C}}_{12}}\left( {u,v} \right)\mathit{\boldsymbol{p}} $ (45)
$ {\mathit{\boldsymbol{r}}_{21}} = {\rm{vec}}\left( {{\mathit{\boldsymbol{R}}_{21}}} \right) = {\mathit{\boldsymbol{C}}_{21}}\left( {u,v} \right)\mathit{\boldsymbol{p}} $ (46)

其中

$ \begin{array}{l} {\mathit{\boldsymbol{C}}_{12}}\left( {u,v} \right) = {\left( {{\mathit{\boldsymbol{B}}_x} \odot {\mathit{\boldsymbol{B}}_y}} \right)^ * } \odot \left( {{\mathit{\boldsymbol{A}}_x} \odot {\mathit{\boldsymbol{A}}_y}} \right) = \left[ {\mathit{\boldsymbol{b}}_x^ * \left( {{u_1}} \right) \otimes } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{b}}_y^ * \left( {{v_1}} \right) \otimes {\mathit{\boldsymbol{a}}_x}\left( {{u_1}} \right) \otimes {\mathit{\boldsymbol{a}}_y}\left( {{v_1}} \right), \cdots ,\mathit{\boldsymbol{b}}_x^ * \left( {{u_k}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { \otimes \mathit{\boldsymbol{b}}_y^ * \left( {{v_k}} \right) \otimes {\mathit{\boldsymbol{a}}_x}\left( {{u_k}} \right) \otimes {\mathit{\boldsymbol{a}}_y}\left( {{v_k}} \right)} \right] \end{array} $ (47)
$ \begin{array}{l} {\mathit{\boldsymbol{C}}_{21}}\left( {u,v} \right) = {\left( {{\mathit{\boldsymbol{A}}_x} \odot {\mathit{\boldsymbol{A}}_y}} \right)^ * } \odot \left( {{\mathit{\boldsymbol{B}}_x} \odot {\mathit{\boldsymbol{B}}_y}} \right) = \left[ {\mathit{\boldsymbol{a}}_x^ * \left( {{u_1}} \right) \otimes } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{a}}_y^ * \left( {{v_1}} \right) \otimes {\mathit{\boldsymbol{b}}_x}\left( {{u_1}} \right) \otimes {\mathit{\boldsymbol{b}}_y}\left( {{v_1}} \right), \cdots ,\mathit{\boldsymbol{a}}_x^ * \left( {{u_k}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. { \otimes \mathit{\boldsymbol{a}}_y^ * \left( {{v_k}} \right) \otimes {\mathit{\boldsymbol{b}}_x}\left( {{u_k}} \right) \otimes {\mathit{\boldsymbol{b}}_y}\left( {{v_k}} \right)} \right] \end{array} $ (48)

所以得到的接收信号为

$ \mathit{\boldsymbol{r}} = {\left[ {\mathit{\boldsymbol{r}}_{12}^{\rm{T}},\mathit{\boldsymbol{r}}_{21}^{\rm{T}}} \right]^{\rm{T}}} = \mathit{\boldsymbol{C}}\left( {u,v} \right)\mathit{\boldsymbol{p}} $ (49)

由式(19)中的过差分集合形式,文献[22]推导了和阵列集合形式,与式(19)类似,可以得到和阵列形式如下

$ \begin{array}{l} {{S'}_{{M_1}{M_2}}} = \left\{ {\left( {{M_1}{m_2} + {M_2}{m_1}} \right),0 \le {m_2} \le {M_2} - 1,} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {0 \le {m_1} \le 2{M_1} - 1} \right\} \end{array} $ (50)

从式(47,48)中可以得到的是,利用差分集合形式或者和阵列集合形式都无法表示其相位上的距离信息,故文献[22]用和-差分集合来表示互质面阵的矢量化后的位置信息

$ \begin{array}{l} S = \left\{ {\left( {{M_{1x}}{m_{2x}} + {M_{1y}}{m_{2y}}} \right) - } \right.\\ \;\;\;\;\;\;\left. {\left( {{M_{2x}}{m_{1x}} + {M_{2y}}{m_{1y}}} \right)} \right\} \end{array} $ (51)

式中:0≤m2xM2x-1, 0≤m2yM2y-1, 0≤m1x≤2M1x -1, 0≤m1y≤2M1y-1。从得到的式(51)和互质线阵的差分集合相比可知,可以将其看作是两个差分集合相加的形式,类似于互质线阵,文献[22]给出的互质面阵中连续均匀的紧凑面阵部分从-(M1xM2xuk + M1yM2yvk)d到(M1xM2xuk + M1yM2yvk)d。所以相对于原有的面阵,互质面阵的空间自由度有所增加。与互质线阵相同的是利用得到的连续的部分的均匀面阵就能使用传统DOA估计算法进行测向,不同的是,互质面阵因为二维角度的限制,所以文献[22]中推导出互质面阵中在有限阵元下可利用的最大空间自由度为

$ {\rm{DO}}{{\rm{F}}_M} = \left[ {\frac{{{{\left( {{S_{UM}} + 5} \right)}^2}}}{8} - 1} \right] $ (52)

式中:SUM=4M1x M1y+M2xM2y-1,[·]为向下取整函数。

得到式(49)后同互质线阵一样,利用空间平滑算法对连续阵元部分进行空间平滑,文献[24]中指出的二维空间平滑算法是利用小块面阵在整个面阵部分进行平滑,此处可选取0d到(M1xM2xuk + M1yM2yvk)d处部分大小空间平滑,假定x(t)为此部分接收信号,Rm1m2 =ARssAH + σn2I为协方差矩阵。那么其它子块协方差可表示为

$ \begin{array}{l} {{R'}_{{m_1}{m_2}}} = \mathit{\boldsymbol{AD}}_x^{\left( {{m_1} - 1} \right)}\mathit{\boldsymbol{D}}_y^{\left( {{m_2} - 1} \right)}{\mathit{\boldsymbol{R}}_{ss}}{\left( {\mathit{\boldsymbol{D}}_y^{\left( {{m_2} - 1} \right)}} \right)^{\rm{H}}}\\ \;\;\;\;\;\;\;\;\;\;\;\;{\left( {\mathit{\boldsymbol{D}}_x^{\left( {{m_1} - 1} \right)}} \right)^{\rm{H}}}{\mathit{\boldsymbol{A}}^{\rm{H}}} + \sigma _n^2\mathit{\boldsymbol{I}} \end{array} $ (53)

式中:Dx=diag(v1, v2, …, vk),Dy=diag(u1, u2, …, uk),1≤m1, m2M1M2。与一维空间平滑算法类似的是,对每个子阵列得到的协方差矩阵相加取平均可得[24]

$ \mathit{\boldsymbol{\dddot R}} = \frac{1}{{M_1^2M_2^2}}\sum\limits_{{m_1} = 1}^{{M_1}{M_2}} {\sum\limits_{{m_2} = 1}^{{M_1}{M_2}} {{{\mathit{\boldsymbol{R'}}}_{{m_1}{m_2}}}} } $ (54)

利用式(54)中得到的协方差矩阵,就能消除信源间的相干性,可知上述协方差矩阵也满足与类似于式(22)的协方差分解式,从而可以利用子空间类算法实现DOA估计。

以子空间类算法的ESPRIT算法为例,通过对$\boldsymbol{\dddot R} $特征分解得到信号子空间Espv,得到Espv的前M1M2(M1M2-1)行和后M1M2(M1M2-1)行的构成的矩阵Espv1Espv2。构造矩阵φpv = Espv1+Espv2后,采用同面阵解模糊方法中子阵i的DOA估计方法,得到x轴仰角与方位角信息,再用同样方法用Espv右乘得到的满秩矩阵估计值构造新信号子空间矩阵Espv′,无噪声情况下Espv′ = Espv,利用Espv′的前M1M2(M1M2-1)行和后M1M2(M1M2-1)行的旋转不变性得到y轴仰角与方位角信息,并且最终仰角与方位角自动匹配。此时的虚拟阵列均匀且阵元间距为半波长,所以得到的估计值即为准确的测向角。

下面给出子空间类算法ESPRIT算法下几种方法的DOA估计性能比较。本文中所有的仿真均采用蒙特卡洛仿真来评估算法性能。假设空间远场信号源入射角为(θ1, φ1)=(5°, 15°),(θ2, φ2)=(5°, 15°)。图 6的仿真条件为普通均匀面阵阵元数M×N=4×6, 互质面阵解模糊方法阵元数M1=4, M2=3, 虚拟化阵元方法互质线阵阵元数M1=2, M2=3。保证总阵元数相同,快拍数L=300条件下的几种方法随着信噪比条件下的DOA估计性能图。

图 6 互质面阵中不同方法下的DOA估计性能 Figure 6 DOA estimation performance under different methods in coprime plan ar array

可以从图 6中看出,在一般情况下,解模糊方法的DOA估计性能优于虚拟化阵元方法,优于均匀面阵方法,但是随着信噪比的增高,虚拟化阵元方法的DOA估计性能变差。

图 7给出了在图 6相同的阵元数情况下,在信噪比SNR=10 dB条件下,几种方法DOA估计性能随着快拍数的变化情况。

图 7 互质面阵中不同快拍数下算法性能 Figure 7 Algorithm performance under different snapshots in coprime plan ar array

综上所述,互质面阵不仅在阵列天线孔径与空间自由度上相对于均匀面阵有着无可比拟的优势,而且在算法性能上,互质面阵一定程度要优于均匀面阵。

互质面阵利用虚拟化阵元的方法在空间自由度上相比于解模糊方法有较大的提高,但是矢量化带来的新的接收信号变成单快拍信源信号,所以虚拟化阵元方法对于接收信号的快拍数具有较高要求,而且算法复杂度相对于解模糊方法要高。

3 结束语

本文介绍了当今互质阵列空间谱估计的研究进展,主要讲述了互质线阵、互质面阵拓扑结构与数据模型,并介绍了目前互质阵列的两种主流DOA估计方法,包括解模糊方法和虚拟化阵元方法,给出了它们的推导过程。在几种阵列形式下,解模糊方法下的互质阵的子空间类算法下的算法性能要优于普通均匀阵列。但是解模糊方法下的互质阵列DOA估计算法牺牲了空间自由度,而虚拟化阵元方法下的互质阵列DOA估计算法可以提高空间自由度,但是在相同条件下,其算法复杂度要高于解模糊方法,并且对于接收信号的快拍数具有较高要求。当然,不论是解模糊方法还是虚拟化阵元的方法下的互质阵DOA估计,阵列天线的孔径得到了有效扩展,阵列的测向精度与测向能力均得到提高。解模糊方法的鲁棒性与虚拟阵元方法的低复杂度算法将是我们进一步的研究方向。对互质阵列新型结构的探索[25-26]以及将互质阵列用于MIMO雷达的高精度测向方法研究[27-30]在今后一段时间内也必是信号处理领域的研究热点。

参考文献
[1] 张小飞, 汪飞, 陈伟华. 阵列信号处理的理论和应用[M]. 北京: 国防工业出版社, 2013.
ZHANG Xiaofei, WANG Fei, CHEN Weihua. Theory and application of ar ray signal processing[M]. Beijing: National Defense Industry Press, 2013.
[2] 单志勇. 空间谱估计技术及其若干应用[D]. 上海: 上海交通大学, 2005.
SHAN Zhiyong. Spatial spectrum estimation technique and its application[D]. Shanghai: Shanghai Jiao Tong University, 2005.http: //d. wanfangdata. com. cn/Thesis/Y856779
[3] KRIM H, VIBERG M. Two decades of array signal processing research: The parametric approach[J]. IEEE Signal Processing Magazine, 1996, 13(4): 67–94. DOI:10.1109/79.526899
[4] ZHANG Xiaofei, CAO Renzheng. Direction of arrival estimation: Introduction[M]. New Jersey: John Wiley & Sons, Inc, 2017.
[5] STOICA P, ARYE N. MUSIC, maximum likelihood, and Cramer-Rao bound[J]. IEEE Transactions on Acoustics Speech & Signal Processing, 1989, 37(5): 720–741.
[6] ROY R, KAILATH T. ESPRIT-estimation of signal parameters via rotational invariance techniques[J]. IEEE Trans ASSP, 1986, 37(7): 984–995.
[7] 黄蕾. 基于阵列孔径扩展的稳健性DOA算法[J]. 信息通信, 2013(8): 40–41.
HUANG Lei. Robust DOA algorithm based on array aperture extension[J]. Information and Communication, 2013(8): 40–41.
[8] DOGAN M C, MENDEL J M. Applications of cumulates to array processing Part Ⅰ: Aperture extension and array calibration[J]. IEEE Transactions on Signal Processing, 1995, 43(5): 1200–1216. DOI:10.1109/78.382404
[9] PORAT B, FRIEDLANDER B. Direction finding algorithm based on high-order statistics[J]. IEEE Transactions on Signal Processing, 1991, 39(9): 2016–2024. DOI:10.1109/78.134434
[10] YIN H, LIU H. Performance of space-division multiple-access (SDMA)with scheduling[J]. IEEE Transactions on Wireless Communications, 2002, 1(4): 611–618. DOI:10.1109/TWC.2002.804188
[11] HUANG K, ANDREWS J G, HEATH R W. Performance of orthogonal beamforming for SDMA with limited feedback[J]. IEEE Transactions on Vehicular Technology, 2009, 58(1): 152–164. DOI:10.1109/TVT.2008.925003
[12] VAIDYANATHAN P P, PAL P. Sparse sensing with co-prime samplers and arrays[J]. IEEE Transactions on Signal Processing, 2011, 59(2): 573–586. DOI:10.1109/TSP.2010.2089682
[13] ZHOU C, SHI Z, GU Y, et al. DECOM: DOA estimation with combined MUSIC for c oprime array[C]// International Conference on Wireless Communications & Signal Processing.[S.l.]:IEEE, 2013:1-5.http://ieeexplore.ieee.org/document/6677080/
[14] PAL P, VAIDYANATHAN P P. Coprime sampling and the music algor ithm[C]//Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop. [S.l.]:IEEE, 2011:289-294.http://ieeexplore.ieee.org/document/5739227/
[15] SHEN Q, LIU W, CUI W, et al. Low-complexity compressive sens ing based DOA estimation for co-prime arrays[C]//International Conference on Digital Signal Processing.[S.l.]: IEEE, 2014:754-758.http://ieeexplore.ieee.org/document/6900765/
[16] JIA T, WANG H, SHEN X, et al. Direction of arrival estimation with co-prime arrays via compressed sensing methods[C]// Oceans.[S.l.]: IEEE, 2016:1-5.http://ieeexplore.ieee.org/document/7485484/
[17] LIU A, YANG Q, ZHANG X, et al. Direction-of-arrival estimation for coprime array using compressive sensing based array interpolation[J]. Inter national Journal of Antennas and propagation, 2017, 2017(3): 1–10.
[18] HU Bin, LV Weihua, ZHANG Xiaofei. 2D-DOA estimation for coprime L-shaped arrays with propagator method[C]//4th National Conference on Electrical Electronics and Computer Engineering (NCEECE 2015). Xi'an:[s.n.], 2015:1551-1556.https://www.atlantis-press.com/proceedings/nceece-15/25847147
[19] LIU S, YANG L S, WU D C, et al. Two-dimensional DOA estimation using a co-prime symmetric cross array[J]. Progress in Electromagnetics Research C, 2014, 54: 67–74. DOI:10.2528/PIERC14081005
[20] WU Q, SUN F, LAN P, et al. Two-dimensional direction-of-arrival estimation for co-prime planar arrays: A partial spectral search approach[J]. IEEE Sensors Journal, 2016, 16(14): 5660–5670. DOI:10.1109/JSEN.2016.2567422
[21] ZHENG W, ZHANG X, SHI Z. Two-dimensional direction of arrival estimation for coprime planar arrays via a computationally efficient one-dimension partial spectral search[J]. IET Radar Sonar and Navigaion, 2017.
[22] SHI Junpeng, HU Guoping, ZHANG Xiaofei, et al. Sparsity-based 2-D DOA estimation for co-prime array: From sum-difference co-array viewpoint[J]. IEEE Transactions on Signal Processing, 2017, 65(21): 5591–5604. DOI:10.1109/TSP.2017.2739105
[23] VAIDYANATHAN P P, PAL P. Sparse sensing with coprime arrays[J]. Signals, Systems and Computers, 2010, 45(2): 1405–1409.
[24] CHEN Y M. On spatial smoothing for two-dimensional direction -of-arrival estimation of coherent signals[J]. IEEE Transactions on Signal Processing, 1997, 45(7): 1689–1696. DOI:10.1109/78.599939
[25] QIN S, ZHANG Y D, AMIN M G. Generalized coprime array configurations for direction-of-arrival estimation[J]. IEEE Transactions on Signal Processing, 2015, 63(6): 1377–1390. DOI:10.1109/TSP.2015.2393838
[26] ZHENG W, ZHANG X, ZHAI H. A generalized coprime planar array geometry for two-dimensional DOA estimation[J]. IEEE Communications Letters, 2017, 21(5): 1072–1078.
[27] WANG X, YANG S, WANG X, et al. Concurrent exploration of MIMO radar and co-prime array for faster and higher resolution scanning[C]//2014 IEEE 48th Asilomar Conference on Signals, Systems and Computers.[S.l.]: IEEE, 2014:153-156.http://ieeexplore.ieee.org/document/7094417/
[28] LI J, SHEN M, DING J. Direction of arrival estimation for co-prime MIMO radar based on unitary root-MUSIC[C]//2nd International Conference on Wireless Communication and Sensor Network.[S.l.]: World Scientific, 2015:307-315.http://www.worldscientific.com/doi/10.1142/9789813140011_0037
[29] LI J, JIANG D, ZHANG X. DOA estimation based on combined unitary ESPRIT for coprime MIMO radar[J]. IEEE Communications Letters, 2017, 21(1): 96–99. DOI:10.1109/LCOMM.2016.2618789
[30] JIA Y, ZHONG X, GUO Y, et al. DOA and DOD estimation based on bistatic MIMO radar with co-prime array[C]// Radar Conference.[S.l.]:IEEE, 2017:0394-0397.http://ieeexplore.ieee.org/document/7944234/