CN111239679B - 一种用于互质面阵下相干信源doa估计的方法 - Google Patents

Abstract

本发明公开了一种互质面阵下相干信源的二维波达方向(Direction of arrival,DOA)估计方法。该方法首先引入空间平滑技术对接收数据的协方差矩阵进行预处理，使其恢复列满秩。在此基础上，该方法基于阵列的旋转不变性对相干信源的波达方向进行估计，同时借助阵列中两个子面阵的互质关系消除模糊值。本发明能实现空间相干信源的有效定位，且DOA估计性能优异。此外本发明无需额外的角度配对方法就能实现信号二维参数的自动配对，在参数估计中也无需进行二维谱峰搜索，因此算法复杂低。

Description

一种用于互质面阵下相干信源DOA估计的方法

技术领域

本发明属于阵列信号处理技术领域，是一种相干信号源的定位方法，特别涉及一种用于互质面阵下相干信源DOA估计的方法。

背景技术

信号空间谱估计作为阵列信号处理领域的一个重要分支，在雷达、通信、声纳、勘测、射电天文及生物医学工程等众多军事及国民经济领域都有广泛的应用。经典的超分辨率DOA估计算法有子空间分解类算法、子空间拟合类算法等，这些算法在理想情况下都具有良好的空间谱估计性能。在实际环境中，由于多径传播等因素的影响，存在大量的相干信号源。由于相干信号源信号子空间与噪声子空间相互渗透，因此传统的DOA估计算法，如MUSIC，ESPRIT等信号子空间类算法，无法对相干信号源进行有效分辨或测向。因此，相干信号源的估计问题始终是一个十分棘手的问题。目前对相干信号源的处理方法一般分为空间平滑类算法、矩阵重构类算法和非降维处理算法。另一方面，近年来，互质阵因其相比于传统均匀阵列拥有更大的阵列孔径和更高的测向精度而受到越来越多专家学者的关注，俨然已成为阵列信号领域的一个研究热点，而目前对互质阵列下相干信源的谱估计算法的研究还存在明显的不足。

发明内容

为了解决现有技术中相干信源在信号参数估计中接收信号协方差矩阵的秩亏损的问题，本发明提供一种用于互质面阵下相干信源DOA估计的方法，该方法实现互质面阵下相干信源的二维DOA估计。

为实现上述目的，本发明采用的技术方案为：

一种用于互质面阵下相干信源DOA估计的方法，首先，利用空间平滑技术对接收数据的协方差矩阵进行预处理，恢复接收信号协方差矩阵的满秩特性，然后基于旋转不变性，运用二维ESPRIT算法对互质面阵下相干信源的波达方向进行估计，并利用阵列中两个子面阵的互质特性消除模糊值，从而实现相干信源的高分辨DOA估计。

进一步的，利用空间平滑技术解相干，包括以下步骤：

互质面阵分成两个子阵，分别记为子阵1和子阵2，子阵阵元只在原点处重合，其中，子阵1是阵元数为M₁×M₁的均匀面阵，子阵2是阵元数为M₂×M₂的均匀面阵；子阵1在X轴和Y轴方向上的阵元间距为d₁＝M₂λ/2，其中：λ为入射信号波长，子阵2在X轴和Y轴方向上的阵元间距为d₂＝M₁λ/2；因此，互质面阵的阵元总数为

根据互质面阵的结构特点，分别对互质面阵的两个子面阵进行空间平滑，将互质面阵的子面阵i(i＝1,2)分为若干个部分重叠的大小为P×Q(1<P≤M_i-1，1<Q≤M_i-1)的子阵，当i＝1时，M_i即M₁，当i＝2时，M_i即M₂，第(p,q)个子阵的阵元位置为L_i,(p,q)＝{((m'+p-1)d_i,(n'+q-1)d_i)|0≤m'≤P-1,0≤n'≤Q-1}，其中：当i＝1时，d_i即为d₁，当i＝2时，d_i即为d₂，其中：m'，n'分别是范围[0,P-1]和[0,Q-1]内的整数，接收信号向量x_i,(p,q)(t)为：

其中，1≤p≤M_i-P+1，1≤q≤M_i-Q+1，i＝1,2；A_i,(1,1)＝A_iy,(1,1)⊙A_ix,(1,1)为子面阵i的第(1,1)个子阵的方向矩阵，A_i,(1,1)在X轴和Y轴方向的方向矢量分别为A_ix,(1,1)和A_iy,(1,1)；

和

分别表示K×K阶对角矩阵Φ_ix的p-1次幂和Φ_iy的q-1次幂，其中，

其中：e为自然常数，j为虚数单位，π为圆周率，θ_k,φ_k分别表示第k(k＝1,2,…,K)条路径的入射俯仰角和方位角；相干信号源矢量为s(t)＝αs₀(t)，α＝[α₁,α₂,…,α_K]^T为由一系列复常数组成的相干系数矩阵，[·]^T表示矩阵的转置，s₀(t)为生成信源；n_i,(p,q)(t)为子面阵i的第(p,q)个子阵的噪声矢量；因此，子面阵i的第(p,q)个子阵接收信号的协方差矩阵R_i,(p,q)表示为

其中，

为x_i,(p,q)(t)的共轭转置；故子面阵i的二维空间平滑的协方差矩阵

由其全部平滑子阵的协方差矩阵的均值

表示：

进一步的，运用二维ESPRIT算法对互质面阵下相干信源的波达方向进行估计，包括以下步骤：

将子面阵i的二维空间平滑协方差矩阵

进行特征分解，得到的K个较大特征值对应的特征矢量构成的信号子空间E_is，将E_is分解成E_ix＝E_is(1:P(Q-1),:)，E_iy＝E_is(P+1:PQ,:)，其中：E_ix为E_is的1到P(Q-1)行，E_iy为E_is的P+1到PQ行；若用D_l(·)表示用矩阵的第l(l＝1,2,…,Q-1)行构造的对角阵，则E_ix、E_iy分别表示为：

E_ix＝A_i1,(1,1)T_i，E_iy＝A_i1,(1,1)Φ_iyT_i (10)

其中，A_i1,(1,1)＝[A_ix,(1,1)D₁(A_iy,(1,1)),A_ix,(1,1)D₂(A_iy,(1,1)),…,A_ix,(1,1)D_Q-1(A_iy,(1,1))]，T_i为K×K阶的满秩矩阵；将式(10)化简得：

E_iy＝E_ixT_i ^-1Φ_iyT_i＝E_ixΨ_i (11)

其中，Ψ_i＝T_i ^-1Φ_iyT_i，T_i ^-1为T_i的逆运算，Ψ_i的特征值为Φ_iy的对角线元素；

由最小二乘法则得到Ψ_i的估计值

其中，

为E_ix的广义逆；

将

特征值分解得到Φ_iy的估计值

并根据

的特征矢量得到T_i ^-1的估计值

在不考虑噪声的情况下有：

其中，Π_i为一个置换矩阵，

表示Π_i的逆运算，由于

和

的特征值一致，因此将

进行特征值分解得到u_ik＝sinθ_k sinφ_k,k＝1,2,…,K的估计值

u_ik是公式sinθ_ksinφ_k,k＝1,2,…,K的值：

其中，

为

的第k个特征值，

表示取

的相角；

对信号子空间进行重构E_is得

其中：

为对矩阵T_i的估计值

求逆；对E'_is进行行变换后得到E″_is，E″_is表示为

将E″_is分解成E'_ix＝E″_is(1:Q(P-1),:)，E'_iy＝E″_is(Q+1:PQ,:)，即E'_ix为E″_is的1到Q(P-1)行，E'_iy为E″_is的Q+1到PQ行；

构造矩阵A_i2,(1,1)＝[A_iy,(1,1)D₁(A_ix,(1,1)),A_iy,(1,1)D₂(A_ix,(1,1)),…,A_iy,(1,1)D_P-1(A_ix,(1,1))]，其中：D_l(·)表示用矩阵的第l(l＝1,2,…,P-1)行构造的对角阵；

则有：

E'_ix＝A_i2,(1,1)Π_i ^-1，E'_iy＝A_i2,(1,1)Φ_ixΠ_i ^-1 (15)

化简得：

(E'_ix)E'_iy＝Π_iΦ_ixΠ_i ^-1 (16)

在不考虑噪声的条件下有：

得到v_ik＝cosθ_k sinφ_k,k＝1,2,…,K的估计值

v_ik是公式cosθ_k sinφ_k,k＝1,2,…,K的值：

其中，

为矩阵(E'_ix)E'_iy的第k个对角元素，

表示取

的相角；由式(13)和(17)，u_ik与v_ik的估计值有相同的列模糊，因此将俯仰角和方位角进行自动配对，利用已配对的

和

得到入射信号的DOA估计。

进一步的，利用阵列中两个子面阵的互质特性消除模糊值，包括以下步骤：

通过对比子阵1和子阵2得到误差最小的

和

再由

和

分别得到真实角度θ_k、φ_k的估计值

与现有技术相比，本发明具有以下有益效果：

(1)、本发明研究了相干信源在互质面阵下的信号二维参数估计问题这一前沿课题，通过借助空间平滑技术解决了相干信源在信号参数估计中接收信号协方差矩阵的秩亏损问题；

(2)、本发明算法在信号参数估计过程中不需要进行二维谱峰搜索，算法复杂度低；

(3)、本发明能实现相干信源二维信号参数(即俯仰角和方位角)的自动配对，不需要额外的角度配对过程

附图说明

图1是互质面阵的阵列结构图。

图2是互质面阵的空间平滑结构示意图。

图3是本发明在SNR＝10dB时的二维空间谱估计结果散布图。

图4是本发明在不同快拍数下的角度估计性能对比图。

图5是本发明在不同阵元数下的角度估计性能对比图。

图6是本发明在不同平滑子阵下的角度估计性能对比图。

图7是本发明在SNR＝15dB时的互质面阵相干信源估计各算法在不同快拍数下性能对比图。

图8是本发明在快拍数J＝400时的互质面阵相干信源估计各算法在不同信噪比下性能对比图。

具体实施方式

下面结合实施例对本发明作更进一步的说明。

一、数据模型

互质面阵的拓扑结构如图1所示。与传统面阵相比，互质面阵可以分成两个子阵，分别记为子阵1和子阵2，子阵阵元只在原点处重合。其中，子阵1为阵元数为M₁×M₁的均匀面阵，子阵2由阵元数为M₂×M₂的均匀面阵构成。子阵1在X轴和Y轴方向上的阵元间距为d₁＝M₂λ/2，其中：λ为入射信号波长，子阵2在X轴和Y轴方向上的阵元间距为d₂＝M₁λ/2。因此，互质面阵的阵元总数为

阵元位置L_s可用如下集合表示：

L_s＝{(zd₁,nd₁)|0≤z,n≤M₁-1}∪{(pd₂,qd₂)|0≤p,q≤M₂-1} (1)

其中：z和n为范围[0,M₁-1]内的整数，p和q为范围[0,M₂-1]内的整数；

假设远场窄带信源s₀(t)产生的K条多径入射到图1所示的互质面阵上，θ_k,φ_k分别表示第k(k＝1,2,…,K)条路径的入射俯仰角和方位角。考虑互质面阵的子阵i(i＝1,2)，其在X轴和Y轴上阵元的方向矢量a_ix(θ_k,φ_k)和a_iy(θ_k,φ_k)可以分别表示为：

其中，e为自然常数，j为虚数单位，π为圆周率，当i＝1时，d_i即为d₁，当i＝2时，d_i即为d₂。因此，子阵i在X轴和Y轴上的方向矩阵可分别表示为A_ix＝[a_ix(θ₁,φ₁),a_ix(θ₂,φ₂),…,a_ix(θ_K,φ_K)]，A_iy＝[a_iy(θ₁,φ₁),a_iy(θ₂,φ₂),…,a_iy(θ_K,φ_K)]。则子阵i在X轴上阵元的接收信号x_i1(t)为：

x_i1(t)＝A_ixs(t)+n_i1(t) (3)

其中，相干信号源矢量为s(t)＝αs₀(t)，α＝[α₁,α₂,…,α_K]^T为由一系列复常数组成的相干系数矩阵，[·]^T表示矩阵的转置，s₀(t)为生成信源。n_i1(t)是子阵i在X轴上阵元的产生的噪声矢量，s(t)为信源矢量，其矩阵维度为K×1。因此，子阵i中平行于X轴的第m(m＝1,2,…,M_i)个分阵列接收信号x_im(t)表示为：

x_im(t)＝A_ixΦ_iy ^(m-1)s(t)+n_im(t) (4)

其中，

diag(·)代表使用向量的元素作为对角元素的对角矩阵。n_im(t)为子面阵i中第m个子阵列的噪声。故整个子阵i的接收信号矩阵x_i(t)为：

将上式改写为：

其中，A_i＝A_iy⊙A_ix为子面阵i的方向矩阵，⊙表示Khatri-Rao积，

表示Kronecker积，n_i(t)为整个子面阵i上的噪声。

定义子阵i的接收信号的协方差矩阵为R_xi＝E[x_i(t)x_i ^H(t)]＝A_iSA_i ^H+σ²I，其中，σ²为噪声方差，E[·]表示求数学期望，x_i ^H(t)为x_i(t)的共轭转置，A_i ^H为A_i的共轭转置，I表示单位阵。S为矩阵维度为K×K的信源协方差矩阵，S表示如下：

S＝E[s(t)s^H(t)]＝E[|s₀(t)|²]αα^H (7)

其中，s^H(t)为s(t)的共轭转置，α^H为α的共轭转置。

由上式可知，协方差矩阵R_xi不再满足满秩特性，导致传统的基于子空间的空间谱估计方法失效。为解决这类难题，本发明将采用空间平滑技术解相干。

二、空间平滑

根据互质面阵的结构特点，分别对互质面阵的两个子面阵进行空间平滑。将互质面阵的子面阵i(i＝1,2)分为若干个图2所示的部分重叠的大小为P×Q(1<P≤M_i-1，1<Q≤M_i-1)的子阵，当i＝1时，M_i即M₁，当i＝2时，M_i即M₂，第(p,q)个子阵的阵元位置为L_i,(p,q)＝{((m'+p-1)d_i,(n'+q-1)d_i)|0≤m'≤P-1,0≤n'≤Q-1}，其中：m'，n'分别是范围[0,P-1]和[0,Q-1]内的整数，接收信号向量x_i,(p,q)(t)为：

其中，1≤p≤M_i-P+1，1≤q≤M_i-Q+1，i＝1,2。A_i,(1,1)＝A_iy,(1,1)⊙A_ix,(1,1)为子面阵i的第(1,1)个子阵的方向矩阵，其在X轴和Y轴方向的方向矢量分别为A_ix,(1,1)和A_iy,(1,1)。

和

分别表示K×K阶对角矩阵Φ_ix的p-1次幂和Φ_iy的q-1次幂，其中，

n_i,(p,q)(t)为子面阵i的第(p,q)个子阵的噪声矢量。因此，子面阵i的第(p,q)个子阵接收信号的协方差矩阵R_i,(p,q)可表示为

其中，

为x_i,(p,q)(t)的共轭转置。故子面阵i的二维空间平滑的协方差矩阵

可由其全部平滑子阵的协方差矩阵的均值

表示：

三、二维ESPRIT算法

将子面阵i的二维空间平滑协方差矩阵

进行特征分解，得到的K个较大特征值对应的特征矢量构成的信号子空间E_is，将E_is分解成E_ix＝E_is(1:P(Q-1),:)，E_iy＝E_is(P+1:PQ,:)，其中：E_ix为E_is的1到P(Q-1)行，E_iy为E_is的P+1到PQ行。若用D_l(·)表示用矩阵的第l(l＝1,2,…,Q-1)行构造的对角阵(字母D单独而言是没有意义的，数学中是diag的缩写)，则E_ix、E_iy可表示为：

E_ix＝A_i1,(1,1)T_i，E_iy＝A_i1,(1,1)Φ_iyT_i (10)

其中，A_i1,(1,1)＝[A_ix,(1,1)D₁(A_iy,(1,1)),A_ix,(1,1)D₂(A_iy,(1,1)),…,A_ix,(1,1)D_Q-1(A_iy,(1,1))]，T_i为K×K阶的满秩矩阵。将上式化简得：

E_iy＝E_ixT_i ^-1Φ_iyT_i＝E_ixΨ_i (11)

其中，Ψ_i＝T_i ^-1Φ_iyT_i，T_i ^-1为T_i的逆运算，Ψ_i的特征值为Φ_iy的对角线元素。

由最小二乘法则可得到Ψ_i的估计值

其中，

为E_ix的广义逆；

将

特征值分解得到Φ_iy的估计值

并根据

的特征矢量得到T_i ^-1的估计值

在不考虑噪声的情况下有：

其中，Π_i为一个置换矩阵，

表示Π_i的逆运算。由于

和

的特征值一致，因此将

进行特征值分解可得u_ik＝sinθ_k sinφ_k,k＝1,2,…,K的估计值

(本发明中仅仅是为了方便描述将sinθ_k sinφ_k,k＝1,2,…,K的值用u_ik代替)：

其中，

为

的第k个特征值，

表示取

的相角。

对信号子空间进行重构E_is得

其中：

为对矩阵T_i的估计值

求逆。对E'_is进行行变换后得到E″_is，E″_is表示为

类似地，将E″_is分解成E'_ix＝E″_is(1:Q(P-1),:)，E'_iy＝E″_is(Q+1:PQ,:)，即E'_ix为E″_is的1到Q(P-1)行，E'_iy为E″_is的Q+1到PQ行。

构造矩阵A_i2,(1,1)＝[A_iy,(1,1)D₁(A_ix,(1,1)),A_iy,(1,1)D₂(A_ix,(1,1)),…,A_iy,(1,1)D_P-1(A_ix,(1,1))]，其中，D_l(·)表示用矩阵的第l(l＝1,2,…,P-1)行构造的对角阵(字母D单独而言是没有意义的，数学中是diag的缩写)；

则有：

化简得：

在不考虑噪声的条件下有：

与u_ik类似，我们可以得到v_ik＝cosθ_k sinφ_k,k＝1,2,…,K的估计值

(本发明中仅仅是为了方便描述将cosθ_k sinφ_k,k＝1,2,…,K的值用v_ik代替)：

其中，

为矩阵(E'_ix)E'_iy的第k个对角元素，

表示取

的相角。由式(13)和(17)，u_ik与v_ik的估计值有相同的列模糊，因此可以将俯仰角和方位角进行自动配对，利用已配对的

和

得到入射信号的DOA估计。

四、解模糊

由于互质面阵的阵元间距大于λ/2，因此存在角度模糊问题。大量研究与文献已证明，互质面阵的角度模糊问题可以通过两个子面阵间的互质关系消除，即两个子面阵相同的角度估计值即为真实信源入射角的DOA估计。根据此理论，通过对比子阵1和子阵2得到误差最小的

和

再由下式分别得到真实角度θ_k、φ_k的估计值

下面利用MATLAB仿真对本发明的算法性能进行分析。其中，采用求根均方误差(Root Mean Square Error,RMSE)来评估算法DOA估计性能，RMSE定义如下：

其中，K表示相干信源的入射路径数，L表示Monte Carlo试验次数，本次仿真中L＝200。

和

分别表示第l次Monte Carlo试验时第k条路径入射信号的仰角θ_k和方位角φ_k的估计值，θ_k,l和φ_k,l分别表示其精确值。

图3为本发明在SNR＝10dB时的二维空间谱估计结果散布图。其中，相干信源分别从角度(5°，10°)、(15°，20°)和(25°，30°)入射到阵元数为M₁＝5、M₂＝6的互质面阵，快拍数J＝200。可以看出，本发明提出的算法可以有效估计出相干信源的DOA信息。

图4是本发明在不同快拍数下的角度估计性能对比图。其中，互质面阵的阵元数M₁＝4，M₂＝5，相干信源的入射方向分别为(15°,10°)、(25°,25°)和(35°,40°)。仿真结果显示，本发明的角度估计性能随着快拍数的增加而变好。

图5是本发明在不同阵元数下的角度估计性能对比图。其中，阵元数M₂＝8，快拍数J＝200，相干信源的入射方向分别为(15°,10°)、(25°,25°)和(35°,40°)。仿真结果显示，俯仰角和方位角估计值的RMSE值随阵元数的增加而降低，算法性能渐优。

图6是本发明在不同平滑子阵下的角度估计性能对比图。其中，阵元数M₁＝7,M₂＝8，快拍数J＝200，相干信源的入射方向分别为(15°,10°)、(25°,25°)和(35°,40°),且规定互质面阵的两个子面阵在仿真中采用同样大小的平滑子阵，P＝Q。仿真结果显示，当P＝5时角度估计性能最佳。

图7-8为互质面阵下相干信源DOA估计的不同算法的性能对比图。其中相干信源分别从角度(15°,10°)、(25°,25°)、(35°,40°)入射到阵元数为M₁＝4,M₂＝5的互质面阵。仿真结果显示，本发明(即SS-ESPRIT)的角度估计性能与SS-PM方法的性能相近，不及SS-MUSIC方法的角度估计性能。但是由于SS-MUSIC方法需要进行二维谱峰搜索，算法复杂度高，不便于硬件实现。

本发明首先引入空间平滑技术对信号协方差矩阵进行预处理，解决了信源相干性导致的信号协方差的秩亏损的问题。在此基础上，本发明基于阵列的旋转不变性实现相干信源的高精度定位。

本发明针对互质阵在进行信号参数过程中会出现模糊值的问题，先对阵列中两个子面阵独立运算，最后利用两个子面阵间的互质关系消除模糊值，得到目标DOA信息的真实估计。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种用于互质面阵下相干信源DOA估计的方法，其特征在于，首先，利用空间平滑技术对接收数据的协方差矩阵进行预处理，恢复接收信号协方差矩阵的满秩特性，然后基于旋转不变性，运用二维ESPRIT算法对互质面阵下相干信源的波达方向进行估计，并利用阵列中两个子面阵的互质特性消除模糊值，从而实现相干信源的高分辨DOA估计；

利用空间平滑技术解相干，包括以下步骤：

互质面阵分成两个子阵，分别记为子阵1和子阵2，子阵阵元只在原点处重合，其中，子阵1是阵元数为M₁×M₁的均匀面阵，子阵2是阵元数为M₂×M₂的均匀面阵；子阵1在X轴和Y轴方向上的阵元间距为d₁＝M₂λ/2，其中：λ为入射信号波长，子阵2在X轴和Y轴方向上的阵元间距为d₂＝M₁λ/2；因此，互质面阵的阵元总数为

根据互质面阵的结构特点，分别对互质面阵的两个子面阵进行空间平滑，将互质面阵的子面阵i(i＝1,2)分为若干个部分重叠的大小为P×Q(1<P≤M_i-1，1<Q≤M_i-1)的子阵，当i＝1时，M_i即M₁，当i＝2时，M_i即M₂，第(p,q)个子阵的阵元位置为L_i,(p,q)＝{((m'+p-1)d_i,(n'+q-1)d_i)|0≤m'≤P-1,0≤n'≤Q-1}，其中：当i＝1时，d_i即为d₁，当i＝2时，d_i即为d₂，其中：m'，n'分别是范围[0,P-1]和[0,Q-1]内的整数，接收信号向量x_i,(p,q)(t)为：

其中，1≤p≤M_i-P+1，1≤q≤M_i-Q+1，i＝1,2；A_i,(1,1)＝A_iy,(1,1)⊙A_ix,(1,1)为子面阵i的第(1,1)个子阵的方向矩阵，A_i,(1,1)在X轴和Y轴方向的方向矢量分别为A_ix,(1,1)和A_iy,(1,1)；

和

分别表示K×K阶对角矩阵Φ_ix的p-1次幂和Φ_iy的q-1次幂，其中，

其中：e为自然常数，j为虚数单位，π为圆周率，θ_k,φ_k分别表示第k(k＝1,2,…,K)条路径的入射俯仰角和方位角；相干信号源矢量为s(t)＝αs₀(t)，α＝[α₁,α₂,…,α_K]^T为由一系列复常数组成的相干系数矩阵，[·]^T表示矩阵的转置，s₀(t)为生成信源；n_i,(p,q)(t)为子面阵i的第(p,q)个子阵的噪声矢量；因此，子面阵i的第(p,q)个子阵接收信号的协方差矩阵R_i,(p,q)表示为

其中，

为x_i,(p,q)(t)的共轭转置；故子面阵i的二维空间平滑的协方差矩阵

由其全部平滑子阵的协方差矩阵的均值

表示：

运用二维ESPRIT算法对互质面阵下相干信源的波达方向进行估计，包括以下步骤：

将子面阵i的二维空间平滑协方差矩阵

进行特征分解，得到的K个特征值对应的特征矢量构成的信号子空间E_is，将E_is分解成E_ix＝E_is(1:P(Q-1),:)，E_iy＝E_is(P+1:PQ,:)，其中：E_ix为E_is的1到P(Q-1)行，E_iy为E_is的P+1到PQ行；若用D_l(·)表示用矩阵的第l(l＝1,2,…,Q-1)行构造的对角阵，则E_ix、E_iy分别表示为：

E_ix＝A_i1,(1,1)T_i，E_iy＝A_i1,(1,1)Φ_iyT_i (10)

其中，

T_i为K×K阶的满秩矩阵；将式(10)化简得：

E_iy＝E_ixT_i ^-1Φ_iyT_i＝E_ixΨ_i (11)

其中，Ψ_i＝T_i ^-1Φ_iyT_i，T_i ^-1为T_i的逆运算，Ψ_i的特征值为Φ_iy的对角线元素；

由最小二乘法则得到Ψ_i的估计值

其中，

为E_ix的广义逆；

将

特征值分解得到Φ_iy的估计值

并根据

的特征矢量得到T_i ^-1的估计值

在不考虑噪声的情况下有：

其中，Π_i为一个置换矩阵，

表示Π_i的逆运算，由于

和

的特征值一致，因此将

进行特征值分解得到u_ik＝sinθ_ksinφ_k,k＝1,2,…,K的估计值

u_ik是公式sinθ_ksinφ_k,k＝1,2,…,K的值：

其中，

为

的第k个特征值，

表示取

的相角；

对信号子空间进行重构E_is得

其中：

为对矩阵T_i的估计值

求逆；对E'_is进行行变换后得到E”_is，E”_is表示为

将E”_is分解成E'_ix＝E”_is(1:Q(P-1),:)，E'_iy＝E”_is(Q+1:PQ,:)，即E'_ix为E”_is的1到Q(P-1)行，E'_iy为E”_is的Q+1到PQ行；

构造矩阵A_i2,(1,1)＝[A_iy,(1,1)D₁(A_ix,(1,1)),A_iy,(1,1)D₂(A_ix,(1,1)),…,A_iy,(1,1)D_P-1(A_ix,(1,1))]，其中：D_l(·)表示用矩阵的第l(l＝1,2,…,P-1)行构造的对角阵；

则有：

化简得：

在不考虑噪声的条件下有：

得到v_ik＝cosθ_ksinφ_k,k＝1,2,…,K的估计值

v_ik是公式cosθ_ksinφ_k,k＝1,2,…,K的值：

其中，

为矩阵(E'_ix)E'_iy的第k个对角元素，

表示取

的相角；由式(13)和(17)，u_ik与v_ik的估计值有相同的列模糊，因此将俯仰角和方位角进行自动配对，利用已配对的

和

得到入射信号的DOA估计；

利用阵列中两个子面阵的互质特性消除模糊值，包括以下步骤：

通过对比子阵1和子阵2得到误差最小的

和

再由

和

分别得到真实角度θ_k、φ_k的估计值

Images