到达方向跟踪的方法和装置及其在基站的波束成形的应用

摘要

到达方向跟踪的方法和装置及其在基站的波束成形的应用。提出了一种方向跟踪方法，该方法即使在到达方向的轨迹交叉的情况下，也可以快速地跟踪无线电波的到达方向。通过该跟踪方法，使用前一个方向更新时刻处的状态向量，利用观测器来预测当前方向更新时刻的到达方向的状态向量以计算该到达方向的临时估计值，该状态向量由到达方向以及该到达方向的速度和加速度构成，此外，根据该临时估计值所预测的状态向量，以及该到达方向的预测值来计算当前方向更新时刻的状态向量，以通过该状态向量计算该到达方向的估计值。

说明书

技术领域

本发明涉及无线电波到达方向跟踪方法及其装置。具体地，本发明涉及一种到达方向跟踪方法和装置，用于估计由自适应阵列天线接收到的无线电波的到达方向。

背景技术

近年来，对于使用自适应阵列天线的移动通信系统的研究和开发受到关注。典型的阵列天线包括多个天线振子，这些天线振子设置在不同的空间位置，以使得这些天线振子的轮廓具有特定的几何形状。对到达阵列天线的无线电波(下文中，从信号处理的角度，将无线电波称为信号)的到达方向进行估计的技术是与自适应阵列天线相关的最重要的基本技术之一。

由于其估计精度和计算负荷，基于子空间的方法成为公知的用来估计信号到达方向的方法，在该方法中，采用了信号子空间和噪声子空间之间的正交性。多重信号分类(MUSIC)技术是一种典型的基于子空间的方法(例如，参见R.O.Schmidt，“Multiple emitter location andsignal parameter estimation，”IEEE Trans.Antennas andPropagation，vol.34，no.3，pp.267-280(1986))。利用空间平滑的基于子空间的方法也是公知的用于估计具有完全相关性的相干信号的到达方向的方法。利用空间平滑的基于子空间的方法的典型示例是基于空间平滑的MUSIC技术(例如，参见T.J.Shan，M.Wax and T.Kailath，“On spatial smoothing for direction-of-arrival estimation ofcoherent signals，”IEEE Trans.Acoust.，Speech，Signal Processing，vol.33，no.4，pp.806-811(1985)以及S.U.Pillai and B.H.Kwon，“Forward/backward spatial smoothing techniques for coherentsignals identification，”IEEE Trans.Acoust.，Speech，SignalProcessing，vol.37，no.1，pp.8-15(1989))。

在用于估计非相关信号的到达方向的基于子空间的方法中，首先根据噪声接收阵列数据来估计阵列协方差矩阵，然后通过该阵列协方差矩阵的特征值分解(EVD)来获得信号子空间和噪声子空间。此后，采用该信号子空间和噪声子空间之间的正交性，来估计入射信号的到达方向。另一方面，对于相关信号(包括完全相关信号，即，相干信号)的到达方向的估计，为了抑制入信号的相关性，将具有沿直线以相同的间隔设置在不同的空间位置的M个天线振子的天线(下文中，可以将这种天线称为均匀线性阵列(ULA))分为多个重叠的子阵列，然后对这些子阵列的协方差矩阵求平均，以将空间平均的协方差矩阵的信号子空间的维数恢复为入射信号的数量。由此，可以采用信号子空间和噪声子空间之间的正交关系，按照与用于估计非相关信号的到达方向的基于子空间的方法相同的方式，来估计相关信号的到达方向。

以下对S.U.Pillai and B.H.Kwon，“Forward/backward spatialsmoothing techniques for coherent signals identification，”IEEETrans.Acoust.，Speech，Signal Processing，vol.37，no.1，pp.8-15(1989)中提出的用于估计相干信号的到达方向的基于空间平滑的MUSIC方法进行详细说明。

现假定p窄带信号{S_i(k)}沿角度{θ_i}入射到均匀线性阵列(ULA)上。如果T_S表示抽样间隔，则各个振子所接收到的信号可以表示为公式(1)。

y(k)＝[Y₁(k)，y₂(k)，…，y_M(k)]^T＝A(θ(k))s(t)+w(k)

$>>A>>=>\Delta >>[>a>>(>>\theta >1>>>(>k>)>>)>>,>a>>(>>\theta >2>>>(>k>)>>)>>,>.>.>.>,>a>>(>>\theta >p>>>(>k>)>>)>>]>>>$

$>>a>>(>>\theta >i>>>(>k>)>>)>>>=>\Delta >>[>1>,>>e>>j>>\omega >0>>\tau >>(>>\theta >i>>>(>k>)>>)>>>>,>.>.>.>,>>e>>j>>\omega >0>>>(>M>->1>)>>\tau >>(>>\theta >i>>>(>k>)>>)>>>>>]>T>>>>$

s(k)＝[s₁(k)，s₂(k)，…，s_p(k)]^T

w(k)＝[w₁(k)，w₂(k)，…，w_M(k)]^T

$>>>\omega >0>>=>2>\pi >>f>0>>,>\tau >>(>>\theta >i>>>(>k>)>>)>>>=>\Delta >>>(>d>/>c>)>>sin>>\theta >i>>>(>k>)>>>>\dots (1)$

其中f₀、c和d分别表示载波频率和载波的传播速度，以及振子间隔(半波长度)。(·)^T表示转置(transposition)，而a(θ_i(k))和A分别与阵列响应向量和矩阵相对应。w_i(k)表示具有零均值和方差σ²的空间和时间上不相关的复白高斯噪声。

首先，我们考虑信号的到达方向对于时间为恒定的情况，即，θ_i(k)＝θ_i。在这种情况下，阵列协方差矩阵表示为公式(2)。

$>>R>>=>\Delta >>E>\{>y>>(>k>)>>>y>H>>>(>k>)>>\}>=>A>>R>S>>>A>H>>+>>\sigma >2>>>I>M>>->->->>(>2>)>>>>$

其中E{·}和(·)^H分别表示期望运算和复共轭转置，R_S＝E{s(k)s^H(k)}是入射信号的协方差矩阵，而I_M是M×M单位矩阵。此外，如果将所接收到的数据y_i(k)和y_m(k)之间的相关性r_im定义为r_im＝E{y_i(k)y^*_m(k)}，则关系r_im＝r*_mi成立，其中(·)^*表示复共轭。公式(2)中的阵列协方差矩阵R可以表示为公式(3)。

为了估计相干信号的到达方向{θ_k}，基于空间平滑的MUSIC方法将整个阵列分为L个交叠的子阵列，每一个子阵列都包括m(1≤m≤M)个振子。

图14示出了均匀线性阵列(ULA)中的子阵列。

如图14所示，阵列天线100由以相同的间隔d设置的M个天线振子101构成，并被分成L个交叠的子阵列。在这种结构中，m和L分别表示子阵列大小和子阵列数量，其中关系L＝M-m+1成立。根据公式(1)，通过公式(4)给出第1子阵列的接收向量y_l(k)。

y_l(k)＝[Y_l(k)，y_l+1(k)，…，y_l+M-1(k)]^T＝A_mD^l-1s(k)+w_l(k)

A_m＝[a_m(θ₁)，a_m(θ₂)，…，a_m(θ_P)]

$>>>a>m>>>(>>\theta >i>>)>>=>[>1>,>>e>>j>>\omega >0>>\tau >>(>>\theta >i>>)>>>>,>.>.>.>,>>e>>j>>\omega >0>>>(>m>->1>)>>\tau >>(>>\theta >i>>)>>>>>]>T>>>>$

w_l(k)＝[w_l(k)，w_l+1(k)，…，w_l-m+1(k)]^T    …(4)

其中D是包括e^jω0τ(θ1)，e^jω0τ(θ2)，…，e^{jω0(m-1)τ(θp)}作为元素的对角矩阵，并且l＝1，2，…L。

此外，a_m(θ_i)和A_m分别表示子阵列的响应向量和矩阵。因此，通过公式(5)给出该子阵列的协方差矩阵。

$>>>R>l>>=>E>\{>>y>l>>>(>k>)>>>>y>l>>H>>>(>k>)>>\}>=>>A>m>>>D>>l>->1>>>>R>s>>>>(>>D>>l>->1>>>)>>H>sup>>A>m>Hsup>>+>>\sigma >2>>>I>m>>->->->>(>5>)>>>>$

此外，通过对(5)中的L个子阵列的协方差矩阵{R_l}进行空间平均，来获得公式(6)中给出的协方差矩阵。

$>>>R>‾>>=>>1>L>>>\Sigma >>l>=>1>>L>>>R>l>>->->->>(>6>)>>>>$

然后，可以将该进行了空间平均的协方差矩阵R的特征值分解(EVD)写为公式(7)。

$>>>R>‾>>=>>\Sigma >>i>=>1>>m>>>\lambda >i>>>e>i>>>>e>i>>H>>=>E\Lambda >>E>H>>->->->>(>7>)>>>>$

其中e_i和λ_i分别表示特征向量和特征值，E是具有列{e_i}的矩阵，而Λ是具有元素{λ_i}的对角矩阵。此外，信号向量{e₁，e₂，…，e_p}和噪声向量{e_P+1，e_P+2，…，e_m}所跨越的空间被分别称为信号子空间和噪声子空间。可以使用阵列响应向量来表示信号子空间。基于信号子空间和噪声子空间之间的正交关系的到达方向估计方法被称之为基于子空间的方法。

根据公式(7)中的协方差矩阵的特征值分析，在噪声向量{e_P+1，e_P+2，…，e_m}和属于信号子空间内的子阵列的响应向量a_m(θ_i)之间建立了由公式(8)限定的正交关系。

$>sup>>e>k>Hsup>>>a>m>>>(>>\theta >i>>)>>=>0>->->->>(>8>)>>>>$

其中k＝p+1，…，m。根据这种正交关系，可以计算出如公式(9)所示的频谱P_ssmusic(θ)。

$>>>>P>‾>>ssmusic>>>(>\theta >)>>=>>1>sup>>\Sigma >>k>=>p>+>1>>msup>>|sup>>e>k>Hsup>>>a>m>>>(>\theta >)>>>|>2>>>>->->->>(>9>)>>>>$

其中 $>>>a>m>>>(>\theta >)>>=>[>1>,>>e>>j>>\omega >0>>\tau >>(>\theta >)>>>>,>.>.>.>,>>e>>j>>\omega >0>>>(>m>->1>)>>\tau >>(>\theta >)>>>>>]>T>>.>>>在该基于空间平滑的MUSIC方法中，根据由公式(9)给出的频谱的最高p峰值的位置，来估计入信号的到达方向。$

如公式(7)所示，该基于子空间的方法(包括(空间平滑)MUSIC)需要阵列协方差矩阵的EVD，以获得信号子空间或噪声子空间，从而估计到达方向。然而，在某些实际应用中，特别是在阵列元素的数量很大的情况下，EVD或者单值分解(SVD)变得复杂，并且在应该以实时的方式来估计入射信号随时间变化的方向时，非常耗时。因此，计算强度很大的特征分解处理导致这种基于子空间的到达方向估计方法受到特征分解(EVD或SVD)的限制。在实际的移动通信系统的许多情况下，由于来自呼叫方(移动终端)的信号通过直达路径或者通过信号在诸如建筑物的物体上反射而产生的反射路径，到达基站中的阵列天线，所以用于估计多径传播环境下的相干信号的到达方向的技术起着非常重要的作用。然而，在上述方向估计方法中，由于无法从干涉信号中分辨出目标信号，所以必须计算所有信号的到达方向。因此，为了处理大量的入波，必须在阵列天线中使用许多振子。这导致与阵列天线相关的尺寸和成本的增加。此外，如果目标信号的到达角度例如由于呼叫方(单个源)的移动而随时间变化，则无法使用普通的基于子空间的方法对到达阵列的信号的方向进行高速和高精度的估计，或者无法在基站处形成精确的接收和发送波束。这导致了基站处的接收和发送系统的性能下降。

近来，对无需特征分解的自适应到达方向估计和跟踪方法进行了研究，例如无需特征分解的自适应的基于子空间的方法(SWEDE)(例如，参见A.Eriksson，P.Stoica，and T.Sderstrm，“On-line subspacealgorithms for tracking moving sources，”IEEE Trans.SignalProcessing，vol.42，no.9，pp.2319-2330(1994))。然而，这些方法在相干信号、低信噪比(SNR)，或者少量快照(snapshot)的情况下表现出非常差的性能。此外，SWEDE中所涉及的最小平方(LS)技术需要很高的计算复杂度。

本发明人提出了一种基于通信信号的循环性稳态(cyclostationarity)的到达方向估计和跟踪方法(例如，参见J.Xinand A.Sano，“Direct ions-of-arrival tracking of coherentcyclostationary signals in array processing，”IEICE Trans.Fundamentals，vol.E86-A，no.8，pp.2037-2046(2003))。然而，这种使用LS技术的方法需要相当大的阵列数据长度，因为它采用了被称为入射信号的循环稳态的时间特性。

本发明人还提出了被称作无需特征分解的基于子空间的方法(SUMWE)的到达方向估计方法，该方法不需要特征分解，并且计算效率高(例如，参见J.Xin and A.Sano，“Computationally efficientsubspace-based method for direction-of-arrival estimation withouteigendecomposition，”IEEE Trans.Signal Processing，vol.52，no.4，pp.876-893(2004))。然而，该方法没有考虑在线到达方向估计以及对随时间而变化的到达方向进行跟踪的问题。

为了解决自适应到达方向估计和对随时间而变化的到达方向进行跟踪的问题，本发明人提出了采用高计算效率的SUMWE技术的被称为自适应方向估计和跟踪(ABEST)的自适应到达方向估计和跟踪方法(例如，参见J.Xin，Y.Ohashi，and A.Sano，“Efficient subspace-basedalgorithms for adaptive direction estimation and tracking ofnarrowband signals in array processing，”Proc.IFAC 8^th Workshopon Adaptation and Learning in Control and Signal Processing(ALCOSP’04)，pp.535-540，Yokohama，Japan，(2004))。

然而，上述用于跟踪无线电波的到达方向的方法无法精确地跟踪相干信号的交叉方向，其中入射信号的方向的轨迹由于信号源(例如呼叫方)的移动而彼此交叉。

发明内容

因此，本发明的目的是提供一种新的到达方向跟踪方法，用于以在线的方式估计无线电波的方向，而无需诸如特征分解的复杂处理，并且即使在无线电波的方向轨迹交叉的情况下，也能够快速跟踪方向。

这里，提出了一种新的用于阵列天线的无线电波到达方向跟踪方法，其中该阵列天线具有沿直线以相同的间隔设置在不同空间位置上的多个天线振子。该提出的方法包括以下步骤：以抽样间隔来计算某些天线振子所接收的阵列数据之间的瞬时相关性；根据该瞬时相关性来计算瞬时相关性矩阵；根据该瞬时相关性矩阵，通过线性运算来估计噪声子空间；使用观测器以方向更新间隔来预测方向的状态向量；根据所预测的方向和所估计的噪声子空间，以方向更新间隔来计算方向的临时估计值；以及根据该临时方向、所预测的状态向量，以及所预测的方向，以方向更新间隔来计算状态向量，从而根据该状态向量来计算该方向的估计值。

为了实现以上目的，根据本发明，提供了一种到达方向跟踪装置，用于估计到达阵列天线的无线电波的方向，该阵列天线具有沿直线以相同的间隔设置在不同空间位置上的多个天线振子。该无线电波方向跟踪装置包括以下单元：瞬时相关性计算单元，用于以抽样间隔来计算某些阵列数据之间的瞬时相关性；瞬时相关性矩阵计算单元，用于根据该瞬时相关性来计算瞬时相关性矩阵；噪声子空间估计单元，用于根据该瞬时相关性矩阵，通过线性运算来估计噪声子空间；方向预测单元，用于使用观测器以方向更新间隔来预测方向的状态向量；临时值计算单元，用于根据所预测的方向和所估计的噪声子空间，以方向更新间隔来计算方向的临时估计值；以及估计值计算单元，用于根据该临时方向、所预测的状态向量，以及所预测的方向，以方向更新间隔来计算状态向量，从而根据该状态向量来计算该方向的估计值。

为了实现以上目的，根据本发明，提供了一种基站装置，用于利用阵列天线来接收无线电波，该阵列天线具有沿直线以相同的间隔设置在不同空间位置上的多个天线振子。该基站装置包括以下部分：上述方向跟踪装置；以及波束成形部分，用于通过使用朝向所估计的到达方向形成波束来接收目标信号。

通过结合附图来阅读以下说明，本发明的以上和其他目的、特征和优点将变得明了，附图以示例的方式示出了本发明的优选实施例。

附图说明

图1是表示根据本发明实施例的方向跟踪方法的概要的流程图。

图2示出了方向更新间隔和抽样间隔之间的关系。

图3示出了阵列天线和发送源的结构。

图4示出了无线电波方向估计系统的结构。

图5是表示图4所示的无线电波方向跟踪部分的结构的方框图。

图6示出了非相关白噪声环境下的阵列协方差矩阵。

图7示出了瞬时阵列协方差矩阵中估计无线电波方向所需的列元素。

图8示出了瞬时阵列协方差矩阵中估计无线电波方向所需的行元素。

图9示出了通过瞬时阵列协方差矩阵的第一列或最后一列中的元素形成瞬时相关性矩阵，以及将该瞬时相关性矩阵分为上子矩阵和下子矩阵。

图10示出了通过瞬时阵列协方差矩阵的第一行或最后一行中的元素形成瞬时相关性矩阵，以及将该瞬时相关性矩阵分为上子矩阵和下子矩阵。

图11是表示对到达方向进行估计的仿真结果的曲线图。

图12示出了根据一实施例的基站接收装置的结构。

图13示出了根据一实施例的基站发送装置的结构。

图14示出了均匀线性阵列中的子阵列。

具体实施方式

下面将参照附图来详细说明本发明的实施例。

下面概要描述根据本发明实施例的无线电波到达方向跟踪方法。

根据本发明实施例的无线电波到达方向跟踪方法尤其适用于下述的情况：由于单个源(例如呼叫方)的移动而导致相干信号的方向轨迹彼此交叉。

当前的描述假定构成阵列天线的M个天线振子接收p个到达的无线电波(其中M＞2p)。本发明适用于前向子阵列、后向子阵列，以及前向和后向子阵列。

此外，本发明适用于对包括非相关信号、相关信号，以及相干信号(完全相关)信号在内的入射信号的跟踪。

图1是表示根据本发明实施例的方向跟踪方法的概要的流程图。

图2示出了方向更新间隔和抽样间隔之间的关系。

假设与抽样间隔TS相比，方向随时间的变化要更慢，并且方向更新间隔(即，对到达方向的计算进行更新的间隔)T可以表示为T＝NT_S，其中N是周期T期间的快照的数量。另外，还假设抽样时刻k和方向更新时刻n之间的关系可以表示为k＝nN，nN+1，…，(n+1)N-1。

根据图1所示的实施例的方向跟踪方法的步骤如下所述。首先，在抽样时刻k，使用方向更新间隔T期间的N个快照来计算某些天线振子的阵列数据之间的瞬时相关性(步骤S1)。

然后，在抽样时刻k(其中抽样速率为1/T_S)，根据该瞬时相关性来形成1至4个Hankel相关性矩阵(步骤S2)，这些相关性矩阵是瞬时相关性矩阵。

此后，根据在步骤S2中形成的瞬时相关性矩阵，通过使用自适应算法(例如，最小均方(LMS)方法，或者采用固定或随时间变化的步长的归一化最小均方(NLMS)方法)进行线性运算，估计在抽样时刻k的噪声子空间(步骤S3)。

此外，通过使用前一个方向更新时刻n-1的方向状态向量(每一个都包括到达方向和该方向的角速度和加速度)，利用观测器来预测当前方向更新时刻n的方向状态向量(步骤S4)。

通过诸如近似牛顿方法的自适应方法，根据所预测的方向来计算当前方向更新时刻n的信号方向的临时估计值(步骤S5)，其中所预测的方向是根据所预测的状态向量和在抽样时刻k＝(n+1)N-1获得的噪声子空间而获得的。

此外，根据该临时估计值、所预测的状态向量，以及所预测的方向，来计算当前方向更新时刻n的状态向量，然后根据所计算的状态向量计算该方向的估计值(步骤S6)。

通过根据该实施例的上述方向跟踪方法，可以精确地跟踪方向，即使在相干信号的方向轨迹由于单个源(例如呼叫方)的移动而彼此交叉的情况下，也不会增加计算负荷。

可以在任意时刻执行步骤S4中的处理，只要在步骤S5中的处理之前执行即可。

以下详细说明本实施例。

图3示出了发送源和阵列天线的结构。

直达波11是从发送源10直接到达基站的阵列天线20的无线电波。反射波12在经过诸如建筑物BL1和BL2的物体处反射后到达基站。尽管图3中仅示出了两个反射波12，但是本说明书假设从发送源10发射出包括直达波11和反射波12在内的总共p个波。另外还假设入射信号的数量p是已知的。该直达波11和反射波12之间在特定抽样时刻k的关系如公式(10)所示。

s_i(k)＝β_is₁(k)    …(10)

其中β_i表示多径系数，其表示反射波s_i(k)相对于直达波s_l(k)的复衰减。这里，β_i≠0并且β₁＝1。

图4示出了方向估计系统的结构。

该方向估计系统包括：阵列天线20；基带和数字处理部分30，用于执行基带处理和数字信号处理；以及方向跟踪部分40，用于执行根据本实施例的方向跟踪处理。

阵列天线20包括M个天线振子21。这里，需要满足条件M＞2p，其中p是到达阵列天线20的无线电波(直达波加上反射波)的数量。

图5是表示方向跟踪部分40的结构的方框图。

该方向跟踪部分40包括：计算单元41，用于计算瞬时相关性；计算单元42，用于计算瞬时相关性矩阵；计算单元43，用于计算线性算子(operator)；计算单元44，用于计算正交投影算子；计算单元45，用于计算方向的临时估计值；预测单元46，用于利用观测器预测方向；以及计算单元47，用于利用观测器计算方向。

下面将详细说明该方向跟踪。

首先，如公式(11)所示，用于利用观测器来预测方向的预测单元46通过使用前一个方向更新时刻(n-1)的方向状态向量利用观测器来预测在方向更新时刻n的方向状态向量。

$>>>>x>^>>i>>>(>n>|>n>->1>)>>=>F>>>x>^>>i>>>(>n>->1>|>n>->1>)>>,>>>\theta >^>>i>>>(>n>|>n>->1>)>>=>>c>T>>>>x>^>>i>>>(>n>|>n>->1>)>>->->->>(>11>)>>>>$

其中 $>>>x>i>>>(>n>)>>>=>\Delta >>[>>\theta >i>>>(>n>)>>,>>>\theta >·>>i>>>(>n>)>>,>>>\theta >>·>·>>>i>>>(>n>)>>>]>T>>>>表示方向\theta$_i(n)的状态向量，和是方向θ_i(n)的角速度和加速度。F和c分别是转置矩阵和观测向量，以下分别对它们进行定义。

$>>F>>=>\Delta >\left(\begin{array}{}>>>1>>>T>>>>>0.5>T>>2>>>>>>0>>>1>>>T>>>>>0>>>0>>>1>>>\end{array}\right)>,>c>>=>\Delta >>>>[>1,0,0>]>>T>>>>$

该处理与图1的步骤S4中的处理相对应。

接下来，根据天线振子所接收的数据，用于计算瞬时相关性的计算单元41根据从基带和数字处理部分30获得的复数字信号的N个快照{y₁(k)，y₂(k)，…，y_M(k)}_k＝nN^(n+1)N-1，如公式(1)所示生成接收信号向量y(k)。此外，根据公式(12a)和(12b)生成(k)和(k)。该(k)和(k)为信号y(k)和y*_M(k)以及y(k)和y*₁(k)之间在抽样时刻k的相关性向量。

通常，当根据阵列天线20所接收的信号来估计方向时，使用阵列协方差矩阵R，该矩阵的元素是天线振子21的信号向量y(k)(＝y₁(k)，y₂(k)，…，y_M(k))中的所接收信号之间的计算出的相关性r₁₁至r_MM。该阵列协方差矩阵R在图6所示的非相关白噪声环境下给出，其中y^H(k)表示所接收信号向量y(k)的复共轭转置。

图6示出了非相关白噪声环境下的阵列协方差矩阵R。

此外，以下等式成立：

y_i(k)＝x_i(k)+w_i(k)，

$>>E>[>>w>i>>>(>k>)>sup>>w>j>*sup>>>(>k>)>>]>=>>\sigma >2>>>(>i>=>j>)>>,>>>$

$>>E>[>>w>i>>>(>k>)>sup>>w>j>*sup>>>(>k>)>>]>=>0>>(>i>\ne >j>)>>>>$

其中x_i(k)表示无噪声的接收信号，而w_j(k)表示非相关白噪声。简而言之，在非相关白噪声环境下，该阵列协方差矩阵R的对角元素r₁₁、r₂₂、…、r_MM中含有噪声。

图7示出了在瞬时阵列协方差矩阵R(k)中估计方向所需的列元素。

图8示出了在瞬时阵列协方差矩阵R(k)中估计方向所需的行元素。

为了估计到达方向，由于瞬时阵列协方差矩阵R(k)具有共轭对称性，所以对第一列和最后一列(即图7所示的第M列)，或者对第一行和最后一行(即图8所示的第M行)进行计算就足够了。更具体地，仅根据公式(12a)计算第M天线振子所接收的数据与第一、第二、…第(M-1)天线振子所接收的数据之间的相关性，或者根据公式(12b)计算第一天线振子所接收的数据与第二、第三、…第M天线振子所接收的数据之间的相关性就足够了。以下说明使用公式(12a)和(12b)两者的情况。

再次参照图5，用于计算瞬时相关性矩阵的计算单元42使用通过公式(12a)和(12b)获得的相关性向量来计算(M-p)×p的Hankel相关性矩阵，如公式(13)所示，该矩阵是在抽样时刻k的瞬时相关性矩阵。

Φ_f(k)＝Hank{h_c，h_r}，Φ_f(k)＝Hank{h_c，h_r}

$>>>\Phi >b>>>(>k>)>>=>>J>>M>->p>>sup>>\Phi >f>*sup>>>(>k>)>>>J>p>>,>>>\Phi >‾>>b>>>(>k>)>>=>>J>>M>->p>>sup>>\Phi >f>*sup>>>(>k>)>>>J>p>>>>\dots (13)$

其中

$>>>h>c>>=>[>>>r>^>>>1>M>>>>(>k>)>>,>>>r>^>>>2>M>>>>(>k>)>>,>.>.>.>,>>>r>^>>>M>->p>,>M>>>>(>k>)>>>]>T>>,>>h>r>>=>[>>>r>^>>>M>->p>,>M>>>>(>k>)>>,>>>r>^>>>M>->p>+>1>,>M>>>>(>k>)>>,>.>.>.>,>>>r>^>>>M>->1>,>M>>>>(>k>)>>>]>T>>,>>>$

$>>>>h>‾>>c>>=>[>>>r>^>>21>>>(>k>)>>,>>>r>^>>31>>>(>k>)>>,>.>.>.>,>>>r>^>>>L>1>>>>(>k>)>>>]>T>>,>>>h>‾>>r>>=>[>>>r>^>>>L>1>>>>(>k>)>>,>>>r>^>>>L>+>1,1>>>>(>k>)>>,>.>.>.>,>>>r>^>>>M>1>>>>(>k>)>>>]>T>>,>>>L＝M-p+1，并且$

J_m是m×m的逆矩阵。

此外，用于计算瞬时相关性矩阵的计算单元42将各个(M-p)×p的Hankel相关性矩阵都分为上部和下部，如公式(14)所示。

…(14)

下面将更加详细地说明公式(12a)、(12b)、(13)以及(14)。

用于计算瞬时相关性的计算单元41通过执行公式(12a)的计算来获得在抽样时刻k阵列天线20的第M天线振子21所接收的数据与第一、第二、…、第(M-1)天线振子所接收的数据之间的瞬时相关性r_1M(k)、r_2M(k)、…r_M-1，M(k)(参见图7的最后一列)。用于计算瞬时相关性矩阵的计算单元42通过执行如图9所示的公式(13)的计算，来获得Hankel相关性矩阵Φ_f(k)，然后将该瞬时相关性矩阵分为上部和下部。

图9表示通过瞬时阵列协方差矩阵R(k)的第一列或最后一列中的元素形成瞬时相关性矩阵，以及将该瞬时相关性矩阵分为上子矩阵和下子矩阵。

用于计算瞬时相关性矩阵的计算单元42从图9所示的计算出的瞬时阵列协方差矩阵R(k)中的最后一列中的(M-1)个瞬时相关性当中选取(M-p)组瞬时相关性(一次一组并同时将选取位置向下移动一个瞬时相关性)，每一组都包括与信号数量(＝p)相同数量的瞬时相关性。计算单元42随后将所选取的瞬时相关性排列成多行，以从第一行开始生成(M-p)×p的瞬时相关性矩阵Φ_f(k)。此后，计算单元42将Φ_f(k)分为公式(14)所示的两个子矩阵：上部子矩阵是p×p矩阵Φ_f1(k)，而下部子矩阵是(M-2p)×p矩阵Φ_f2(k)。

类似地，用于计算瞬时相关性的计算单元41通过执行公式(12b)的计算来获得在抽样时刻k阵列天线20的第一天线振子21所接收的数据；以及第二、第三、…、第M天线振子所接收的数据之间的瞬时相关性r₂₁(k)、r₃₁(k)、…r_M1(k)(参见图7的第一列)。用于计算瞬时相关性矩阵的计算单元42执行公式(13)的计算以从图9所示的瞬时阵列协方差矩阵R(k)的第一列中的(M-1)个瞬时相关性当中选取(M-p)组瞬时相关性一次一组并同时将选取位置向下移动一个瞬时相关性，每一组都包括与信号数量(＝p)相同数量的瞬时相关性。计算单元42随后将所选取的瞬时相关性排列成多行，以从第一行开始生成(M-p)×p的瞬时相关性矩阵Φ_f(k)。此后，计算单元42将Φ_f(k)分为公式(14)所示的两个子矩阵：上部子矩阵是p×p矩阵Φ_f1(k)，而下部子矩阵是(M-2p)×p矩阵Φ_f2(k)。

还对瞬时阵列协方差矩阵R(k)的多行执行上述处理。更具体地，用于计算瞬时相关性的计算单元41通过执行公式(12b)的计算来获得在抽样时刻k阵列天线20的第一天线振子21所接收的数据与第二、第三、…、第M天线振子所接收的数据之间的瞬时相关性r₁₂(k)、r₁₃(k)、…r_1M(k)(参见图8的第一行)。用于计算瞬时相关性矩阵的计算单元42通过执行公式(13)的计算来获得Hankel相关性矩阵(瞬时相关性矩阵)Φ_b(k)，如图10所示，然后将该瞬时相关性矩阵分为上部和下部。

图10表示通过瞬时阵列协方差矩阵R(k)的第一行或最后一行中的元素形成瞬时相关性矩阵，以及将该瞬时相关性矩阵分为上子矩阵和下子矩阵。

用于计算瞬时相关性矩阵的计算单元42从图10所示的瞬时阵列协方差矩阵R(k)的第一行中的(M-1)个瞬时相关性当中选取(M-p)组瞬时相关性(一次一组并同时将选取位置从右向左移动一个瞬时相关性)，每一组都包括与信号数量(＝p)相同数量的瞬时相关性。计算单元42随后将所选取的瞬时相关性排列成多行，以从第一行开始生成(M-p)×p的瞬时相关性矩阵Φ_b(k)。此后，计算单元42将Φ_b(k)分为公式(14)所示的两个子矩阵：上部子矩阵是p×p矩阵Φ_b1(k)，而下部子矩阵是(M-2p)×p矩阵Φ_b2(k)。

类似地，用于计算瞬时相关性的计算单元41通过执行公式(12a)的计算来获得在抽样时刻k阵列天线20的第一天线振子21所接收的数据与第一、第二、…、第(M-1)天线振子所接收的数据之间的瞬时相关性r_M1(k)、r_M1(k)、…r_M，M-1(k)(参见图8的最后一行)。用于计算瞬时相关性矩阵的计算单元42随后执行公式(13)的计算，以从图10所示的瞬时阵列协方差矩阵R(k)的最后一行中的(M-1)个瞬时相关性当中选取(M-p)组瞬时相关性(一次一组并同时将选取位置从右向左移动一个瞬时相关性)，每一组都包括与信号数量(＝p)相同数量的瞬时相关性。计算单元42随后将所选取的瞬时相关性排列成多行，以从第一行开始生成(M-p)×p的瞬时相关性矩阵Φ_b(k)。此后，计算单元42将Φ_b(k)分为公式(14)所示的两个子矩阵：上部子矩阵是p×p矩阵Φ_b1(k)，而下部子矩阵是(M-2p)×p矩阵Φ_b2(k)。

通过这种方式，生成了四个Hankel相关性矩阵及其八个分割子矩阵。

下面将说明估计噪声子空间的处理。用于计算线性算子的计算单元43使用公式(14)中所示的八个矩阵来生成表示为公式(15)的两个矩阵Φ₁(k)和Φ₂(k)，以计算表示为公式(16)的估计误差矩阵E(k)。

$>>>\Phi >1>>>(>k>)>>>=>\Delta >>[>>\Phi >>f>1>>>>(>k>)>>,>>>\Phi >‾>>>f>1>>>>(>k>)>>,>>\Phi >>b>1>>>>(>k>)>>,>>>\Phi >‾>>>b>1>>>>(>k>)>>]>>>$

$>>>\Phi >2>>>(>k>)>>>=>\Delta >>[>>\Phi >>f>2>>>>(>k>)>>,>>>\Phi >‾>>>f>2>>>>(>k>)>>,>>\Phi >>b>2>>>>(>k>)>>,>>>\Phi >‾>>>b>2>>>>(>k>)>>]>>> \dots (15)$

$>>E>>(>k>)>>>=>\Delta >sup>>\Phi >2>Hsup>>>(>k>)>>-sup>>\Phi >1>Hsup>>>(>k>)>>P>>(>k>->1>)>>->->->>(>16>)>>>>$

用于计算线性算子的计算单元43随后使用LMS算法来获得由公式(17)给出的线性算子P(k)。

p(k)＝p(k-1)+μΦ₁(k)E(k)        …(17)

其中μ表示步长，并且是满足条件表达式0＜μ＜2/tr{Ψ₁}的正值。这里， $>>>>\psi >‾>>1>>>=>\Delta >>\{>>\Phi >1>>>(>k>)>sup>>\Phi >1>Hsup>>>(>k>)>>\}>>>和tr\{·\}是轨迹运算。$

接下来，用于计算正交投影算子的计算单元44通过使用从公式(17)获得的P(k)来执行如公式(18)所示的QR分解。

$>>P>>(>k>)>>>P>H>>>(>k>)>>+>>I>P>>>=>\Delta >>>p>‾>>=>>QR>‾>>->->->>(>18>)>>>>$

此外，通过公式(19)获得正交投影算子П(k)，以估计噪声子空间。

П(k)＝Q(k)(I_M-2p-P^H(k)R^-1Q^HP(k))Q^H(K)    …(19)

接下来，用于计算该方向的临时估计值的计算单元45使用从公式(11)获得的预测方向以及从公式(19)获得的П(n)＝П(k)|_k＝(n+1)N-1，根据公式(20)使用近似牛顿方法来计算抽样时刻k＝(n+1)N-1(即，方向更新时刻n)处的到达方向的临时估计值

$>>>>\theta >~>>i>>>(>n>)>>=>>>\theta >^>>i>>>(>n>|>n>->1>)>>->>>Re>\{>>>d>‾>>H>>>(>\theta >)>>\Pi >>(>n>)>>>a>‾>>>(>\theta >)>>\}>>>>>d>‾>>H>>>(>\theta >)>>\Pi >>(>n>)>>>d>‾>>>(>\theta >)>>>>>|>>\theta >=>>>\theta >~>>i>>>(>n>|>n>->1>)>>>>->->->>(>20>)>>>>$

其中

$>>>d>‾>>>(>\theta >)>>=>j>>\omega >0>>>(>d>/>c>)>>cos>>\theta >k>>[>0>,>>e>>j>>\omega >0>>\tau >>(>\theta >)>>>>,>2>>e>>j>2>>\omega >0>>\tau >>(>\theta >)>>>>,>.>.>.>,>>(>L>->2>)>>>e>>j>>\omega >0>>>(>L>->2>)>>\tau >>(>\theta >)>>>>>]>T>>>>$

最后，用于利用观测器来计算到达方向的计算单元47根据公式(21)，利用观测器来计算在方向更新时刻n的到达方向的状态向量，以通过使用从公式(20)获得的该方向的临时估计值以及从公式(11)获得的和给出该到达方向的估计值

$>>>>x>^>>i>>>(>n>|>n>)>>=>>>x>^>>i>>>(>n>|>n>->1>)>>+>>g>i>>>(>>>\theta >~>>i>>>(>n>)>>->>>\theta >^>>i>>>(>n>|>n>->1>)>>)>>,>>>\theta >^>>i>>>(>n>|>n>)>>=>>c>T>>>>x>^>>i>>>(>n>|>n>)>>->->->>(>21>)>>>>$

其中g_i表示观测器增益。应该注意，观测器增益g_i被设定为使得F-g_iC^T的所有特征值都位于单位圆中。

如上所述，方向跟踪部分40可以通过使用某一固定步长μ的LMS算法以及近似牛顿方法来在线跟踪信号的到达方向。另外，可以精确地估计出方向，而即使在直达波和反射波的到达方向的轨迹由于单个源(例如呼叫方)的移动而彼此交叉的情况下，也不会增大计算负荷。

在上述实施例中，获得了四个Hankel相关性矩阵Φ_f(k)、Φ_f(k)、Φ_b(k)、Φ_b(k)，并且每一个Hankel相关性矩阵都被分为两个子矩阵Φ_f1(k)、Φ_f2(k)；Φ_f1(k)、Φ_f2(k)；Φ_b1(k)、Φ_b2(k)；Φ_b1(k)、Φ_b2(k)，以确定两个矩阵Φ₁(k)和Φ₂(k)。然而，为了估计该到达方向，由于阵列协方差矩阵R具有共轭对称性，所以使用图7所示的第一列或最后一列中的瞬时相关性、图8所示的第一行或最后一行中的瞬时相关性，或者任意两个或更多行或列中的瞬时相关性就足够了。简言之，使用四个Hankel相关性矩阵中的任意1至4个就足以确定矩阵Φ₁(k)和Φ₂(k)。

例如，如果要使用四个Hankel相关性矩阵中的任意一个，则根据公式(22a)至(22d)中的任意一个来确定矩阵Φ₁(k)和Φ₂(k)。

Φ₁(k)＝Φ_f1(k)，Φ₂(k)＝Φ_f2(k)                                                  …(22a)

Φ₁(k)＝Φ_f1(k)，Φ₂(k)＝Φ_f2(k)                                                …(22b)

Φ₁(k)＝Φ_b1(k)，Φ₂(k)＝Φ_b2(k)                                                  …(22c)

Φ₁(k)＝Φ_b1(k)，Φ₂(k)＝Φ_b2(k)                                                …(22d)

如果要使用四个Hankel相关性矩阵中的任意两个，则根据公式(23a)至(23f)中的任意一个来确定矩阵Φ₁(k)和Φ₂(k)。

Φ₁(k)＝[Φ_f1(k)，Φ_f1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_f2(k)]                        …(23a)

Φ₁(k)＝[Φ_f1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_b2(k)]                          …(23b)

Φ₁(k)＝[Φ_f1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_b2(k)]                        …(23c)

Φ₁(k)＝[Φ_f1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_b2(k)]                        …(23d)

Φ₁(k)＝[Φ_f1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_b2(k)]                      …(23e)

Φ₁(k)＝[Φ_b1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_b2(k)，Φ_b2(k)]                        …(23f)

如果要使用四个Hankel相关性矩阵中的任意三个，则根据公式(24a)至(24d)中的任意一个来确定矩阵Φ₁(k)和Φ₂(k)。

Φ₁(k)＝[Φ_f1(k)，Φ_f1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_f2(k)，Φ_b2(k)]    …(24a)

Φ₁(k)＝[Φ_f1(k)，Φ_f1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_f2(k)，Φ_b2(k)]  …(24b)

Φ₁(k)＝[Φ_f1(k)，Φ_b1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_b2(k)，Φ_b2(k)]    …(24c)

Φ₁(k)＝[Φ_f1(k)，Φ_b1(k)，Φ_b1(k)]，Φ₂(k)＝[Φ_f2(k)，Φ_b2(k)，Φ_b2(k)]  …(24d)

已对使用四个Hankel相关性矩阵的情况进行了说明。

另外，尽管在上述实施例中，线性算子P(k)是通过具有特定的固定步长μ的LMS算法从公式(17)获得的，但是该步长μ可以随时间而变化。更具体地，通过根据公式(25)来确定公式(17)中的步长μ，以使步长μ可以在抽样时刻k随瞬时相关性矩阵Φ₁(k)而变化，可以实现步长μ随时间而变化的LMS算法。

$>>\mu >=>2>/>tr>\{>>\Phi >1>>>(>k>)>sup>>\Phi >1>Hsup>>>(>k>)>>\}>->->->>(>25>)>>>>$

此外，可以通过表示为公式(26)的NLMS算法，而不是LMS算法，来计算该线性算子P(k)。

$>>P>>(>k>)>>=>P>>(>k>->1>)>>+>>\mu >‾>>>>R>~>>>->1>>>>>Q>~>>H>>>\Phi >1>>>(>k>)>>E>>(>k>)>>->->->>(>26>)>>>>$

其中和是通过 $>>>\Phi >1>>>(>k>)>sup>>\Phi >1>Hsup>>>(>k>)>>=>>Q>~>>>R>~>>>>获得的，其表示\Phi$₁(k)Φ₁^H(k)的QR分解。这里，NLMS技术的步长μ应该满足稳定条件0＜μ＜2。

下面说明根据本实施例的方向跟踪方法的仿真结果。

图11表示对到达方向进行估计的仿真结果的曲线图。

在图11的曲线图中，纵轴与角度相对应，而横轴与时间刻度相对应。在该曲线图中，虚线表示所接收信号的实际到达方向，而实线表示估计值的曲线。

均匀线性阵列(ULA)中的元素数量为M＝9。四个相干信号从SNR为15dB、10dB、13dB和13dB的角度θ₁(n)、θ₂(n)、θ₃(n)和θ₄(n)入射到阵列天线上。在方向更新期间间隔T＝1s期间观测到了N(＝100)个快照，并且通过公式(26)中的NLMS算法来计算在方向更新时刻n的正交投影算子П(n)＝П(k)|_k＝(n+1)N-1，以从公式(20)获得到达方向的估计值在该示例中，NLMS方法的步长为μ＝0.88。对通过100次算术运算而得到的估计值进行绘图，如图11所示，该图表示，即使入射信号的到达方向的轨迹彼此交叉，这些估计值也与实际的到达方向精确地一致。如上所述，根据本实施例的方向跟踪方法可以快速并精确地估计随时间而变化的、完全相关的信号的到达方向，而即使这些到达方向彼此交叉，也无需复杂的特征分解。

尽管在上述实施例中，对非相关白噪声环境下的方向跟踪方法进行了说明，但是本发明还可以应用于空间相关的噪声环境。然而，在空间相关的噪声环境中，需要执行以下过程以生成瞬时相关性矩阵。

在空间相关的噪声环境下，需要根据公式(27)来生成Hankel相关性矩阵，而不是根据公式(13)来获得(M-p)×p的Hankel相关性矩阵，在公式(27)中，q表示噪声的空间相关性长度(即，如果|k|＞q，则E{w_i(n)w^*_i+k(n)}＝0)。

更具体地，在空间相关的噪声环境中，可以根据在抽样时刻k阵列天线20中的第M天线振子21所接收到的信号与第一、第二、…、第(M-p-q)天线振子所接收到的信号之间的相关性，以及在抽样时刻k阵列天线20中的第一天线振子21所接收到的信号与第(q+2)、(q+3)、…、第M天线振子所接收到的信号之间的相关性，以与图9和图10中相同的方式生成如方程(27)所示的在抽样时刻k的(M-p-q)×p的Hankel相关性矩阵。

如果M＞2p+q，则可以将各个Hankel相关性矩阵都分为如公式(27)右手侧而非公式(14)所示的两个子矩阵：p×p矩阵和(M-2p-q)×p矩阵。此后，根据上述矩阵，可以按照与非相关白噪声环境相同的方式，通过从公式(15)生成矩阵Φ₁(k)和Φ₂(k)，来估计方向。

可以通过组合用于跟踪方向的上述方向跟踪部分40和波束成形单元来构成基站接收装置，该波束成形单元用于生成接收波束模式，以使得峰值朝向所估计的到达方向。

图12示出了根据本实施例的基站接收装置的结构。

该基站接收装置例如设置在基站中，并且包括阵列天线20、基带和数字处理部分30、方向跟踪部分40、瞬时波束成形器50，以及信道接收部分60。

下面简要说明该基站接收装置的操作。

当阵列天线20接收到信号时，基带和数字处理部分30对每一个天线振子执行信号处理，以输出复数字接收数据。方向跟踪部分40接收该复数字接收数据，并执行上述跟踪方向的处理，以计算在方向更新时刻n的到达方向的估计值。瞬时波束成形器(接收波束成形器)50根据所估计的到达方向来形成波束，以使得峰值朝向信号源。简言之，瞬时波束成形器50在抑制干扰和噪声的同时提取目标信号，以将该目标信号发送至信道接收部分60。信道接收部分60利用公知的方法来执行接收处理，以对接收到的数据进行解调和输出。使用与根据本实施例的上述方向跟踪方法获得的到达方向相关的信息的瞬时波束成形器50可以有多种结构。为了接收信号，可以利用例如以下文献中所公开的波束成形方法使波束朝向目标信号到达方向：O.L.Frost，“An algorithm for linearlyconstrained adaptive array processing，”Proc.IEEE，vol.60，no.8，pp.926-935(1975)；以及J.Xin，H.Tsuji，Y.Hase，and A.Sano，“Array beamforming based on cyclic signal detection，”Proc.IEEE48^th Vehicular Technology Conference，pp.890-894，Ottawa，Canada(1998)。

可以通过组合用于跟踪方向的上述方向跟踪部分40和波束成形单元来构成基站发送装置，该波束成形单元用于生成发送波束模式，以使得峰值朝向所估计的到达方向。

图13示出了根据该实施例的基站发送装置的结构。图13包括图12中所示的基站接收装置。

当从发送部分80输入发送数据后，瞬时波束成形器(发送波束成形器)70形成发送波束模式，使得峰值朝向由方向跟踪部分40估计的到达方向，以向基带和数字处理部分30发送复数字发送信号。基带和数字处理部分30将该复数字发送数据转换为无线电信号，并将其发送至阵列天线20a中的各个天线振子21a。结果，向接收站发射了波束，以降低误码率。图13中的阵列天线20和20a可以集成为单个阵列天线。

根据本发明，通过使用前一个方向更新时刻的到达方向的状态向量，利用观测器来预测当前方向更新时刻的到达方向的状态向量，以计算该到达方向的临时估计值，此外，根据该临时估计值、所预测的状态向量，以及所预测的到达方向来计算当前方向更新时刻的状态向量，以根据该状态向量估计该到达方向。因此，可以计算出该方向的估计值，由此可以精确地实时跟踪该方向，而即使在直达波和反射波的到达方向的轨迹由于单个源(例如呼叫方)的移动而彼此交叉的情况下，也不会增加计算负荷。

前述内容仅是对本发明原理的说明。此外，由于本领域的技术人员很容易想到多种修改和变化，所以并不旨在将本发明限定为所示出和所说明的确切结构和应用，因此可以认为，所有适当的修改和等价物均落入所附权利要求及其等价物所限定的本发明的范围之内。