阵列信号处理的基本原理是利用收集到的有用的空间信号，对其进行信号处理，同时对噪声信号进行抑制，从而得到源信号的特征及其所包含的有用的信息^[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]中提出了互质线阵的概念。互质线阵由阵元数分别为M和N的两个均匀线阵组成，其中，M和N为互质数，子阵的阵元间距分别为Nλ/2和Mλ/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(·)表示矢量，X_{i, j}表示矩阵的第(i, j)个元素，I_p表示为P×P的单位矩阵，⊗表示Kronecker积，⊙表示Khatri-Rao积，diag(x)表示以x中元素构成的对角矩阵。

1 互质线阵空间谱估计

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

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

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

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

$ \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个不相关窄带信号入射到上述互质线阵上，其入射角度分别为[θ₁, θ₂, …, θ_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)

式中：a₂₁^T (θ_k)为a₂^T(θ_k)去除第一行元素后构成的新方向矢量。

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

以子空间类算法中的ESPRIT算法为例，解模糊方法通过联合互质线阵中的两个均匀子阵DOA估计结果获得测向角，假定子阵i (i=1, 2)得到接收信号x_i(t)，协方差矩阵R_xi = E[x_i(t)·x_i^H(t)]，对协方差矩阵特征分解后得到信号子空间E_{s, i}，再利用信号子空间的前M_i－1行和后M_i－1行得到两个新矩阵E_{s1, i}和E_{s2, 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)}}} $)，d_i为子阵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)可知，当子阵中阵元间距为d_i=M_iλ/2时，因为一维测向角度范围是(－90 °, 90°)，根据|sinθ－sinθ^a|＜2，可以从式(9)中推导出k的取值范围是k∈[－(M_i－1), …, (M_i－1)]。又因为正弦函数的范围，所以k实际只有M_i个取值。当M_i等于1时，k只能取值为0，此时估计结果只有一个，也就是传统阵列形式下的DOA估计结果。当M_i的值大于1时，互质线阵的子阵会产生M_i个模糊值，其中包含对应的真实角度估计值。

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

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

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

(10)

式中：k₁=－(M₂－1), …, －1, 1, …, (M₂－1)。

同理，阵元数为M₂的子阵2阵元间距为d₂=M₁λ/2有关系式

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

(11)

式中：k₂=－(M₁－1), …, －1, 1, …, (M₁－1)

根据式(10，11)可得

$ \frac{{2{k_1}}}{{{M_2}}} = \frac{{2{k_2}}}{{{M_1}}} $

(12)

由于M₁, M₂互为质数，所以不存在k₁, k₂使式(12)成立。

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

1.2.2 虚拟化阵元方法

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

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

$ \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)

式中：x₁(t)和x₂(t)分别为子阵一和子阵二的接收信号矩阵，互质线阵方向矢量表示为a(θ_k) = [a_M1^T(θ_k), a_M2^T(θ_k)]^T，接收信号的协方差矩阵为R_xx=E[x(t)x^H(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 = [σ₁², σ₂², …, σ_K²]^T，σ_k²表示为第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)

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

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

文献[14]利用空间平滑技术与MUSIC算法进行DOA估计。假设z₁是z消除重复行并且按距离位置排列好后的数据矩阵，z_1i是包含有M₁M₂+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)中得出R_ss= $\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(σ₁², σ₂², …, σ_K²)，A₁₁表示空间平滑算法中第一个参考子阵列的方向矩阵。上述表达式符合基本子空间类DOA估计算法的协方差分解形式，故能利用它得到正确的DOA估计。

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

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

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

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

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

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

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

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

考虑上述的互质面阵，子阵阵元只在原点位置重合，故互质面阵的阵元总数为T=M₁²+M₂²－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)中可以得到子阵i在x轴和y轴上的方向矩阵分别表示为A_ix=[a_ix(θ₁, ϕ₁), a_ix(θ₂, ϕ₂), …, a_ix(θ_k, ϕ_k)]，A_iy=[a_iy(θ₁, ϕ₁), a_iy(θ₂, ϕ₂), …, a_iy(θ_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)

式中：n_i1(t)为子阵i在x轴上阵元的噪声。s(t)∈C^K×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}}} $)，n_im(t)为第m个分阵列的噪声。则可得整个子阵i的接收信号x_i(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)

式中：A_i=A_iy⊙A_ix，n_i(t)为整个子面阵i上的噪声。

2.2 DOA估计方法

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

2.2.1 解模糊方法

以子空间类算法中的ESPRIT算法为例，互质面阵中解模糊方法也是通过联合两个子阵估计结果获得DOA估计。通过子阵i得到的式(28)中的接收信号构造协方差矩阵R_xi = E[x_i (t)x_i^H (t)]，特征分解后得到信号子空间E_sp=(A_iy⊙A_ix)T，T为满秩矩阵，通过信号子空间E_sp的前M_i(M_i－1)行和后M_i(M_i－1)行得到E_sp1和E_sp2，构造类似于式(7)中的方程Φ_p = E_sp1⁺E_sp2，特征分解后通过特征值就得到y轴方向上的仰角与方位角的估计值。

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

文献[20]中详细推导了互质面阵的解模糊方法的DOA估计理论。与互质线阵一致，模糊产生是由于阵元间距大于半波长引起的。假设有空间角度为(θ_k, ϕ_k)(k=1, 2, …, K)的K个非相干窄带信源入射到互质面阵中，θ_k, ϕ_k分别代表第k个信源的仰角和方位角。考虑面阵中的子阵列i，其阵元间距假设为d_i=M_iλ/2，那么子阵i在x轴方向和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)

式中：Δφ_x为x轴方向上相邻阵元的相位差，Δφ_y为y轴方向上相邻阵元间的相位差。k_x和k_y为整数，又因为θ_k∈[0, π]，ϕ_k∈[0, π/2]，所以有－1≤sinϕ_kcosθ_k≤1，0≤sinϕ_ksinθ_k≤1。则k_x和k_y的取值范围分别为

$ {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)

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

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

根据式(29，30)，当阵元间距为d_i=M_jλ/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;i≠j，k_{i, x}和k_{i, y}都是整数，且取值范围为(－M_j, M_j)和(－M_j/2, M_j/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)

式中：k_{1, x}和k_{1, y}为整数，且其取值范围分别为(－M₂, M₂)和(－M₂/2, M₂/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)

式中：k_{2, x}和k_{2, y}为整数，且其取值范围分别为(－M₁, M₁)和(－M₁/2, M₁/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)

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

2.2.2 虚拟化阵元方法

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

类比于互质线阵，文献[22]提出的虚拟化方法是：假设面阵x轴和y轴方向上分布着2M_1x×2M_1y和M_2x×M_2y的两个均匀面阵，记为子阵1和子阵2。其阵元间距分别为M_2xd，M_2yd和M_1xd，M_1yd，d=λ/2，M_1x和M_2x互为质数，M_1y和M_2y互为质数，且M_1x＜M_2xM_1y＜M_2y。根据文献[22]中的互质面阵，假设M_1x = M_1y = M₁, M_2x = M_2y = M₂。

根据文献[22]，假定空间有K个不相关的窄带信源入射到上述面阵中，并且θ_k, ϕ_k代表第k个信源的仰角和方位角，定义u_k =sinθ_ksinϕ_k，v_k=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)

式中：方向矢量可以分别表示为：a_x(θ_k, ϕ_k) =a_x(v_k) = [1, …, e^{-jπ(2M_1x-1)M_2xv_k}]^T，a_y(θ_k, ϕ_k) = a_y (u_k) = [1, …, e^{-jπ(2M_1y-1)M2yu_k}]^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)

式中：A_x=[a_x(v₁), a_x(v₂), …, a_x(v_k)]，A_y=[a_x(u₁), a_x(u₂), …, a_x(u_k)]。同样的对子阵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)

式中：B_x=[b_x(v₁), b_x(v₂), …, b_x(v_k)]，B_y=[b_x(u₁), b_x(u₂), …, b_x(u_k)]分别为子阵2在平面方向的方向矩阵，同样的可以将子阵2方向矩阵中的方向矢量表示为b_x (θ_k, ϕ_k) = b_x(v_k) = [1, …, e^{-jπ(M_2x-1)M_1x v_k}]^T，b_y (θ_k, ϕ_k) = b_y (u_k) = [1, …, e^{-jπ(M_2y-1)M_1y u_k}]^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)

式中：R_s = E[s(t)s^H(t)] = diag(σ₁², σ₂², …, σ_k²)为信源矩阵，再对式(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≤m_2x≤M_2x-1, 0≤m_2y≤M_2y-1, 0≤m_1x≤2M_1x -1, 0≤m_1y≤2M_1y-1。从得到的式(51)和互质线阵的差分集合相比可知，可以将其看作是两个差分集合相加的形式，类似于互质线阵，文献[22]给出的互质面阵中连续均匀的紧凑面阵部分从-(M_1xM_2xu_k + M_1yM_2yv_k)d到(M_1xM_2xu_k + M_1yM_2yv_k)d。所以相对于原有的面阵，互质面阵的空间自由度有所增加。与互质线阵相同的是利用得到的连续的部分的均匀面阵就能使用传统DOA估计算法进行测向，不同的是，互质面阵因为二维角度的限制，所以文献[22]中推导出互质面阵中在有限阵元下可利用的最大空间自由度为

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

(52)

式中：S_UM=4M_1x M_1y+M_2xM_2y－1，[·]为向下取整函数。

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

$ \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)

式中：D_x=diag(v₁, v₂, …, v_k)，D_y=diag(u₁, u₂, …, u_k)，1≤m₁, m₂≤M₁M₂。与一维空间平滑算法类似的是，对每个子阵列得到的协方差矩阵相加取平均可得^[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} $特征分解得到信号子空间E_spv，得到E_spv的前M₁M₂(M₁M₂－1)行和后M₁M₂(M₁M₂－1)行的构成的矩阵E_spv1和E_spv2。构造矩阵φ_pv = E_spv1⁺E_spv2后，采用同面阵解模糊方法中子阵i的DOA估计方法，得到x轴仰角与方位角信息，再用同样方法用E_spv右乘得到的满秩矩阵估计值构造新信号子空间矩阵E_spv′，无噪声情况下E_spv′ = E_spv，利用E_spv′的前M₁M₂(M₁M₂－1)行和后M₁M₂(M₁M₂－1)行的旋转不变性得到y轴仰角与方位角信息，并且最终仰角与方位角自动匹配。此时的虚拟阵列均匀且阵元间距为半波长，所以得到的估计值即为准确的测向角。

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

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

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

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

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

3 结束语

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