一种基于六角形阵列的相干信号二维DOA估计方法
技术领域
本发明属于天线阵列信号处理领域，它特别涉及一种高精度、去相干、二维信号来波方向(DOA)的子空间估计方法。
背景技术
超分辨波达方向估计(DOA)技术是空间谱估计的一个重要研究点，在高分辨阵列信号处理中得到广泛的应用。其角度估计的高精度和其应用前景引起人们极大的关注。多重信号分类算法(MUSIC,multiple signal classification)和基于旋转不变技术的信号参数估计算法(ESPRIT,estimation of signal parameter via rotational inviance techniques)是两种经典的超分辨子空间类算法。其中，MUSIC算法可以应用于任意形状的阵列，而且能估计多个参数的源信号。但在实际应用中,信噪比及快拍数有限，MUSIC算法的测角精度及分辨力因此受到限制,且可测向信号数小于阵元数。在实际环境中，由于多径效应的影响使经典超分辨DOA估计算法的性能迅速下降，甚至完全失效。相干信号的波达方向估计是空间谱估计亟需解决的一个实际问题。
目前，针对相干信号解相干的方法包括：降维处理算法中的空间平滑类算法、矩阵重构算法和非降维处理算法中的频域平滑法和虚拟阵列变化法等。降维处理算法存在阵列孔径损失的缺点，而非降维处理算法往往针对特定的环境。针对相干信号源的估计问题，已提出了许多解决的途径。基于空间平滑技术改进的MUSIC解相干算法具有很好的性能，但需要将阵列分成多个子阵，不仅减小了阵元数和阵列孔径，还对阵列结构有一定的要求，在实际应用中无法得到很好的去相关性能，针对任意形状阵列，如六角形阵列并未有较好的去相干的算法。
本发明是在六角形阵列结构模型下，去解决相干信号的二维DOA估计问题。六角形平面天线阵是由一个位于中心的天线单元和以它为圆心的若干个同心的六边形环阵（每个含阵元个数为6n，n=1时为最靠近圆心的六边形环阵）组成，它也可看为一个中心存在阵元的改进型圆环阵列。此种阵列使主瓣变的更窄，副瓣更低，对主瓣以外的干扰得到有效的抑制。本发明是基于最简单的六角阵列结构进行研究，由于其阵列结构的特点，去相关的降维处理算法已不再适用。本发明根据六角阵自身的结构特点，利用实际阵列输出信号为复信号的特性，利用共轭矩阵的思想来虚拟扩展阵列，并利用基于矩阵分解的算法中合并扩展矩阵来增加相关矩阵的秩的思想，以实现解相干目的，在没有减少阵元数目和阵列孔径的基础上提出了一种基于扩展重构相关矩阵的去相干的二维DOA估计算法。
发明内容
本发明所要解决的技术问题是提供一种在不减少阵元数目和阵列孔径的前提下，利用扩展重构相关矩阵去相干的二维DOA估计方法，该方法能够针对非均匀线性阵列，能够精确 的对相干信号源的来波方向进行估计。
本发明利用了共轭矩阵的思想来扩展阵列，对阵列接收数据求其共轭矩阵为基础，扩展重构相关矩阵，提供了一种基于六角形阵列的、去相干的、二维的DOA估计。快拍接收的原始数据矩阵及其共轭矩阵，等效于可利用的阵列长度为原来的2倍，充分利用了原始数据矩阵和其新矩阵之间的自相关和互相关信息，提高DOA估计的性能。
一种基于六角形阵列的相干信号二维DOA估计方法，包括如下步骤：
步骤一：六角形天线阵列阵元数为N，输出数据快拍数为K，则各阵元输出数据矩阵X(k)=AS(k)+N(k)为N×K的复矩阵，求其共轭矩阵X*(k)，并求出原阵列输出数据的相关矩阵为RXX=E{X(k)·X(k)H},H表示其共轭转置，A为阵列的导向矢量矩阵，S(k)表示源信号矢量矩阵；N(k)表示阵列输出的噪声均值为零方差为σ2的加性高斯白噪声，且与源信号不相关；
步骤二：共轭数据矩阵对于天线阵列来说，相当于其共轭对称虚拟阵列的数据输出，针对共轭虚拟扩展后的阵列，可分为三个子六角形阵列，第一个子六角阵列即为真实的六角形阵列，第二个子六角阵列为以参考阵元为中心，子阵上阵元其中一半属于真实阵列的部分，另一半为共轭虚拟阵列的部分，第三个子六角形阵列为共轭虚拟阵列；
步骤三：根据扩展后阵列的每个子六角形阵列所含的阵元，得出子六角形阵列的输出数据矩阵X1(k)，X2(k)，X3(k)；
步骤四：分别求出三个子六角形阵列的自相关矩阵和互相关矩阵： R X 1 X 1 = E { X 1 ( k ) · X 1 ( k ) H } , ]]> R X 2 X 2 = E { X 2 ( k ) · X 2 ( k ) H } , ]]> R X 3 X 3 = E { X 3 ( k ) · X 3 ( k ) H } , ]]> R X 2 X 1 = E { X 2 ( k ) · X 1 ( k ) H } , ]]> R X 3 X 1 = E { X 3 ( k ) · X 1 ( k ) H } , ]]> R X 3 X 2 = E { X 3 ( k ) · X 2 ( k ) H } ; ]]>
步骤五：构造出新的相关矩阵 R = R X 1 X 1 R X 2 X 2 R X 3 X 3 R X 2 X 1 R X 3 X 1 R X 3 X 2 ; ]]>
步骤六：对R进行奇异值分解，分解出N个特征值为λ1≥λ2≥…≥λP≥λP+1=…=λN=σ2，通过判断大特征值的个数来估计信号源个数，并根据所对应的特征向量分别得到信号子空间Es和噪声子空间矩阵En；
步骤七：利用二维MUSIC算法构建空间谱函数， ![]()
θ和![]()
分别为源信号的俯仰角和方位角，![]()
是随θ和![]()
变化的阵列对接收信号的导向矢量。使俯仰角θ在(0°,90°)和方位角![]()
在(0°,360°)范围内变化时，找出空间谱![]()
极大值点所对应的角度即为源信号的DOA。
本发明的积极效果：本发明通过对阵列接收数据求共轭来虚拟扩展阵列，得到三个子阵列的六角阵，然后根据三个子阵列的自相关矩阵和互相关矩阵来扩展重构子空间算法的相关矩阵，实现了去相干的目的。对新的相关矩阵进行奇异值分解得到信号子空间和噪声子空间，再利用二维MUSIC谱估计方法对相干的来波信号进行DOA估计。本发明在解除源信号间相干 性的同时进一步提高了二维DOA的估计精度性能。
附图说明
图1是六角形阵列的两种阵列结构：7阵元和19阵元的六角形平面阵列结构图；
图2是最简单的六角形阵列接收信号的来波方法的图；
图3是以六角形阵列最小阵元数的阵列为例，得出的共轭扩展后的阵列，并画分为3个子阵列的结构图；
具体实施方式
下面就结合附图进一步详细说明该发明的工作过程。去相干的二维DOA估计过程可以分为三大部分。
第一部分主要是针对六角形阵列接收数据建模、共轭虚拟扩展阵列、根据阵列结构分割为子阵列形式，主要为步骤一和步骤二。
假设有D(D≤N)个远场窄带信号源入射到如附图2所示的六角形平面阵列上，阵元个数为N，阵元间距为r，入射信号的俯仰角和方位角分别为![]()
阵列输出的噪声为零均值的加性高斯白噪声，方差为σ2，且与源信号不相关。
针对含最少阵元数的六角阵为例子来研究阵列的输出数据模型。
取附图2中平面阵列最左边的阵元作为坐标原点即参考阵元，则该阵列的接收信号为：
X(t)=AS(t)+N(t) (1)
其中A表示阵列流形矩阵；S(t)=[s1(t),s2(t),…sD(t)]T表示信号矢量矩阵；N(t)=[n1(t),n2(t),…nN(t)]T表示阵列的噪声矢量矩阵，T表示矩阵的转置运算。
![]()
其中：![]()
c为光速，c=3×108m/s，d=1,2,…D；
![]()
实际应用中，阵列的输出信号为复信号，求其共轭矩阵，来扩展阵列，如图3所示。对于阵列快拍接收的数据X(k)，则其共轭阵列的接收数据为X*(k)，阵列的输出数据可表示为：
X(k)=[x1(k),…xn(k),…xN(k)]T (3)
阵列的相关矩阵表示为：
RXX=E{X(k)·X(k)H} (4)其中，H表示矩阵的共轭转置运算。
则共轭矩阵的数据X*(k)为：
X*(k)=[x1*(k),…xn*(k),…xN*(k)]T (5)
步骤二即如附图2所示，把虚拟扩展后的阵列划分为3个相交的六角形阵列，阵元间隔距离分别为r。
第二部分是针对扩展后阵列的处理，主要包括步骤三到步骤五，包含两部分内容，首先是针对扩展后阵列的每个子六角形阵列所含的阵元，得出3个子六角形阵列的输出数据矩阵X1(k)，X2(k)，X3(k)，其次是分别求出三个子六角形阵列的自相关矩阵和互相关矩阵，扩展重构出新的相关矩阵，达到去相干的目的。
（1）子阵列输出数据的构造
每个子阵列的输出数据可以表达为：
X 1 ( k ) = [ x 1 ( k ) , · · · x k ( k ) , · · · x 7 ( k ) ] T X 2 ( k ) = [ x 5 ( k ) , x 1 ( k ) , x 4 ( k ) , x 4 * ( k ) , x 1 * ( k ) , x 7 * ( k ) , x 6 ( k ) ] T X 3 ( k ) = [ x 1 * ( k ) , x 5 * ( k ) , x 4 * ( k ) , x 3 * ( k ) , x 2 * ( k ) , x 7 * ( k ) , x 6 ( k ) ] T - - - ( 6 ) ]]>
从扩展后的阵列模型看出，子阵2和子阵3可由真实的六角形阵列沿x轴依次向左平移r的距离得到，根据式(3)和(5)，上式又可表述为：
X1(k)=AS(k)+N(k)
X2(k)=A·Φ0S(k)+N(k) (7)
X3(k)=A·Φ02S(k)+N(k)
式中，Φ0为对角矩阵：![]()
（2）扩展重构相关矩阵的构造
根据（1）中，可以得到各个子阵列的空间自相关矩阵和互相关矩阵：
R X 1 X 1 = E { X 1 ( k ) · X 1 ( k ) H } = AR S A H + R N R X 2 X 2 = E { X 2 ( k ) · X 2 ( k ) H } = AΦ 0 R S Φ 0 H A H + R N R X 3 X 3 = E { X 3 ( k ) · X 3 ( k ) H } = A ( Φ 0 2 ) R S ( Φ 0 2 ) H A H + R N R X 2 X 1 = E { X 2 ( k ) · X 1 ( k ) H } = A Φ 0 R S A H + R N R X 3 X 1 = E { X 3 ( k ) · X 1 ( k ) H } = A ( Φ 0 2 ) R S A H + R N R X 3 X 2 = E { X 3 ( k ) · X 2 ( k ) H } = A ( Φ 0 2 ) R S Φ 0 H A H + R N - - - ( 8 ) ]]>
式中：E{□}表示统计期望；H表示共轭转置运算；RS=E{S(k)SH(k)}表示信号协方差矩阵；
RN=σ2I表示阵列噪声的协方差矩阵，I为单位矩阵。
根据上述自相关矩阵和互相关矩阵，合并得到扩展的相关矩阵为：
R = R X 1 X 1 R X 2 X 1 R X 2 X 2 R X 3 X 1 R X 3 X 2 R X 3 X 3 - - - ( 9 ) ]]>
根据式(8)和(9)又可表述为：
R=A[RS Φ0RS Φ0RSΦ0H (Φ02)RS (Φ02)RSΦ0H (Φ02)RS(Φ02)H]AH (10)
+[RN RN RN RN RN RN]
令
G=[RS Φ0RS Φ0RSΦ0H (Φ02)RS (Φ02)RSΦ0H (Φ02)RS(Φ02)H] ,
Rn=[RN RN RN RN RN RN]，则：
R=AGAH+Rn (11)
定理1：设Q是L×M维无零行向量的矩阵（L≤M），P为L×L维对角矩阵，其对角元素互不相等，若rank(Q)=r<L，则rank[Q PQ]=r+1。
定理2：设H是M×L维无零行向量的矩阵（L≤M），P为L×L维对角矩阵，其对角元素互不相等，若rank(H)=r<L，则rank[H HP]=r+1。
假设达到阵列的D(D≤N)个信号都是相干的，则rank(RS)=1，由于各信号角度都不相同，当不存在任意信源的方位角和另一信源的俯仰角之和为90°（即![]()
）和任意两个俯仰角相等时它们的方位角之和为360°（即θi＝θj，![]()
）和信源的方位角不为90°和270°时，Φ0的对角元素是各不相同的，为满秩矩阵。由上面的理论推导可知：rank(RS)+n≤rank(G)≤D，其中n为G中矩阵的个数。则，根据式(11)可知R可以分解出D维的信号子空间，再利用DOA估计算法，如MUSIC算法，获得D个信号的角度信息，即实现了去相关的目的。该方法根据其扩展重构后新相关矩阵的秩，最大能估计的相干信号个数至少为N-1个，增加了可估计相干源数目。
第三部分主要是在扩展重构相关矩阵的基础上，利用子空间算法进行DOA估计的处理，包含步骤六和步骤七。根据第二部分的理论分析可以知道，对于相干源信号，重构的新相关矩阵的秩增加可解除相干，在此基础上，步骤六利用了奇异值分解对矩阵R进行分解，由于噪声为零均值方差为σ2的加性高斯白噪声，矩阵G为满秩，可知能分解出N个特征值为：
λ1≥λ2≥…≥λP≥λP+1=…=λN=σ2 (12)
根据上式关系，可以看出大的特征值的个数即为信号源个数。其对应的特征矢量分别为e1,e2,…,eP,eP+1,…eN则
R = Σ i N λ i e i e i H = EΣ E H - - - ( 13 ) ]]>
其中E=[e1,e2,…,eP,eP+1,…eN],Σ=diag{λ1,λ2,…,λP,σ2,…,σ2}。定义信号子空间为Es=[e1,e2,…,eP]，噪声子空间为En=[eP+1,…,eN]。由此，根据特征值对应的特征向量得到了信号子空间Es和噪声子空间矩阵En；步骤七选择了子空间类算法中的其中一种MUSIC算法，利用噪声子空间矩阵En与导向矢量![]()
的正交性，求其二维空间谱函数，如式(14)所示：
![]()
其中，θ和![]()
分别为源信号的俯仰角和方位角，![]()
是随θ和![]()
变化的阵列对接收信号的导向矢量。使θ和![]()
变化，找出空间谱![]()
的极大值点所对应的角度即为源信号的DOA。